SEM²: Introducing mechanics in cell and tissue modeling using coarse-grained homogeneous particle dynamics

Morphogenesis, the process by which tissues and organs change shape to acquire novel structures and functions during embryonic development, has long fascinated life scientists^1,2 and physicists.³ More recently, engineers became interested in how tissues build themselves to design better preclinical models, including organoids^1,4,5 and organs-on-chips.^6–8 New tools in molecular biology^9,10 and optical engineering^11–13 have shown us that mechanical forces across spatial scales modulate morphogenesis.¹⁴ Yet, we have a limited ability to computationally model these interactions.¹⁵ 

From a modeling perspective, morphogenetic events are characterized by complex tissue rearrangements driven by the proliferation and migration of millions of cells.^2,9,16 In other words, we must model large non-linear deformations driven by discrete, local events in a continuously growing system. To model these events, continuous approaches treat developing tissues as deforming solids, flowing liquids and/or visco-poro-elastic combinations.¹⁷ Alternatively, cells can be modeled as discrete objects interacting on a lattice,^15,18,19 or freely in space.^3,20 All modalities face tradeoffs. We can simulate cell mechanics with continuous models, but it is hard to incorporate discrete biology. Discrete models more naturally capture biological processes but are costly, and cell mechanics is lost (e.g., Potts models¹⁹) or limited to cell membranes (e.g., vertex modeling³). Yet, the mechanical forces at play in the intracellular space, or cytoplasm, can influence the mechanical behavior of entire tissues.^21,22 To assess mechanics during morphogenetic events, we propose an extension to the class of methods known as subcellular element modeling (SEM), which we termed subcellular element modeling and mechanics or SEM².

SEM is a discrete, off-lattice framework that treats cells as mechanically competent agents^23,24 capable of biologically relevant behaviors, such as flowing,²⁵ migration,²⁶ proliferation,^27,28 and epithelia shaping.^29,30 In SEM, each cell is formed by N_p particles that provide a coarse-grained representation of the subcellular space [Fig. 1(a)]. Two key results make SEM relevant for cell and tissue mechanics. First, Newman established the link with cell mechanics using experimentally measured cell rheology to define from first principles a family of mildly attractive and mildly repulsive potentials that enforce adhesion and volume exclusion [Fig. 1(b)].^23,24 Later, Koumoutsakos introduced SEM++,²⁷ a C++ implementation of SEM that leverages the open-source and strongly supported Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) library³¹ to speed up 3D computations. Furthermore, algorithmic manipulations of the particles have been used to model cell migration and proliferation,²⁷ while multiple particle types were introduced to model cell nuclei or membranes, explicitly.^27,32 Together these results have made SEM the framework of choice for modeling cell shape as an emergent property of subcellular interactions.³³ Yet, a treatment of mechanics across spatial scales in SEM is missing.

In SEM², we realized that known potentials and computed particle displacement could be used to estimate sub-cellular mechanics via per-particle stress and strain and relate them to cell and tissue level deformation during morphogenesis. In this paper, we updated SEM++ and combined it with new analyses and visualizations to (i) simulate particle ensembles [Fig. 1(c)] during key cellular behaviors, such as cell division [Fig. 1(d)]; and, (ii) study how stress and strain propagate across spatial scales. To demonstrate this approach, we present in silico studies of cells undergoing creep, migration, and proliferation in traditional and engineered cell culture environments. We chose these simulation scenarios because they have been performed and analyzed for cell-level behaviors in the SEM literature,^23,27 which enabled us to validate SEM² against the previous results and demonstrate SEM² capabilities for more detailed analysis of multiscale mechanics. Then, we used SEM² to simulate cell proliferation in unconstrained and constrained conditions to showcase how SEM² could be used for designing preclinical models such as organoids and microfluidic organs-on-chips, respectively. SEM² and the input scripts for all presented simulations are available on our GitHub at https://github.com/Synthetic-Physiology-Lab/sem2.

In SEM², we wanted to verify that these relations hold for arbitrary $N p$⁠. To demonstrate it, we resorted to modeling a classical SEM experiment, cell creep. In this setup, cells are loaded between two glass pipettes attached to force transducers.^34,35 As one pipette is kept stationary, the other is displaced, thus applying controlled stress on the cell and generating measurable strain via cell lengthening [Fig. 2(a)]. To simulate this experiment, we set up a simulation with a constrained cell composed of $N p$ particles, with $N p$ ranging from 250 to 10 000 (all other parameters are listed in Table I). To repeat these simulations as a function of particle number, one must convert the constant stress applied via the movable pipette in all simulations to the appropriate external force experienced by the particles in contact with the glass pipettes. However, the surface area in contact with the glass pipette and the precise number of particles in that area vary non-linearly with N_p based on cell geometry and particle stacking (supplementary material Fig. S1). To solve this problem, we first created two parallel, flat surfaces, known as slabs in LAMMPS [black rectangle, Fig. 2(b)].^23,24 Then, we calculated the surface area [red shaded area, Fig. 2(b)] of the slabs and used it to convert the applied stress to the total force that we allocated to the particles in the top slab (see Methods).

Having standardized force application, we asked whether the relationships in Eqs. (3a)–(3d) could be regarded as scaling laws: that is, do the cells composed of 250 and 10 000 particles undergo the same strain, under the same applied stress? To answer this question, we first considered a cell composed of a single particle type [no nucleus, Fig. 2(b)] and subjected it to a short-duration stress [5 Pa for 7 s, Fig. 2(a)]. To account for the effect of stochastic thermal fluctuations, we performed three different simulations per value of N_p and plotted the average (solid curve) and the standard deviation [shaded area, Fig. 2(c)] for each condition. Our results show that SEM² predicts a similar axial strain time course for all simulated $N p$ [Fig. 2(c)]. Second, following Milde et al.,²⁷ we introduced the cell nucleus as a more massive particle of greater radius [Fig. 2(d), supplementary material video SV1] which is accomplished by using a shifted version of the potential in Eq. (2) to accommodate the nuclear radius. Then, we repeated the simulations to demonstrate that the scalability of the approach can be extended to multiple particle types [Fig. 2(e), a detailed plot of axial and lateral strain can be seen in Fig. S4].

While small-strain regimes are typical of a healthy physiology,³⁶ we reasoned that high-strain situations occur in biology and are likely to produce more interesting subcellular mechanics. Therefore, we performed another round of in silico creep experiments for a fixed $N p$ = 1000 without [Fig. 3(a)] and with the nucleus [Fig. 3(b), supplementary material video SV2]. To explore the large strain regime, we let the cell equilibrate for 10 s, stepped the applied stress to 20 Pa, and held it for 30 s before allowing it to relax for another 30 s [blue dashed line, Figs. 3(c) and 3(d)]. We assessed the resulting cell strain [black solid line, Figs. 3(c) and 3(d)] and calculated the per-particle shear strain based on the reference configuration immediately before stress application [color-coded in Figs. 3(a) and 3(b)]. Both with and without the nucleus, the strain increased with the applied stress and decreased upon stress removal. Notably, and differently from the small strain regime, the cell strain did not reach a steady state and was ∼2.5 times smaller in the presence of the nucleus (note that the reference configuration for the cell with the nucleus was 22% longer than without it). To analyze inter-particle distance statistically, we calculated the radial distribution function (rdf) [Figs. 3(e) and 3(f)]. The rdf plot shows the relative abundance (Z-axis) of pairs of particles that are located at a given distance (X-axis) at various time points (Y-axis) during the simulation. In the equilibrium phase [dark green distributions, Figs. 3(e) and 3(f)], most particle pairs were near their equilibrium distance, d_eq ∼1.8 μm, as expected based on the rheology-preserving applied forces [Fig. 1(b) and Eqs. (1) and (2)].²⁷ Instead, many particles were strained away from their equilibrium locations during the stretch [light green distributions, Figs. 3(e) and 3(f)]. Finally, upon stress removal, particles could locally return to their equilibrium distance [gray distributions, Figs. 3(e) and 3(f)]. Yet, due to the overall different configuration, the total number of particle pairs at d_eq decreased, providing a subcellular rationale for the plastic deformation observed at the cell level [Figs. 3(c) and 3(d), Fig. S5 in the supplementary material]. These observations demonstrate that SEM² can be applied to simulate and analyze cell-level and particle-level mechanics, thus offering a novel way to study subcellular mechanics in classical cell creep experiments.

In Eq. (4), $F BIO y i , j , k$ is a general per-particle force and is intended to simulate accelerations imparted to certain particles by active biological mechanisms. For example, $F BIO$ can simulate reactive force being exchanged between the nucleus and cytoskeletal elements during translocation (Fig. S6) or a static force capable of directing nuclear movement toward the leading or trailing edge directly (Fig. S7). To balance accuracy and computational time (Fig. S6), we implemented a spring-like constraint. First, we introduced the variable $Δ n c$⁠, which measures the distance between the centroid of the cell and the centroid of the nucleus; then, we specified F_BIO as a force in a virtual spring between these centroids whose equilibrium length is zero [Fig. 4(a), Methods]. By adjusting the spring stiffness of this biasing mechanism, we could arbitrarily reduce $Δ n c$ to simulate various degrees of subcellular active transport ranging from it having no-bias [Fig. 4(a)], a weak bias [Fig. 4(b)], or a strong bias toward the cell centroid [Fig. 4(c), also supplementary material video SV4]. For these simulations, we chose a T_m that would result in nuclear drift in the original SEM++ framework (T_m =100 and T_sim = 100 s) and showed it being mitigated with different spring stiffnesses in SEM². These simulations clearly showed that during migration, the centroids of each cell and its nucleus (depicted as a black sphere and black cube, respectively) would drift apart if not counteracted. To investigate how much nuclear drift affects cell migration, we performed more in silico experiments altering T_m from 50 to 250 s while keeping the total duration to T_sim = 100 s. In the absence of F_BIO, a degree of nuclear drift was observed for all T_m [Fig. 4(d)], but the dependence of $Δ n c$ on T_m was reduced or completely negated by increasing the nuclear biasing strength [Fig. 4(e)]. Finally, the simulations showed that for the same choice of T_m and T_sim, cells moved a smaller distance in the absence of biasing forces [see trajectories in Fig. 4(f) and supplementary material video SV4], suggesting uncompensated nuclear drift affects cell motility.

To demonstrate how SEM^23,27 could be used to study the dynamics of specific subcellular elements, we considered two particles, one near and another away from the nucleus, in the simulations with and without biasing forces (Fig. 4). We observed that the two particles moved much farther apart in the no-bias condition than they did under strong bias [Fig. 5(a)]. In fact, the analysis of the full trajectories of the two particles during migration with no bias [Fig. 5(b)] and with a strong bias [Fig. 5(c)] indicated that the particle near the nucleus barely moved in the direction of migration.

Next, we asked whether the potential between cytoplasmic and nuclear particles might also impact subcellular mechanics inside the whole cell. Particle-level strain (Fig. 3) is ill-defined during migration simulations in SEM because particles are continuously added and removed from the ensemble, so there is no single reference particle configuration at the beginning of the experiment. Instead, we focused on per-particle mean volumetric stress (⁠ $σ m= σ x x + σ y y + σ z z 3$⁠), the relative extent of which can be assessed directly from the particle energy at any given time frame (see Methods). When we compared particles near the nucleus and away from it, we found that the magnitude and fluctuations of stress were higher in the particle near the nucleus for both the migration with no bias [Fig. 5(d)] and with a strong bias [Fig. 5(e)]. However, in the case of migration with a biasing force, the particle near the nucleus is mostly under greater stress [Fig. 5(e)] due to the push and pull of the nucleus, which is under a spring-like force toward the cell center.

We then asked whether this behavior was isolated to these two particles or existed across the cell. To answer this question, we color-coded the cytoplasmic particles of the migrating cells according to their stress normalized to the maximum stress computed during the simulation. Since the nucleus is a single rigid particle, we excluded it from this stress-dependent visualization (see Methods). We followed the convention that positive and negative stresses denote tensile and compressive behavior, respectively.³⁸ We noticed that the mean volumetric stresses were mostly tensile (positive) along the cell membrane, whereas the inside was under compressive (negative) stress. To analyze the distribution of subcellular stress, we computed the average per-particle hydrostatic stress in 1D bins and found tensile peaks in correspondence with the cell membrane. For the case of migration with no bias, the nucleus is at a peripheral position resulting in an asymmetric stress distribution [Fig. 5(f)]. However, for the case of migration with biasing force, since the nucleus is at the center of the cell, the particle stress distribution is roughly symmetric around the nucleus [Fig. 5(f)]. Taken together, these results suggest that SEM² can be used to model migration in more general cases than SEM++ and demonstrated how the per-particle analysis could be used to gain insights into subcellular mechanics during cell migration.

Finally, to tackle problems in tissue morphogenesis with SEM, we also need to model cell proliferation. To implement proliferation, SEM++ considers two separate phases: growth and division. In the growth phase, particles are added similarly to migration but without removal (see supplementary material video SV5). When the particle number per cell doubles to 2N_p, SEM++ enters the division mode based on the assumption that the cell grows to twice its normal size before dividing into two daughter cells. In division mode, a second nuclear particle is created and allowed to interact with the original nuclear particle via a mildly repulsive potential, imitating how the poles of the mitotic spindle drift apart during cell division.³⁹ To nucleate a daughter cell, the new nuclear particle is given a new cell index [the kth in the $y i , j , k$ notation introduced in Eq. (1)]. To complete cell division, the cell index of each cytoplasmic particle is re-assigned to one of the daughter cells based on the cell index of the closest nuclear particle, thus ensuring the splitting of the mother cell into two independent daughter cells. However, we preferred not to change the modeling assumptions regarding the nuclear particles to avoid artifacts due to particle size and weight changing during the simulation. Instead, SEM² uses the nuclear-centering bias introduced in the previous section (see Methods). First, it creates the second nuclear particle close to the first at the center of the cell, thus modeling the increased DNA density observed before cell division.⁴⁰ Then, it assigns cytoplasmic particles to the closer nucleus, thus splitting the cell at the center. Finally, the nuclear biasing mechanism recenters the nucleus of each daughter cell [Fig. 6(a), supplementary material video SV6].

To study subcellular mechanics during cell division with SEM,² we resorted again to per-particle calculations and visualization by computing the mean volumetric stress per particle. Figure 6(b) shows a proliferating cell during growth and division with the cytoplasmic particles color-coded according to their normalized stress (see Methods). We noticed that the volumetric stresses are mostly tensile (positive) at the cell–cell boundaries and compressive (negative) inside. This could be better understood when we created spatial bins in two [Fig. 6(c)] and one [Fig. 6(d)] dimensions. In 2D, we observe tensile stresses at the outer rim of particle ensembles and the cell–cell interface [Fig. 6(c)], whereas in 1D we see tensile (positive) peaks [Fig. 6(d) arrows indicate peaks]. Instead, particle volumetric stresses remained mostly compressive (negative) inside the cell [Figs. 6(c) and 6(d)]. In other words, the gradient between compressive and tensile stress could be used to pinpoint the position of the cell membrane. This is particularly interesting because, so far in SEM, information about membranes has been prescribed with dedicated membrane particles or obtained algorithmically.^27,32

To demonstrate that SEM² can be used for tissue engineering, we investigated subcellular mechanics during cell proliferation in simulation scenarios mimicking organoids^1,4,5 and organs-on-chips.^6–8 To do that, we performed three rounds of cell proliferation in unconstrained and constrained environments that model low-adhesion plates and the central section of a microfluidic channel, respectively. During the simulations, cells proliferated freely in the unconstrained environment [Fig. 7(a)] or along the channel in the constrained one [Fig. 7(b), supplementary material video SV7]. We obtained the particle stress distribution in these two mechanical environments and depicted them in Figs. 7(c) and 7(d), respectively (supplementary material video SV8). Again, particle stresses were tensile at the periphery and the cell–cell interfaces, whereas compressive inside. When we binned the particle stress along the unconstrained axis, the tensile peaks occurred when two or more cell–cell interfaces aligned, which was more frequent in the constrained environment [Fig. 7(d)].

Finally, these simulations allowed us to look at the problem from the perspective of multiscale mechanics. At the tissue level, the volumetric stress vs time of the two tissues was the same [Fig. 8(a)], since it was dominated by cell proliferation which is unaffected by mechanical constraints in this simulation. However, the axial strain was higher in the cylindrical channel [Fig. 8(b)], as expected from the mechanical environment that coerces newly formed cells to align axially (Fig. 7). Similarly, at the cell level, volumetric strain followed the cell cycle: it increased as the cell grew and decreased when it divided in both constrained and unconstrained simulations [Fig. 8(c)]. In other words, the growth of cells was not affected by the lateral constraint; instead, cells adapted their shape and orientation to make use of the unconstrained axial direction. This was confirmed by the particle-level stress distribution in 1D bins [Fig. 8(d)], where tensile peaks can be observed in correspondence of cell–cell interfaces, and they were higher and more frequent in the constrained condition [orange vs cyan lines, black arrows in Fig. 8(d)]. Together, these results suggest that SEM² can be used to gain insights into multiscale and subcellular mechanisms of cell proliferation in traditional and engineered cell culture platforms.

We presented SEM² as an updated version of SEM++ that enables the study of multiscale mechanics using coarse-grained particle dynamics. We extended the capability of SEM++ to model cell migration and proliferation²⁷ by adding a general force term, F_BIO, to the Langevin formulation [Eq. (4)]. This force was used to model active nuclear positioning during migration (Fig. 4), as well as nuclear re-centering after cell division (Fig. 6). Furthermore, we demonstrated initial simulations of organoid- and organs-on-chips-like platforms using constrained and unconstrained proliferation (Fig. 7). Finally, in SEM², we introduced particle-level strain (Fig. 3) and stress calculations (Figs. 5–7) downstream of the particle dynamics simulations, which enabled new statistical analyses and visualizations, such as the radial distribution function (Fig. 3) or the particle trajectories (Figs. 4 and 5).

Our goal here was not to demonstrate that SEM² can solve all problems in multiscale mechanobiology, but to introduce a framework and a software package that can be used to tackle these problems. Starting from a simple set of potentials previously derived from experimental cell stiffness and viscosity values in Eq. (3),^23,27 SEM² generated several interesting emergent mechanical behaviors that are qualitatively consistent with recent experimental findings. For example, when exposed to high stress during creep experiments, cells deformed plastically in good agreement with many experimental observations, as reviewed recently.⁴¹ Interestingly, we observed a ∼60% residual plastic deformation after high-stress creep [Figs. 3(c) and 3(d)], which is in remarkable quantitative agreement with recent data from magnetic twisting cytometry where visco-elasto-plastic deformation of 3T3 fibroblasts was analyzed.⁴² Similarly, in cell migration simulations, our modeling work showed that nuclear positioning can affect the overall motility of the cells. When the nucleus lags at the trailing edge (our no-bias condition in Fig. 4), its positioning inside the cell is mediated by the pushing from the cell rear, as reviewed before.^43,44 Yet, we also saw that the cell accelerates when the force bias is used to push the nucleus actively toward the leading edge. This is consistent with recent evidence that cells relocate their nucleus toward the leading edge to overcome obstacles or constrictions.^44,45 However, our current simulations feature a rigid nucleus that cannot be deformed, so even in the presence of reactive F_BIO (Fig. S6), we cannot model a “squeezing force” as seen in experiments.^46,47 During proliferation, we saw that the cell–cell interfaces are characterized by tensile stresses with opposite signs from the compressive stress observed inside the cell (Figs. 6 and 7). This behavior has also been studied experimentally, as much research groups have seen tensile forces at the cell–cell junction in vitro^8,48 and in vivo during embryogenesis in zebrafish.⁴⁹ Finally, the cells rearrange in response to their pre-programmed growth rate differently in confined vs unconfined environments (Figs. 7 and 8). This observation is consistent with many experimental studies in engineered cell culture platforms, including hepatocytes in a liver-on-a-chip^50,51 and cardiac muscle cells in 3D microfluidic⁵² and 2D micropatterned⁵³ cell culture platforms. Importantly, for these simulations we used only parameter sweeps (e.g., the stiffness constant in F_BIO) or parameter values that were published previously and obtained via first-principles calculations, such as top-down coarsening of bulk mechanical properties.^23,27 Instead of seeking quantitative agreement with experimental datasets, we were encouraged to see these experimentally relevant behaviors emerge spontaneously from our simulation. For example, tensile stresses measured at the cell–cell junctions are due to cortical actin interactions, which take place at a smaller scale than the particles we are currently using. Yet in SEM,² they emerge from interacting coarse-grained particles that do not explicitly represent actin or any intracellular component. This suggests that the cellular effect of the subcellular interactions is sufficiently captured by coarse-grained particles of the size we employed.

This study focuses on introducing the framework for subcellular mechanics into SEM++ and leveraged previous work for the choice of potential and particle types.^23,27 As such, it has some inherent limitations. First, it models cells as ensembles of cytoplasmic particles with one larger nuclear defect-like particle. While Newman's rheology-preserving potential provides support for using this framework in the context of cell mechanics,²³ the Langevin formulation will need to be refined to accurately recapitulate subcellular mechanics. In particular, as the particle size decreases to model smaller organelles with good resolution (endoplasmic reticulum, cytoskeletal components, nuclear envelope) we will not be able to assume that neighboring particles will be of the same type and subject to the same pairwise potentials. However, new particle types and heterogeneous potentials can be added in SEM², and their mechanical interactions can be described and analyzed with the tools presented here. In fact, the subcellular mechanics approach we propose is agnostic to the choice of potentials,²³ which opens up the possibility of using different formulations, including data-driven approaches.⁵⁴ This creates opportunities to use SEM² to integrate different experimental sources, such as live cell imaging of organelles,¹³ structural biology,⁵⁵ atomic⁵⁶ and traction force microscopy,⁵⁷ etc. To that end, it will be important to obtain models with a consistent particle size in the cell, nucleus (and other organelles), as well as the extracellular matrix⁵⁸ as these are often manipulated and analyzed in modern experiments.

A major downside of particle-based methods is their computational cost. In SEM², it took us ∼34 min on four cores (AMD EPYC 7763, 2.45 GHZ) to perform the proliferation simulations with 8 cells and ∼10,000 particles. Yet, simulating a tissue consisting of ∼28000 cells with ∼33 × 10⁶ particles took ∼2 months (supplementary video SV9). Having built SEM++ and SEM² on top of LAMMPS, we have a very scalable architecture for basic operations like updating the particles' velocity field. However, every time we add or remove a particle from the particle list (e.g., proliferation), we introduce the intrinsically serial operation of updating the particle list, which increases overhead and decreases scalability. To reduce the impact of these limitations, dedicated parallelization schemes should be developed.⁵⁹ At present, SEM² is best suited for those scenarios where SEM was also the modeling framework of choice: namely, simulating biological systems in which morphology arises as an emergent property.¹⁵ This is a particularly important point for the engineering of tissue models for pre-clinical applications (Figs. 7 and 8) as well as for future simulations of in vitro and in vivo embryonic development. We also believe SEM² can be useful in studying mechanobiology when validated frameworks for heterogeneous particle types are developed. Instead, classical continuum mechanics¹⁷ or faster discrete methods (Potts,¹⁹ vertex³) should be preferred if the biological domain does not change significantly (e.g., cardiac electrophysiology) or if subcellular mechanics is of less interest, respectively.

In conclusion, we provide SEM,² an updated implementation of SEM++ that features calculations, visualizations, and analyses of per-particle stress and strain thus enabling the study of multiscale mechanics during cell creep, migration, and proliferation experiments. The full code is available on GitHub (https://github.com/Synthetic-Physiology-Lab/sem2) for biologists, physicists, and engineers interested in modeling multiscale mechanics in shape-changing tissues.

We have performed the simulations with LAMMPS, utilizing the package SEM++. We have created a cell consisting of N_p particles in a cylindrical region of space using the “fix sem_proliferate” functionality of SEM++. We have included walls on the top and bottom of the cell along the z axis, with “fix wall/lj93” functionality of LAMMPS. For example, we applied the 9–3 Lennard-Jones potential between walls and particles to model the adhesion of particles with the fixed and movable plates in the experiments. Additionally, there is a cylindrical constraint around the z-axis using the “fix indent” functionality of LAMMPS. We have employed Brownian dynamics time integration to update the position and velocity of particles, using the “fix bd” functionality of SEM++ along with the “fix langevin” functionality of LAMMPS. Upper and lower slabs of particles have been created considering all particles in the geometric regions of thickness 1.5 μm (∼10% of the cell height) at the cell top and bottom.

To make the bottom slab fixed, the velocity of the particles of this layer was always fixed as 0. The input stress, applied to the top layer, has three different segments: (1) hold, (2) load, and (3) unload, as shown in Fig. 2. The stress, $σ= F ext A slab$⁠, applied to the top layer is 0, $σ$⁠, and 0, respectively (Fig. 2). The force to be applied on each particle during the loading segment is computed as $f ext= F ext N slab$ where $N slab$ is the number of particles of the movable slab, as depicted in Fig. 2. To quantify the slab surface area for computing the applied stress we generated a surface mesh on selected subcellular particles, employing the Gaussian density method in OVITO.⁶⁰ To obtain the total contact surface area, we cut out that part of the surface mesh that belongs to the particles forming the movable slab. We then calculated the sum of the projections of those mesh facets onto the plane perpendicular to the stretch direction (see supplementary material Fig. S2). The nominal stress can then be computed as the ratio of the total force acting on the slab and the contact surface area. The axial strain is calculated as $ε z(t)= z ( t ) − z 0 z 0$⁠, where $z 0$ is the initial height of the cell before the load segment, and $z(t)$ is the height at time t. Two more samples have been generated by uniformly rotating the sample around the x-axis.

We created a cell of N_p particles in free space using the “fix sem_PMN” functionality of SEM++ without the walls perpendicular to the z-axis or the cylindrical shell centered around the z axis. We implemented migration with the fix sem_PMN functionality of SEM++. This is a polarity-induced migration, where the direction of migration depends on the shape of the cell and is computed based on the polarity of the cell.²⁷ A migration event occurs when particles are added near the leading edge and removed from the trailing edge. Migration events can be observed in supplementary material video SV3. The probability of a migration event is inversely proportional to the parameter migration time, T_m, and proportional to N_p. In between migration events, Brownian Dynamics time integration is carried out using the “fix bd” functionality of SEM++ and the “fix langevin” functionality of LAMMPS. Similarly, a proliferation event occurs when particles are added,²⁷ as shown in supplementary material video SV5. Proliferation events are directly proportional to N_p and inversely proportional to the total proliferation time.²⁷ When the number of particles in a cell reaches 2N_p, the nucleus is replicated as explained in the main text (supplementary material video SV6).

To model active transport of the nucleus during migration, we incorporated a biasing force, which acts like a spring force between the cell center of mass (COM) and the nucleus. At each time step, the distance between the cell center of mass and nucleus, $Δ n c$⁠, is computed. We utilized the “fix addforce” functionality of LAMMPS to incorporate a force on the nuclear particle proportional to the square of the distance, $Δ n c$⁠. The constant of proportionality can be tuned to model nuclear position in various cell types, ranging from those in which the nucleus is at the cell center to those close to the cell periphery. We have incorporated this functionality of adding a $Δ n c$-dependent force bias fully from the LAMMPS input script, without necessitating any change in the source code.

We computed each particle's mean volumetric stress (in energy units) as $σ m= σ x x + σ y y + σ z z 3$⁠. As we assume homogenous particle distribution and thus particle volume, the stress derived from LAMMPS could provide for an estimate of the spatial distribution of stress inside the cell. To distinctly visualize the distribution of stresses across the subcellular material, we chose a lookup table that emphasizes contrast for most computed values. We provided an analysis of this in supplementary material Fig. S3, and the particle stress distribution of a migrating cell can be seen in Fig. 5(f). It should be noted that in Fig. 5(f), the radius of the cytoplasmic particles is smaller than d_eq/2 to visualize the stress distribution among all particles inside the cell.

We utilized the “spatial binning” modifier in OVITO to further analyze the spatial distribution of stresses. We computed the mean particle stress inside 1D bins and normalized it by the maximum magnitude of stress across all bins and plotted it in Fig. 5(f).

Instead of modeling this biasing mechanism as a mechanical constraint on the nucleus in the form of a tunable force toward the cell center, there can be alternate implementations of $F BIO$⁠. For example, an equal and opposite force can be applied on the neighboring cytoplasmic particles, as described in the supplementary material (Fig. S6, video SV10). Alternatively, we can implement a static force on the nucleus in the direction of the leading edge from the trailing edge. This can be utilized to position the nucleus near the trailing edge, at the center, or near the leading edge, as described in the supplementary material, Fig. S7. It would be also interesting to explore if an explicit squeezing force is necessary to model nuclear movement during cell migration through constricted mechanical environments.^46,47 This will be a future direction to work on when we model a deformable nucleus with multiple particles.

See the supplementary material figures: Fig. S1. Single cell creep: illustration of pre-stretched configuration. Fig. S2. Schematic illustrating the projected surface area computation. Fig. S3. Particle stresses for all particles during the proliferation simulation. Fig. S4. Axial and lateral strain versus time during the uniaxial creep experiments. Fig. S5. Plastic strain during uniaxial creep experiments with 20 Pa stress. Fig. S6. Cell migration with reaction force (RF) on cytoplasmic particles. Fig. S7. Nuclear transport and cell migration with constant FBIO. Fig. S8. Failure conditions in constant stress simulations. Supplementary material videos: Video S1. Small stress creep experiment with cells without and with the nucleus. Video S2. Large stress creep experiment with cells without and with the nucleus. Video S3. Particle polymerization and depolymerization during cell migration. Video S4. Migrating cells with no-, weak-, and strong-nuclear biasing mechanisms. Video S5. Particle polymerization during one cell division Video S6. One cell division with various visualizations of the cell and its subcellular particles. Video S7. Cell proliferation in unconstrained and constrained environments. Video S8. Cell proliferation in unconstrained and constrained environments, color-coded. Video S9. Proliferation simulation resulting in a tissue consisting of ∼ 28 000 cells with ∼ 33 × 106 particles. Particles were colored according to their cell number. Video S10. Cell migration videos without and with reaction force on cytoplasmic particles neighboring the nucleus.

We would like to acknowledge Dr. Axel Kohlmeyer and the LAMMPS community for helping and supporting the authors via the online forum, the workshop, and other initiatives.

This work was supported by the European Research Council (ERC) Starting Grant No. 852560 to F.S.P. and by the Italian Ministry of University and Research (MUR) through the PRIN project COSMIC (No. 2022A79M75) to A.R.

Alexander Stukowski and Constanze Kalcher are employed by OVITO GmbH which sells the software we used in the preparation of some of the figures.

Ethics approval is not required.

Sandipan Chattaraj: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Project administration (equal); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Michele Torre: Conceptualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Constanze Kalcher: Formal analysis (equal); Software (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Alexander Stukowski: Software (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Simone Morganti: Conceptualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Alessandro Reali: Conceptualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Francesco Silvio Pasqualini: Conceptualization (lead); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article and its supplementary material.