Modeling multiscale mechanics in shape-shifting engineered tissues, such as organoids and organs-on-chip, is both important and challenging. In fact, it is difficult to model relevant tissue-level large non-linear deformations mediated by discrete cell-level behaviors, such as migration and proliferation. One approach to solve this problem is subcellular element modeling (SEM), where ensembles of coarse-grained particles interacting via empirically defined potentials are used to model individual cells while preserving cell rheology. However, an explicit treatment of multiscale mechanics in SEM was missing. Here, we incorporated analyses and visualizations of particle level stress and strain in the open-source software SEM++ to create a new framework that we call subcellular element modeling and mechanics or SEM2. To demonstrate SEM2, we provide a detailed mechanics treatment of classical SEM simulations including single-cell creep, migration, and proliferation. We also introduce an additional force to control nuclear positioning during migration and proliferation. Finally, we show how SEM2 can be used to model proliferation in engineered cell culture platforms such as organoids and organs-on-chip. For every scenario, we present the analysis of cell emergent behaviors as offered by SEM++ and examples of stress or strain distributions that are possible with SEM2. Throughout the study, we only used first-principles literature values or parametric studies, so we left to the Discussion a qualitative comparison of our insights with recently published results. The code for SEM2 is available on GitHub at https://github.com/Synthetic-Physiology-Lab/sem2.

Morphogenesis, the process by which tissues and organs change shape to acquire novel structures and functions during embryonic development, has long fascinated life scientists1,2 and physicists.3 More recently, engineers became interested in how tissues build themselves to design better preclinical models, including organoids1,4,5 and organs-on-chips.6–8 New tools in molecular biology9,10 and optical engineering11–13 have shown us that mechanical forces across spatial scales modulate morphogenesis.14 Yet, we have a limited ability to computationally model these interactions.15 

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.17 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 models19) or limited to cell membranes (e.g., vertex modeling3). 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 SEM2.

SEM is a discrete, off-lattice framework that treats cells as mechanically competent agents23,24 capable of biologically relevant behaviors, such as flowing,25 migration,26 proliferation,27,28 and epithelia shaping.29,30 In SEM, each cell is formed by Np 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++,27 a C++ implementation of SEM that leverages the open-source and strongly supported Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) library31 to speed up 3D computations. Furthermore, algorithmic manipulations of the particles have been used to model cell migration and proliferation,27 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.33 Yet, a treatment of mechanics across spatial scales in SEM is missing.

FIG. 1.

Subcellular element modeling and mechanics (SEM2). (a) A particle-based representation of two cells where the solid lines and dashed lines denote subcellular and intercellular interactions, respectively. (b) The plot of intracellular potential as a function of interparticle distance normalized to the equilibrium distance (potential minimum). (c) Particle-based and surface-based representations of a single cell with its nucleus. (d) A representation of a dividing cell in which the cytoplasmic particles are color-coded based on the per-particle volumetric stress (see Methods).

FIG. 1.

Subcellular element modeling and mechanics (SEM2). (a) A particle-based representation of two cells where the solid lines and dashed lines denote subcellular and intercellular interactions, respectively. (b) The plot of intracellular potential as a function of interparticle distance normalized to the equilibrium distance (potential minimum). (c) Particle-based and surface-based representations of a single cell with its nucleus. (d) A representation of a dividing cell in which the cytoplasmic particles are color-coded based on the per-particle volumetric stress (see Methods).

Close modal

In SEM2, 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 SEM2 against the previous results and demonstrate SEM2 capabilities for more detailed analysis of multiscale mechanics. Then, we used SEM2 to simulate cell proliferation in unconstrained and constrained conditions to showcase how SEM2 could be used for designing preclinical models such as organoids and microfluidic organs-on-chips, respectively. SEM2 and the input scripts for all presented simulations are available on our GitHub at https://github.com/Synthetic-Physiology-Lab/sem2.

Simulating cells and tissues with SEM implies solving the Langevin equation for an arbitrary number of subcellular particles ( N p),23,
η y ̇ i , j , k = ξ i , j , k + F C y i , j , k .
(1)
In Eq. (1), η is the viscous drag coefficient and y i , j , k, y ̇ i , j , k are the position and velocity of each particle. We chose the indexes i, j, and k to denote particle number, particle type, and cell number, respectively. That is, if we assume two cells are modeled using a single nuclear particle and many cytoplasmic ones per cell [two types, Fig. 1(a)], the positions of the 100th cytoplasmic particle in cell 1 and the 300th cytoplasmic particle in cell 2 would be denoted as y 100 , 1 , 1 and y 300 , 1 , 2. With this notation, ξ i , j , k is the thermal fluctuation felt by each particle and F C y i , j , k is the net pairwise force acting on it. In particular, Newman used the following empirically defined, mild potential [ V d where d is the inter-particle distance] to enforce adhesion and volume exclusion [Fig. 1(b) and Eq. (2)],23,27
V d = u 0 e 2 ρ 1 d 2 d e q 2 α u 0 e ρ 1 d 2 d e q 2 .
(2)
In Eq. (2), u 0 is the potential's well depth, ρ and α are the scaling and shifting factors, and d e q is the equilibrium distance between particles. To limit the computational cost, these potentials focus on short-range interactions and cut off at 2.5 times the equilibrium distance. Unless otherwise noted, we used parameter values as in SEM++.23 To assign the forces in Eq. (1) based on cell-level rheology, Newman assumed N p densely packed particles connected by springs and recovered the following relationships:23,
d e q = 2 R cell ( p d N p ) 1 3 ,
(3a)
η = η 0 N p ,
(3b)
κ = κ 0 N p 1 3 ( 1 λ N p 1 3 ) ,
(3c)
u 0 = κ d e q 2 ( 8 ρ 2 ) .
(3d)
In Eqs. (3a)–(3d), p d is the sphere close packing density. R cell, κ 0, and η 0 are the radius, stiffness, and viscosity of a cell, respectively; finally, λ = 0.75 is a tuning coefficient.23,27

In SEM2, 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 Np 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).

FIG. 2.

SEM2 scalably captures cell rheology. (a) Illustration of a single-cell creep experiment. Constant stress (5 Pa) is applied to a cell sandwiched between a fixed and a movable plate, resulting in an axial strain. (b) The cell in the pre-stretched configuration is plotted using a Gaussian surface mesh (white) that encloses all particles (see also 1 °C and Methods). The top red (shaded) surface indicates the area considered for stress calculation. The black rectangles encompass the particles (shown in red) in the fixed and movable slabs. (c) Creep strain percentage vs time for a single cell with a single particle type. Solid lines indicate mean time courses across three different simulation runs, shaded areas are standard deviations. (d) The cell (white) with a single nuclear particle (blue) in the pre-stretched configuration, as in Fig. 2(b). (e) Creep strain percentage vs time for a single cell with cytoplasmic and nuclear particles, as in Fig. 2(c).

FIG. 2.

SEM2 scalably captures cell rheology. (a) Illustration of a single-cell creep experiment. Constant stress (5 Pa) is applied to a cell sandwiched between a fixed and a movable plate, resulting in an axial strain. (b) The cell in the pre-stretched configuration is plotted using a Gaussian surface mesh (white) that encloses all particles (see also 1 °C and Methods). The top red (shaded) surface indicates the area considered for stress calculation. The black rectangles encompass the particles (shown in red) in the fixed and movable slabs. (c) Creep strain percentage vs time for a single cell with a single particle type. Solid lines indicate mean time courses across three different simulation runs, shaded areas are standard deviations. (d) The cell (white) with a single nuclear particle (blue) in the pre-stretched configuration, as in Fig. 2(b). (e) Creep strain percentage vs time for a single cell with cytoplasmic and nuclear particles, as in Fig. 2(c).

Close modal
TABLE I.

Parameters used in the model.

Parameter Description Value
R cell  Cell radius  10 μ
N p  Number of particles in a cell  1000 (unless specified otherwise) 
p d  Sphere close packing density  0.74 
ρ  Scaling factor 
α  Shifting factor 
κ 0  Cell stiffness  5 × 10−3 N m−1 
η 0  Cell viscosity  5 × 10−3 N s m−1 
λ  Tuning coefficient  0.75 
m p  Mass of cytoplasmic particle  3.1 ng (for Np = 1000) 
m nuc  Mass of the nucleus  50* m p 
R nuc  Radius of nucleus (10% of cell volume)  4.64 μ
Parameter Description Value
R cell  Cell radius  10 μ
N p  Number of particles in a cell  1000 (unless specified otherwise) 
p d  Sphere close packing density  0.74 
ρ  Scaling factor 
α  Shifting factor 
κ 0  Cell stiffness  5 × 10−3 N m−1 
η 0  Cell viscosity  5 × 10−3 N s m−1 
λ  Tuning coefficient  0.75 
m p  Mass of cytoplasmic particle  3.1 ng (for Np = 1000) 
m nuc  Mass of the nucleus  50* m p 
R nuc  Radius of nucleus (10% of cell volume)  4.64 μ

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 Np and plotted the average (solid curve) and the standard deviation [shaded area, Fig. 2(c)] for each condition. Our results show that SEM2 predicts a similar axial strain time course for all simulated N p [Fig. 2(c)]. Second, following Milde et al.,27 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,36 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, deq ∼1.8 μm, as expected based on the rheology-preserving applied forces [Fig. 1(b) and Eqs. (1) and (2)].27 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 deq 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 SEM2 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.

FIG. 3.

Subcellular mechanics with SEM2. (a) and (b) Representations of cells (Np = 1000) without (a) and with (b) nuclear particles, as they undergo a high-stress (20 Pa) creep experiment. Cell membranes and particles (excluding the nucleus) have been colored based on the von-Mises shear strain invariant of the per-particle Green–Lagrangian strain tensor. The strain has been computed with respect to the unloaded configuration and the color coding was normalized to the simulation maximum to facilitate comparisons. Large and small strains are indicated by hot colors and cold colors, respectively. (c) and (d) The applied stress and the resultant axial strain are plotted for a single cell without (c) and with the nucleus (d). (e) and (f) The radial distribution function (rdf) of cytoplasmic particles in the cell without (e) and with the nucleus (f) stacked for various time points during the simulations (dark green, light green, and gray indicates rdf traces of time points in the equilibration, stretch, and relaxation phases, respectively).

FIG. 3.

Subcellular mechanics with SEM2. (a) and (b) Representations of cells (Np = 1000) without (a) and with (b) nuclear particles, as they undergo a high-stress (20 Pa) creep experiment. Cell membranes and particles (excluding the nucleus) have been colored based on the von-Mises shear strain invariant of the per-particle Green–Lagrangian strain tensor. The strain has been computed with respect to the unloaded configuration and the color coding was normalized to the simulation maximum to facilitate comparisons. Large and small strains are indicated by hot colors and cold colors, respectively. (c) and (d) The applied stress and the resultant axial strain are plotted for a single cell without (c) and with the nucleus (d). (e) and (f) The radial distribution function (rdf) of cytoplasmic particles in the cell without (e) and with the nucleus (f) stacked for various time points during the simulations (dark green, light green, and gray indicates rdf traces of time points in the equilibration, stretch, and relaxation phases, respectively).

Close modal
The creep and scaling validations are fundamental to the claim that SEM is a rheology-preserving framework for studying tissue morphogenesis. However, to actually model morphogenesis, we must also model cell migration and proliferation. One recent approach involves modeling 2D cell migration with mechano-chemically activated actuator elongation and substrate adhesion.26 Here, we leveraged the 3D simulation software SEM++,27 where the particle ensembles move as we add a particle on one side of the cell (leading edge) and remove a particle from the other (trailing edge, supplementary material video SV3). Importantly, the force balance in Eq. (1) [Fig. 1(b)] generates rapid re-adjustments of particle positions and velocities in the ensemble that restore cell rheology. The particle addition/removal frequency determines how fast a cell moves, which is captured in the model parameter migration time Tm.27 When Tm is large with respect to the simulation duration (Tsim), the number of particle addition events is small, and the cell moves a small distance from its initial position (and vice versa). However, we noticed that in SEM++ the nucleus position within the cell changed with varying Tm and could drift close to the cell trailing edge [Fig. 4(a)]. The fact that the nucleus moved very little during migration without bias in SEM++ is probably because the interaction between cytoplasmic and nuclear particles is also restricted by the cutoff, and it becomes net positive in the direction of migration only when the nucleus has drifted toward the cell edge. While this behavior is observed in certain diseased conditions,37 most healthy cells exhibit a degree of active nuclear transport that is not captured by simply modeling passive rheology. Therefore, we extended the equation of motion to
η y ̇ i , j , k = ξ i , j , k + F C y i , j , k + F BIO y i , j , k .
(4)
FIG. 4.

Nuclear transport and cell migration. (a)–(c) Cell migration snapshots at three different timepoints in the no-bias (a) weak bias, (b) and strong bias (c) conditions. The cell centroid is depicted as a black sphere, the nucleus centroid as a black cube. The distance between the nucleus and cell centers, Δ n c, is indicated in (a). (D) Δ n c is plotted vs time with no bias for different values of Tm. Note that the nuclear drift saturates at ∼10 μm which is the radius of the cell. (e) Δ n c vs Tm for the three different scenarios. (f) The full trajectories of the cell centers for the three different scenarios.

FIG. 4.

Nuclear transport and cell migration. (a)–(c) Cell migration snapshots at three different timepoints in the no-bias (a) weak bias, (b) and strong bias (c) conditions. The cell centroid is depicted as a black sphere, the nucleus centroid as a black cube. The distance between the nucleus and cell centers, Δ n c, is indicated in (a). (D) Δ n c is plotted vs time with no bias for different values of Tm. Note that the nuclear drift saturates at ∼10 μm which is the radius of the cell. (e) Δ n c vs Tm for the three different scenarios. (f) The full trajectories of the cell centers for the three different scenarios.

Close modal

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 FBIO 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 Tm that would result in nuclear drift in the original SEM++ framework (Tm =100 and Tsim = 100 s) and showed it being mitigated with different spring stiffnesses in SEM2. 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 Tm from 50 to 250 s while keeping the total duration to Tsim = 100 s. In the absence of FBIO, a degree of nuclear drift was observed for all Tm [Fig. 4(d)], but the dependence of Δ n c on Tm was reduced or completely negated by increasing the nuclear biasing strength [Fig. 4(e)]. Finally, the simulations showed that for the same choice of Tm and Tsim, 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 SEM23,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.

FIG. 5.

Subcellular mechanics during cell migration. (a) The locations of two particles, one near and another away from the nucleus, are shown before migration and after migrating for 50 s with no bias and strong bias. (b) and (c) The trajectories of the two particles during migration with no bias (b) and strong bias (c). (d) and (e) Normalized mean volumetric stress of a particle near the nucleus and another away from the nucleus during migration with no bias (d) and strong bias (e). (f) Cytoplasmic particles are color-coded per the normalized hydrostatic stress per particle in the cell before migration and after migrating for 50 s with no bias or strong bias. Mean particle stress in 1D bins is also plotted.

FIG. 5.

Subcellular mechanics during cell migration. (a) The locations of two particles, one near and another away from the nucleus, are shown before migration and after migrating for 50 s with no bias and strong bias. (b) and (c) The trajectories of the two particles during migration with no bias (b) and strong bias (c). (d) and (e) Normalized mean volumetric stress of a particle near the nucleus and another away from the nucleus during migration with no bias (d) and strong bias (e). (f) Cytoplasmic particles are color-coded per the normalized hydrostatic stress per particle in the cell before migration and after migrating for 50 s with no bias or strong bias. Mean particle stress in 1D bins is also plotted.

Close modal

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.38 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 SEM2 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 2Np, 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.39 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, SEM2 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.40 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].

FIG. 6.

Cell proliferation. (a) Depiction of the cell proliferation simulation with nuclear particle and the membrane visualized with a Gaussian surface. (b) Depiction of the same proliferation simulations with cytoplasmic particles color-coded with normalized volumetric stress. The cell grows (from time, t0) with an increase in the number of particles, followed by duplication of the nucleus (t1), resulting in two daughter cells (t2). The interfaces between cells and the peripheral particles are under tensile (positive) stress, whereas the internal particles are under compressive (negative) stress. (c) Mean particle stress in 2D spatial bins on a projected x–y plane is color-coded. The cell membrane and the interface can be distinguished based on particle stress. Colormap is a per-particle in panel B and represents the average particle-level stress in the 2D bins in panel (c). (d) Plot of mean per-particle hydrostatic stress in 1D spatial bins along the x-axis. Black arrows indicate the tensile peak corresponding to a cell–cell interface.

FIG. 6.

Cell proliferation. (a) Depiction of the cell proliferation simulation with nuclear particle and the membrane visualized with a Gaussian surface. (b) Depiction of the same proliferation simulations with cytoplasmic particles color-coded with normalized volumetric stress. The cell grows (from time, t0) with an increase in the number of particles, followed by duplication of the nucleus (t1), resulting in two daughter cells (t2). The interfaces between cells and the peripheral particles are under tensile (positive) stress, whereas the internal particles are under compressive (negative) stress. (c) Mean particle stress in 2D spatial bins on a projected x–y plane is color-coded. The cell membrane and the interface can be distinguished based on particle stress. Colormap is a per-particle in panel B and represents the average particle-level stress in the 2D bins in panel (c). (d) Plot of mean per-particle hydrostatic stress in 1D spatial bins along the x-axis. Black arrows indicate the tensile peak corresponding to a cell–cell interface.

Close modal

To study subcellular mechanics during cell division with SEM,2 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 SEM2 can be used for tissue engineering, we investigated subcellular mechanics during cell proliferation in simulation scenarios mimicking organoids1,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)].

FIG. 7.

Traditional and engineered cell culture platforms. (a) and (b) Simulations of unconstrained (a) and constrained (b) proliferation simulations during which a single mother cell undergoes three rounds of cell division. (c) and (d) Color-coded particle stress distribution and binned mean particle stresses along the vertical axis are shown for unconstrained (c) and constrained (d) proliferation. Black arrows indicate the tensile peaks corresponding to the aligned cell–cell interfaces.

FIG. 7.

Traditional and engineered cell culture platforms. (a) and (b) Simulations of unconstrained (a) and constrained (b) proliferation simulations during which a single mother cell undergoes three rounds of cell division. (c) and (d) Color-coded particle stress distribution and binned mean particle stresses along the vertical axis are shown for unconstrained (c) and constrained (d) proliferation. Black arrows indicate the tensile peaks corresponding to the aligned cell–cell interfaces.

Close modal

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 SEM2 can be used to gain insights into multiscale and subcellular mechanisms of cell proliferation in traditional and engineered cell culture platforms.

FIG. 8.

Stresses and strains in cell proliferation. (a) Unconstrained and constrained volumetric strain vs simulation time of the proliferating cells shown in Fig. 7. (b) Unconstrained and constrained axial strain along the unconstrained axis vs simulation time of the proliferating cells. (c) Unconstrained and constrained volumetric strain of a single cell vs time during the same proliferation simulation. (d) Mean particle stresses in 1D bins along the unconstrained axis for unconstrained and constrained proliferation. Black arrows indicate the tensile peaks corresponding to the aligned cell–cell interfaces.

FIG. 8.

Stresses and strains in cell proliferation. (a) Unconstrained and constrained volumetric strain vs simulation time of the proliferating cells shown in Fig. 7. (b) Unconstrained and constrained axial strain along the unconstrained axis vs simulation time of the proliferating cells. (c) Unconstrained and constrained volumetric strain of a single cell vs time during the same proliferation simulation. (d) Mean particle stresses in 1D bins along the unconstrained axis for unconstrained and constrained proliferation. Black arrows indicate the tensile peaks corresponding to the aligned cell–cell interfaces.

Close modal

We presented SEM2 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 proliferation27 by adding a general force term, FBIO, 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 SEM2, 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 SEM2 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 SEM2 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.41 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.42 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 FBIO (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 vitro8,48 and in vivo during embryogenesis in zebrafish.49 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-chip50,51 and cardiac muscle cells in 3D microfluidic52 and 2D micropatterned53 cell culture platforms. Importantly, for these simulations we used only parameter sweeps (e.g., the stiffness constant in FBIO) 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,2 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,23 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 SEM2, 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,23 which opens up the possibility of using different formulations, including data-driven approaches.54 This creates opportunities to use SEM2 to integrate different experimental sources, such as live cell imaging of organelles,13 structural biology,55 atomic56 and traction force microscopy,57 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 matrix58 as these are often manipulated and analyzed in modern experiments.

A major downside of particle-based methods is their computational cost. In SEM2, 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 × 106 particles took ∼2 months (supplementary video SV9). Having built SEM++ and SEM2 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.59 At present, SEM2 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.15 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 SEM2 can be useful in studying mechanobiology when validated frameworks for heterogeneous particle types are developed. Instead, classical continuum mechanics17 or faster discrete methods (Potts,19 vertex3) 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,2 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 Np 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.60 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 calculated the per-particle strain tensor in OVITO using the “atomic strain” modifier. The atomic strain computation in OVITO is based on the finite strain theory. The Green–Lagrangian strain tensor, E, is calculated from the deformation gradient tensor of each particle. The von-Mises shear strain invariant, γ, for each particle, is then computed as
γ = E x y 2 + E x z 2 + E y z 2 + 1 6 ( E x x E y y ) 2 + ( E x x E z z ) 2 + ( E y y E z z ) 2 1 2 .
(5)
We computed the strain tensor with reference to the initial configuration before loading, i.e., at t = 10 s, with a cutoff radius of 10 μm. For visualization in Figs. 3(a) and 3(b), we have normalized all per-particle strains with the highest magnitude of shear strain observed across all simulations. Demonstration of failure at higher stresses can be found in supplementary material Fig. S8.

We created a cell of Np 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.27 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, Tm, and proportional to Np. 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,27 as shown in supplementary material video SV5. Proliferation events are directly proportional to Np and inversely proportional to the total proliferation time.27 When the number of particles in a cell reaches 2Np, 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.

For biasing the nucleus toward the center, let us consider the equation of motion after Newton's law,
Δ n c = u n d t + 1 2 a n d t 2 .
(6)
In the absence of other forces, the nucleus would reach the current center of mass of the cell in one time step, d t, if u n and a n are the velocity and acceleration of the nuclear particle, respectively. If m nuc is the mass of the nucleus, S is a dimensionless elastic scaling constant proportional to the stiffness of a virtual spring connecting the nucleus to the cell center, the biasing force on the nucleus is defined by
F nuc _ bias = S m nuc a n = S m nuc 2 ( Δ n c u n d t ) d t 2 ,
(7)
F BIO = F BIAS ( t ) = j = 2 S m nuc 2 i = 1 i = N p y i , j , k N p y i , 2 , k N nuc y ̇ i , 2 , k d t d t 2 .
(8)
The per-particle stress has been determined utilizing the “compute stress/atom” functionality of LAMMPS. This yields per-particle virial stress, which has units of energy. Since this is an overdamped system, the contribution to energy and stress is primarily from potential energy. The stress tensor for atom I is defined as
σ a b = 1 2 n = 1 N n b r 1 a F 1 b + r 2 a F 2 b ,
(9)
where σ a b is the stress; a and b take on values from x, y, and z to generate components of the tensor N n b is the number of neighbors of the Ith atom. r 1, r 2, F 1, and F 2, are positions and forces due to pairwise interactions.

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 deq/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.

1.
N.
Gritti
,
D.
Oriola
, and
V.
Trivedi
, “
Rethinking embryology in vitro: A synergy between engineering, data science and theory
,”
Dev. Biol.
474
,
48
61
(
2021
).
2.
C.
Collinet
and
T.
Lecuit
, “
Programmed and self-organized flow of information during morphogenesis
,”
Nat. Rev. Mol. Cell Biol.
22
(
4
),
245
265
(
2021
).
3.
S.
Kim
,
M.
Pochitaloff
,
G. A.
Stooke-Vaughan
, and
O.
Campàs
, “
Embryonic tissues as active foams
,”
Nat. Phys.
17
(
7
),
859
866
(
2021
).
4.
E.
Bayir
,
A.
Sendemir
, and
Y. F.
Missirlis
, “
Mechanobiology of cells and cell systems, such as organoids
,”
Biophys. Rev.
11
(
5
),
721
728
(
2019
).
5.
H. A.
McCauley
and
J. M.
Wells
, “
Pluripotent stem cell-derived organoids: Using principles of developmental biology to grow human tissues in a dish
,”
Development
144
(
6
),
958
962
(
2017
).
6.
S. N.
Bhatia
and
D. E.
Ingber
, “
Microfluidic organs-on-chips
,”
Nat. Biotechnol.
32
(
8
),
760
772
(
2014
).
7.
G.
Wang
,
M. L.
McCain
,
L.
Yang
,
A.
He
,
F. S.
Pasqualini
,
A.
Agarwal
,
H.
Yuan
,
D.
Jiang
,
D.
Zhang
,
L.
Zangi
,
J.
Geva
,
A. E.
Roberts
,
Q.
Ma
,
J.
Ding
,
J.
Chen
,
D.-Z.
Wang
,
K.
Li
,
J.
Wang
,
R. J. A.
Wanders
,
W.
Kulik
,
F. M.
Vaz
,
M. A.
Laflamme
,
C. E.
Murry
,
K. R.
Chien
,
R. I.
Kelley
,
G. M.
Church
,
K. K.
Parker
, and
W. T.
Pu
, “
Modeling the mitochondrial cardiomyopathy of Barth syndrome with induced pluripotent stem cell and heart-on-chip technologies
,”
Nat. Med.
20
(
6
),
616
623
(
2014
).
8.
Y.
Aratyn-Schaus
,
F. S.
Pasqualini
,
H.
Yuan
,
M. L.
McCain
,
G. J. C.
Ye
,
S. P.
Sheehy
,
P. H.
Campbell
, and
K. K.
Parker
, “
Coupling primary and stem cell-derived cardiomyocytes in an in vitro model of cardiac cell therapy
,”
J. Cell Biol.
212
(
4
),
389
397
(
2016
).
9.
X.
Liu
,
J. P.
Tan
,
J.
Schröder
,
A.
Aberkane
,
J. F.
Ouyang
,
M.
Mohenska
,
S. M.
Lim
,
Y. B. Y.
Sun
,
J.
Chen
,
G.
Sun
,
Y.
Zhou
,
D.
Poppe
,
R.
Lister
,
A. T.
Clark
,
O. J. L.
Rackham
,
J.
Zenker
, and
J. M.
Polo
, “
Modelling human blastocysts by reprogramming fibroblasts into iBlastoids
,”
Nature
591
(
7851
),
627
632
(
2021
).
10.
M.
Matsuda
,
Y.
Yamanaka
,
M.
Uemura
,
M.
Osawa
,
M. K.
Saito
,
A.
Nagahashi
,
M.
Nishio
,
L.
Guo
,
S.
Ikegawa
,
S.
Sakurai
,
S.
Kihara
,
T. L.
Maurissen
,
M.
Nakamura
,
T.
Matsumoto
,
H.
Yoshitomi
,
M.
Ikeya
,
N.
Kawakami
,
T.
Yamamoto
,
K.
Woltjen
,
M.
Ebisuya
,
J.
Toguchida
, and
C.
Alev
, “
Recapitulating the human segmentation clock with pluripotent stem cells
,”
Nature
580
(
7801
),
124
129
(
2020
).
11.
B.
Yang
,
M.
Lange
,
A.
Millett-Sikking
,
X.
Zhao
,
J.
Bragantini
,
S.
VijayKumar
,
M.
Kamb
,
R.
Gómez-Sjöberg
,
A. C.
Solak
,
W.
Wang
,
H.
Kobayashi
,
M. N.
McCarroll
,
L. W.
Whitehead
,
R. P.
Fiolka
,
T. B.
Kornberg
,
A. G.
York
, and
L. A.
Royer
, “
DaXi—High-resolution, large imaging volume and multi-view single-objective light-sheet microscopy
,”
Nat. Methods
19
(
4
),
461
469
(
2022
).
12.
P. J.
Keller
,
A. D.
Schmidt
,
J.
Wittbrodt
, and
E. H. K.
Stelzer
, “
Reconstruction of zebrafish early embryonic development by scanned light sheet microscopy
,”
Science
322
(
5904
),
1065
1069
(
2008
).
13.
E.
Sapoznik
,
B.-J.
Chang
,
J.
Huh
,
R. J.
Ju
,
E. V.
Azarova
,
T.
Pohlkamp
,
E. S.
Welf
,
D.
Broadbent
,
A. F.
Carisey
,
S. J.
Stehbens
,
K.-M.
Lee
,
A.
Marín
,
A. B.
Hanker
,
J. C.
Schmidt
,
C. L.
Arteaga
,
B.
Yang
,
Y.
Kobayashi
,
P. R.
Tata
,
R.
Kruithoff
,
K.
Doubrovinski
,
D. P.
Shepherd
,
A.
Millett-Sikking
,
A. G.
York
,
K. M.
Dean
, and
R. P.
Fiolka
, “
A versatile oblique plane microscope for large-scale and high-resolution imaging of subcellular dynamics
,”
eLife
9
,
e57681
(
2020
).
14.
J.
Howard
,
S. W.
Grill
, and
J. S.
Bois
, “
Turing's next steps: The mechanochemical basis of morphogenesis
,”
Nat. Rev. Mol. Cell Biol.
12
(
6
),
392
398
(
2011
).
15.
A. H. K.
Roeder
,
P. T.
Tarr
,
C.
Tobin
,
X.
Zhang
,
V.
Chickarmane
,
A.
Cunha
, and
E. M.
Meyerowitz
, “
Computational morphodynamics of plants: Integrating development over space and time
,”
Nat. Rev. Mol. Cell Biol.
12
(
4
),
265
273
(
2011
).
16.
S.
Okuda
and
K.
Sato
, “
Polarized interfacial tension induces collective migration of cells, as a cluster, in a 3D tissue
,”
Biophys. J.
121
(
10
),
1856
1867
(
2022
).
17.
M.
Majda
,
N.
Trozzi
,
G.
Mosca
, and
R. S.
Smith
, “
How cell geometry and cellular patterning influence tissue stiffness
,”
Int. J. Mol. Sci.
23
(
10
),
5651
(
2022
).
18.
V. A.
Grieneisen
,
J.
Xu
,
A. F. M.
Marée
,
P.
Hogeweg
, and
B.
Scheres
, “
Auxin transport is sufficient to generate a maximum and gradient guiding root growth
,”
Nature
449
(
7165
),
1008
1013
(
2007
).
19.
I. M. N.
Wortel
,
I.
Niculescu
,
P. M.
Kolijn
,
N. S.
Gov
,
R. J.
De Boer
, and
J.
Textor
, “
Local actin dynamics couple speed and persistence in a cellular Potts model of cell migration
,”
Biophys. J.
120
(
13
),
2609
2622
(
2021
).
20.
P. J.
Brown
,
J. E. F.
Green
,
B. J.
Binder
, and
J. M.
Osborne
, “
A rigid body framework for multicellular modeling
,”
Nat. Comput. Sci.
1
(
11
),
754
766
(
2021
).
21.
V.
Venturini
,
F.
Pezzano
,
F.
Català Castro
,
H.-M.
Häkkinen
,
S.
Jiménez-Delgado
,
M.
Colomer-Rosell
,
M.
Marro
,
Q.
Tolosa-Ramon
,
S.
Paz-López
,
M. A.
Valverde
,
J.
Weghuber
,
P.
Loza-Alvarez
,
M.
Krieg
,
S.
Wieser
, and
V.
Ruprecht
, “
The nucleus measures shape changes for cellular proprioception to control dynamic cell behavior
,”
Science
370
(
6514
),
eaba2644
(
2020
).
22.
A. J.
Lomakin
,
C. J.
Cattin
,
D.
Cuvelier
,
Z.
Alraies
,
M.
Molina
,
G. P. F.
Nader
,
N.
Srivastava
,
P. J.
Sáez
,
J. M.
Garcia-Arcos
,
I. Y.
Zhitnyak
,
A.
Bhargava
,
M. K.
Driscoll
,
E. S.
Welf
,
R.
Fiolka
,
R. J.
Petrie
,
N. S.
De Silva
,
J. M.
González-Granado
,
N.
Manel
,
A. M.
Lennon-Duménil
,
D. J.
Müller
, and
M.
Piel
, “
The nucleus acts as a ruler tailoring cell responses to spatial constraints
,”
Science
370
(
6514
),
eaba2894
(
2020
).
23.
S. A.
Sandersius
and
T. J.
Newman
, “
Modeling cell rheology with the Subcellular Element Model
,”
Phys. Biol.
5
(
1
),
015002
(
2008
).
24.
S. A.
Sandersius
,
C. J.
Weijer
, and
T. J.
Newman
, “
Emergent cell and tissue dynamics from subcellular modeling of active biomechanical processes
,”
Phys. Biol.
8
(
4
),
045007
(
2011
).
25.
C. R.
Sweet
,
S.
Chatterjee
,
Z.
Xu
,
K.
Bisordi
,
E. D.
Rosen
, and
M.
Alber
, “
Modelling platelet–blood flow interaction using the subcellular element Langevin method
,”
J. R. Soc. Interface.
8
(
65
),
1760
1771
(
2011
).
26.
M.
Tarama
,
K.
Mori
, and
R.
Yamamoto
, “
Mechanochemical subcellular-element model of crawling cells
,”
Front. Cell Dev. Biol.
10
,
1046053
(
2022
).
27.
F.
Milde
,
G.
Tauriello
,
H.
Haberkern
, and
P.
Koumoutsakos
, “
SEM++: A particle model of cellular growth, signaling and migration
,”
Comp. Part. Mech.
1
(
2
),
211
227
(
2014
).
28.
M.
Banwarth-Kuhn
,
A.
Nematbakhsh
,
K. W.
Rodriguez
,
S.
Snipes
,
C. G.
Rasmussen
,
G. V.
Reddy
, and
M.
Alber
, “
Cell-based model of the generation and maintenance of the shape and structure of the multilayered shoot apical meristem of Arabidopsis thaliana
,”
Bull. Math. Biol.
81
(
8
),
3245
3281
(
2019
).
29.
A.
Nematbakhsh
,
W.
Sun
,
P. A.
Brodskiy
,
A.
Amiri
,
C.
Narciso
,
Z.
Xu
,
J. J.
Zartman
, and
M.
Alber
, “
Multi-scale computational study of the mechanical regulation of cell mitotic rounding in epithelia
,”
PLoS Comput. Biol.
13
(
5
),
e1005533
(
2017
).
30.
A.
Nematbakhsh
,
M.
Levis
,
N.
Kumar
,
W.
Chen
,
J. J.
Zartman
, and
M.
Alber
, “
Epithelial organ shape is generated by patterned actomyosin contractility and maintained by the extracellular matrix
,”
PLoS Comput. Biol.
16
(
8
),
e1008105
(
2020
).
31.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in't Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
, and
S. J.
Plimpton
, “
LAMMPS—A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
,”
Comput. Phys. Commun.
271
,
108171
(
2022
).
32.
K.
Lykov
,
Y.
Nematbakhsh
,
M.
Shang
,
C. T.
Lim
, and
I. V.
Pivkin
, “
Probing eukaryotic cell mechanics via mesoscopic simulations
,”
PLoS Comput. Biol.
13
(
9
),
e1005726
(
2017
).
33.
J.
Metzcar
,
Y.
Wang
,
R.
Heiland
, and
P.
Macklin
, “
A review of cell-based computational modeling in cancer biology
,”
JCO Clin. Cancer Inf.
3,
1
13
(
2019
).
34.
A.
Micoulet
,
J. P.
Spatz
, and
A.
Ott
, “
Mechanical response analysis and power generation by single-cell stretching
,”
ChemPhysChem
6
(
4
),
663
670
(
2005
).
35.
N.
Desprat
,
A.
Richert
,
J.
Simeon
, and
A.
Asnacios
, “
Creep function of a single living cell
,”
Biophys. J.
88
(
3
),
2224
2233
(
2005
).
36.
P. W.
Alford
,
B. E.
Dabiri
,
J. A.
Goss
,
M. A.
Hemphill
,
M. D.
Brigham
, and
K. K.
Parker
, “
Blast-induced phenotypic switching in cerebral vasospasm
,”
Proc. Natl. Acad. Sci. U. S. A.
108
(
31
),
12705
12710
(
2011
).
37.
G. G.
Gundersen
and
H. J.
Worman
, “
Nuclear positioning
,”
Cell
152
(
6
),
1376
1389
(
2013
).
38.
M. A. N.
Dewapriya
and
R. K. N. D.
Rajapakse
, “
Atomistic and continuum modelling of stress field at an inhomogeneity in graphene
,”
Mater. Des.
160
,
718
730
(
2018
).
39.
N.
Pavin
and
I. M.
Tolić
, “
Mechanobiology of the mitotic spindle
,”
Dev. Cell
56
(
2
),
192
201
(
2021
).
40.
P.
Risteski
,
M.
Jagrić
,
N.
Pavin
, and
I. M.
Tolić
, “
Biomechanics of chromosome alignment at the spindle midplane
,”
Curr. Biol.
31
(
10
),
R574
R585
(
2021
).
41.
K.
Molnar
and
M.
Labouesse
, “
The plastic cell: Mechanical deformation of cells and tissues
,”
Open Biol.
11
(
2
),
210006
(
2021
).
42.
N.
Bonakdar
,
R.
Gerum
,
M.
Kuhn
,
M.
Spörrer
,
A.
Lippert
,
W.
Schneider
,
K. E.
Aifantis
, and
B.
Fabry
, “
Mechanical plasticity of cells
,”
Nat. Mater.
15
(
10
),
1090
1094
(
2016
).
43.
F. J.
Calero-Cuenca
,
C. S.
Janota
, and
E. R.
Gomes
, “
Dealing with the nucleus during cell migration
,”
Curr. Opin. Cell Biol.
50
,
35
41
(
2018
).
44.
P.
Friedl
,
K.
Wolf
, and
J.
Lammerding
, “
Nuclear mechanics during cell migration
,”
Curr. Opin. Cell Biol.
23
(
1
),
55
64
(
2011
).
45.
A. D.
Doyle
,
D. J.
Sykora
,
G. G.
Pacheco
,
M. L.
Kutys
, and
K. M.
Yamada
, “
3D mesenchymal cell migration is driven by anterior cellular contraction that generates an extracellular matrix prestrain
,”
Dev. Cell
56
(
6
),
826
841.e4
(
2021
).
46.
A. L.
McGregor
,
C.-R.
Hsia
, and
J.
Lammerding
, “
Squish and squeeze—The nucleus as a physical barrier during migration in confined environments
,”
Curr. Opin. Cell Biol.
40
,
32
40
(
2016
).
47.
J.
Keys
,
B. C. H.
Cheung
,
M.
Wu
, and
J.
Lammerding
,
Rear Cortex Contraction Aids in Nuclear Transit during Confined Migration by Increasing Pressure in the Cell Posterior
(
Cell Biology
,
2022
).
48.
G.
Charras
and
A. S.
Yap
, “
Tensile forces and mechanotransduction at cell–cell junctions
,”
Curr. Biol.
28
(
8
),
R445
R457
(
2018
).
49.
K.
Bambardekar
,
R.
Clément
,
O.
Blanc
,
C.
Chardès
, and
P.-F.
Lenne
, “
Direct laser manipulation reveals the mechanics of cell contacts in vivo
,”
Proc. Natl. Acad. Sci. U. S. A.
112
(
5
),
1416
1421
(
2015
).
50.
Z.
Yang
,
X.
Liu
,
E. M.
Cribbin
,
A. M.
Kim
,
J. J.
Li
, and
K.-T.
Yong
, “
Liver-on-a-chip: Considerations, advances, and beyond
,”
Biomicrofluidics
16
(
6
),
061502
(
2022
).
51.
Y.
Nakao
,
H.
Kimura
,
Y.
Sakai
, and
T.
Fujii
, “
Bile canaliculi formation by aligning rat primary hepatocytes in a microfluidic device
,”
Biomicrofluidics
5
(
2
),
022212
(
2011
).
52.
A.
Mathur
,
P.
Loskill
,
K.
Shao
,
N.
Huebsch
,
S.
Hong
,
S. G.
Marcus
,
N.
Marks
,
M.
Mandegar
,
B. R.
Conklin
,
L. P.
Lee
, and
K. E.
Healy
, “
Human iPSC-based cardiac microphysiological system for drug screening applications
,”
Sci. Rep.
5
(
1
),
8883
(
2015
).
53.
S. P.
Sheehy
,
F.
Pasqualini
,
A.
Grosberg
,
S. J.
Park
,
Y.
Aratyn-Schaus
, and
K. K.
Parker
, “
Quality metrics for stem cell-derived cardiac myocytes
,”
Stem Cell Rep.
2
(
3
),
282
294
(
2014
).
54.
B.
Basu
,
N. H.
Gowtham
,
Y.
Xiao
,
S. R.
Kalidindi
, and
K. W.
Leong
, “
Biomaterialomics: Data science-driven pathways to develop fourth-generation biomaterials
,”
Acta Biomater.
143
,
1
25
(
2022
).
55.
K.
Beicker
,
E. T.
O'Brien
,
M. R.
Falvo
, and
R.
Superfine
, “
Vertical light sheet enhanced side-view imaging for AFM cell mechanics studies
,”
Sci. Rep.
8
(
1
),
1504
(
2018
).
56.
C. M.
Hobson
,
M.
Kern
,
E. T.
O'Brien
,
A. D.
Stephens
,
M. R.
Falvo
, and
R.
Superfine
, “
Correlating nuclear morphology and external force with combined atomic force microscopy and light sheet imaging separates roles of chromatin and lamin A/C in nuclear mechanics
,”
Mol. Biol. Cell.
31
(
16
),
1788
1801
(
2020
).
57.
V.
Maruthamuthu
,
B.
Sabass
,
U. S.
Schwarz
, and
M. L.
Gardel
, “
Cell-ECM traction force modulates endogenous tension at cell–cell contacts
,”
Proc. Natl. Acad. Sci. U. S. A.
108
(
12
),
4708
4713
(
2011
).
58.
M.
Théry
,
V.
Racine
,
A.
Pépin
,
M.
Piel
,
Y.
Chen
,
J.-B.
Sibarita
, and
M.
Bornens
, “
The extracellular matrix guides the orientation of the cell division axis
,”
Nat. Cell Biol.
7
(
10
),
947
953
(
2005
).
59.
J.
Pahlke
and
I. F.
Sbalzarini
, “
A unifying mathematical definition of particle methods
,”
IEEE Open J. Comput. Soc.
4
,
97
108
(
2023
).
60.
A.
Stukowski
, “
Visualization and analysis of atomistic simulation data with OVITO—The open visualization tool
,”
Modell. Simul. Mater. Sci. Eng.
18
(
1
),
015012
(
2010
).

Supplementary Material