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

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.


INTRODUCTION
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. 37][8] New tools in molecular biology 9,10 and optical engineering [11][12][13] have shown us that mechanical forces across spatial scales modulate morphogenesis. 14Yet, we have a limited ability to computationally model these interactions. 15rom a modeling perspective, morphogenetic events are characterized by complex tissue rearrangements driven by the proliferation and migration of millions of cells. 2,9,16In 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. 17Alternatively, cells can be modeled as discrete objects interacting on a lattice, 15,18,19 or freely in space. 3,20All 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 19 ) or limited to cell membranes (e.g., vertex modeling 3 ).Yet, the mechanical forces at play in the intracellular space, or cytoplasm, can influence the mechanical behavior of entire tissues. 21,22To 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 2 .
SEM is a discrete, off-lattice framework that treats cells as mechanically competent agents 23,24 capable of biologically relevant behaviors, such as flowing, 25 migration, 26 proliferation, 27,28 and epithelia shaping. 29,30In 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,24Later, Koumoutsakos introduced SEMþþ, 27 a Cþþ implementation of SEM that leverages the open-source and strongly supported Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) library 31 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,32Together these results have made SEM the framework of choice for modeling cell shape as an emergent property of subcellular interactions. 33Yet, a treatment of mechanics across spatial scales in SEM is missing.
In SEM 2 , 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 2 against the previous results and demonstrate SEM 2 capabilities for more detailed analysis of multiscale mechanics.Then, we used SEM 2 to simulate cell proliferation in unconstrained and constrained conditions to showcase how SEM 2 could be used for designing preclinical models such as organoids and microfluidic organs-on-chips, respectively.SEM 2 and the input scripts for all presented simulations are available on our GitHub at https://github.com/Synthetic-Physiology-Lab/sem2.

SEM 2 SCALABLY RECAPITULATES CELL RHEOLOGY IN CREEP EXPERIMENTS
Simulating cells and tissues with SEM implies solving the Langevin equation for an arbitrary number of subcellular particles (N p ), 23 g_ In Eq. ( 1), g 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, n 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 In Eq. ( 2), u 0 is the potential's well depth, q and a are the scaling and shifting factors, and d eq 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þþ. 23To 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 In Eqs.(3a)-(3d), p d is the sphere close packing density.R cell , j 0 , and g 0 are the radius, stiffness, and viscosity of a cell, respectively; finally, k ¼ 0:75 is a tuning coefficient. 23,27 SEM 2 , 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,35As 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,24hen, 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. SEM 2 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 prestretched 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 prestretched 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).
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 2 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].

SUBCELLULAR MECHANICS DURING CREEP EXPERIMENTS WITH SEM 2
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, d eq $1.8 lm, as expected based on the rheology-preserving applied forces [Fig.1(b) and Eqs. ( 1) and ( 2)]. 27Instead, many particles were strained away from their equilibrium locations during the stretch [light green distributions, Figs.

MODELING MIGRATION AND ACTIVE NUCLEAR TRANSPORT WITH SEM 2
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. 26Here, 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 T m . 27When T m is large with respect to the simulation duration (T sim ), 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 T m 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 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 4.64 lm spring-like constraint.First, we introduced the variable D nc , 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 D nc 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 2 .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 D nc 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.

SUBCELLULAR MECHANICS DURING CELL MIGRATION AND NUCLEAR TRANSPORT WITH SEM 2
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 (r m ¼ rxx þryyþrzz 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 colorcoded 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. 38We 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 perparticle 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 2 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.

MODELING CELL PROLIFERATION IN SEM 2
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. 39To 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 2 uses the nuclear-centering bias introduced in the previous section (see Methods).First, it creates the second nuclear particle 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 colorcoded 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 to the first at the center of the cell, thus modeling the increased DNA density observed before cell division. 40Then, 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, 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 The cell grows (from time, t 0 ) with an increase in the number of particles, followed by duplication of the nucleus (t 1 ), resulting in two daughter cells (t 2 ).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.
particularly interesting because, so far in SEM, information about membranes has been prescribed with dedicated membrane particles or obtained algorithmically. 27,32

MODELING ORGANOIDS AND ORGAN-CHIPS WITH SEM 2
7][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 2 can be used to gain insights into multiscale and subcellular mechanisms of cell proliferation in traditional and engineered cell culture platforms.

DISCUSSION
We presented SEM 2 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 27 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 recentering 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 2 , 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 2 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 2 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. 41Interestingly, 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. 42Similarly, 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,44Yet, 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,45However, 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,47During 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. 49Finally, 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-achip 50,51 and cardiac muscle cells in 3D microfluidic 52 and 2D micropatterned 53 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 firstprinciples calculations, such as top-down coarsening of bulk mechanical properties. 23,27Instead 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,27As 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 rheologypreserving 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 SEM 2 , 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. 54This creates opportunities to use SEM 2 to integrate different experimental sources, such as live cell imaging of organelles, 13 structural biology, 55 atomic 56 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 matrix 58 as these are often manipulated and analyzed in modern experiments.
A major downside of particle-based methods is their computational cost.In SEM 2 , 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 6 particles took $2 months (supplementary video SV9).Having built SEMþþ and SEM 2 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. 59At present, SEM 2 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. 15This is a particularly important point for the engineering of tissue models for preclinical applications (Figs. 7 and 8) as well as for future simulations of in vitro and in vivo embryonic development.We also believe SEM 2 can be useful in studying mechanobiology when validated frameworks for heterogeneous particle types are developed.Instead, classical continuum mechanics 17 or faster discrete methods (Potts, 19 vertex 3 ) 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 perparticle 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.

Single-cell creep simulations
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 lm ($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, r ¼ Fext A slab , applied to the top layer is 0, r, and 0, respectively (Fig. 2).The force to be applied on each particle during the loading segment is computed as f ext ¼ Fext 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. 60To 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 e z ðtÞ ¼ zðtÞÀz0 z0 , 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, c, for each particle, is then computed as 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 lm.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.

Migration and proliferation
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. 27A 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, 27 as shown in supplementary material video SV5.Proliferation events are directly proportional to N p and inversely proportional to the total proliferation time. 27When 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, D nc , 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, D nc .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 D nc -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, In the absence of other forces, the nucleus would reach the current center of mass of the cell in one time step, dt, 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 X i¼Np i¼1 y i;j;k N p À X y i;2;k N nuc À _ y i;2;k dt 0 @ 1 A dt 2 : (8)

Per particle stress visualizations
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 where r ab is the stress; a and b take on values from x, y, and z to generate components of the tensor N nb is the number of neighbors of the I th 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 r m ¼ rxxþryyþrzz 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,47This will be a future direction to work on when we model a deformable nucleus with multiple particles.

FIG. 1 .
FIG. 1. Subcellular element modeling and mechanics (SEM 2 ).(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 surfacebased 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).
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 2 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 .
FIG. 3. Subcellular mechanics with SEM 2 .(a) and (b) Representations of cells (N p ¼ 1000) without (a) and with (b) nuclear particles, as they undergo a highstress (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. 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, D nc , is indicated in (a).(D) D nc is plotted vs time with no bias for different values of T m .Note that the nuclear drift saturates at $10 lm which is the radius of the cell.(e) D nc vs T m for the three different scenarios.(f) The full trajectories of the cell centers for the three different scenarios.

FIG. 5 .
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 colorcoded 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. 6 .
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 colorcoded with normalized volumetric stress.The cell grows (from time, t 0 ) with an increase in the number of particles, followed by duplication of the nucleus (t 1 ), resulting in two daughter cells (t 2 ).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. 7 .
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. 8 .
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.

TABLE I .
Parameters used in the model.