An agent-based model for elasto-plastic mechanical interactions between cells , basement membrane and extracellular matrix

The basement membrane (BM) and extracellular matrix (ECM) play critical roles in developmental and cancer biology, and are of great interest in biomathematics. We introduce a model of mechanical cell-BM-ECM interactions that extends current (visco)elastic models (e.g. [8, 16]), and connects to recent agent-based cell models (e.g. [2, 3, 20, 26]). We model the BM as a linked series of Hookean springs, each with time-varying length, thickness, and spring constant. Each BM spring node exchanges adhesive and repulsive forces with the cell agents using potential functions. We model elastic BM-ECM interactions with analogous ECM springs. We introduce a new model of plastic BM and ECM reorganization in response to prolonged strains, and new constitutive relations that incorporate molecular-scale effects of plasticity into the spring constants. We find that varying the balance of BM and ECM elasticity alters the node spacing along cell boundaries, yielding a nonuniform BM thickness. Uneven node spacing generates stresses that are relieved by plasticity over long times. We find that plasto-viscoelastic cell shape response is critical to relieving uneven stresses in the BM. Our modeling advances and results highlight the importance of rigorously modeling of cell-BM-ECM interactions in clinically important conditions with significant membrane deformations and time-varying membrane properties, such as aneurysms and progession from in situ to invasive carcinoma.

1.1.Biological background.The system we are going to describe is a tissue made of cells, basement membrane and extracellular matrix.
Epithelium is composed of sheets of cells that cover organ surfaces and often perform specialized functions.Each cell consists of a nucleus and cytoplasm: a gel-like substance containing organelles and the supporting cytoskeleton.The nucleus and cytoplasm are enveloped in a very thin phospholipid bilayer membrane (6 − 10 nm).This membrane maintains the integrity of the cell and performs cellular functions, such as the selective passage of substrates and attachment of the cell to the BM, ECM, and other cells.Due to their (plasto-visco)elasticity, cells can deform to create and maintain bonds to adhesive ligands (a receptor's target molecules).Normal epithelial cells are polarized: integrins on a basal side adhere to the basement membrane; the opposite apical side has relatively few adhesion molecules and is often used to release secretory products; E-cadherin molecules on the lateral sides tightly adhere to E-cadherin on neighboring cells along their lateral sides [18,33].
Epithelium is supported by the stroma: a loose connective tissue mainly composed of extracellular matrix, a scaffolding of fibers (collagen, elastin, fibronectin, etc.) embedded in a mixture of water and glycoproteins.The stroma is interlaced by blood vessels that supply vital substrates to the tissue [18].
Epithelium and stroma are separated by a semi-permeable basement membrane: a thin (< 100 nm) sheet of specialized ECM.All BMs have four major constituents: type IV collagen, laminin, perlecan and nidogen, which form a complex network of fibers.Collagen IV is a triple helix, like a molecular-scale spring.Laminin is a ligand for integrin, and nidogen forms the cross-links between laminin and collagen, to make a strong membrane with ligands for adhesion.BM acts as a mechanical support and as a barrier and filter for cells and molecules, and it regulates cell growth, differentiation and migration.Normal epithelial cells must be anchored to the BM to survive; if contact is lost, they experience programmed death (anoikis).Further, intact basement membrane also prevents tumor cells from invading the surrounding tissue [1,10,17,18].Figure 1 shows an example of tissue.1.2.Simplifications.In this paper, we focus upon the impact of cell-BM-ECM biomechanical interactions on the morphology of the BM and cell arrangement.To make this impact clearer, we generally make the following simplifications: • As in [2,3,20,26], the model can simulate interactions between cells of varying sizes and the BM.In Sections 5.1-5.2,we focus on the impact of BM and ECM elasticity and plasticity on the overall cell and BM arrangement.To facilitate this, we simulate identical cells with fixed nuclear radius R N , equivalent radius R, and maximum adhesion interaction distance R A .We demonstrate the generalized model in Section 5. 3, where we simulate the BM morphological response to two smaller cells (daughter cells following a mitosis event) and a void left after an apoptosis event.
• We assume an isotropic distribution of cell surface receptors: the entire cell surface can equally bind to ligands on neighboring cells and the BM.However, we still effectively model cell polarity by simulating a single layer of cells, similarly to polarity in the recent work by Norton et al. [23].• We do not explicitly model cell morphology, but rather track the nuclear and total cell volumes.We approximate the cells as spherical as necessary.• We simulate a fixed cell population (no proliferation, apoptosis, or necrosis) to eliminate the confounding impact of cell proliferation-induced stress and death-associated relaxation [20,22].However, it is trivial to integrate the fuller model from [20].See Section 6.
2. Model formulation.We simulate a 2-D rectangle where a basement membrane separates a single layer of epithelium from a 50 µm-thick region of ECM.Each side of the ECM is treated as a rigid wall.See Figure 2 for a schematic of the domain.We assume that the ECM is a continuous medium, while the basement membrane is discretized in a finite set of N p points: each of these interacts with the cells and the ECM, and is subject to membrane elasto-plastic forces.Any two successive points identify an edge (whose length is k ) of a rectangular BM segment, which is fully defined by the thickness T k of the membrane; the BM consists of the union of these rectangular segments (see Figure 3).Each BM point has its own coordinate system made of normal and tangential vectors, defined by the following relationships: An example of this coordinate system is shown in Figure 3. Further, each segment is assumed to have a constant area over time (because of the conservation of the material); thus, thickness is inversely proportional to the length.2.1.Forces acting upon the BM.The points of the membrane are subject to competing forces which determine their motion: they interact with the cells and the ECM through adhesion and repulsion forces, they experience elasto-plastic forces due to the membrane itself, and they undergo a drag force.See Figure 4.In order to describe the cell-BM interactions, we extend the approach introduced in [20].The interactions with ECM and elastic forces are modeled using a mass-spring model; a similar approach was used in [11] to describe cell-cell interactions.The plastic response is described in Section 2.3.we approximate the number of ligands B k : Thus, the strength of the cell-BM adhesive force is proportional to its non dimensional integrin surface receptor expression I k,i and to B k .As the cell approaches the BM, there are more adhesion receptors in contact with their ligands on the BM, and the strength of the force increases.Hence, the adhesive force between the i th cell and the k th point of the BM is modeled by: where c cba is a constant, d (x i , x k ) is the displacement between the cell i (with position x i ) and the membrane point x k , and R i A is the cell's maximum adhesion interaction distance; n cba and the potential function ϕ are described in Section 2.5.

BM-Cell repulsion.
The basement membrane resists deformation and penetration by the cells.We model this BM-Cell repulsion by: where c cbr is a constant, d(x i , x k ) is the cell-point displacement, R i N and R i are, respectively, the nuclear and the equivalent radii (as described earlier); M , n cbr and the potential function ψ are detailed in Section 2.5.
In Equations ( 3) and ( 4), the displacement d is given by Thus F ik cba and F ik cbr both act along the line joining the cell i to the membrane point k.We assume that the points of the BM move in the directions given by Equation (1); thus, the interaction forces must be projected upon these vectors.Furthermore, the tangential component is neglected, since it primarily rotates the cell but does not alter the cell's position.

BM-ECM interaction.
To describe BM-ECM interactions, we introduce a system of springs of elastic constant K ECM , which is tied to ECM density and Young's modulus.The interaction force is directed along the direction n k , and is given by the following Hooke's law: where k is the length of the spring associated to the k th point, and k,0 is its equilibrium length; K ECM,k is proportional to the amount of ECM, and is given by where K ECM,unit is the elastic constant of the spring per unit of length.
2.1.4.BM elastic force.The points of the membrane are subject to an elastic force modeled by tying each point to the neighboring ones through two springs of elastic constants K BM .This force is given by: where k+1 is the length of the spring which ties points k and k + 1, and k+1,0 is its equilibrium length; the variables k and k,0 have analogous definitions.The stiffness constant is proportional to the Young's modulus of the BM and to its density, i.e., to the area of the segments: where T k is the thickness of the k th segment.However, since the area of the segments is assumed to be constant over time, K BM,k is constant, too.
2.1.5.Drag force.Because of the viscosity ν of luminal and interstitial fluid, each moving point of the BM also experiences a drag force which is proportional to the velocity of the membrane: We express the balance of forces acting on the k th point of the membrane by Newton's second law: 2.2.Constitutive equations.As we stated in Section 1.1, the basement membrane is composed of many fibers of different materials.Thus, we can imagine each segment as a set of springs, each one representing a fiber; since, at each time step, the segment may be stretched or compressed, the arrangement of the fibers is not fixed.For example, when the segment is stretched, the number of fibers arranged along its length increases, while the thickness is proportionally reduced because the area of each segment is assumed to be constant.Further, the links between fibers are also rearranged (see Figure 5).We will not model each fiber, but we are interested in a relationship between the size of the segment and an equivalent stiffness.
Figure 5.A BM segment is comprised of elastic fibers (blue springs), which are cross-linked (green boxes).When stretched, the strain can cause random cross-links to break (denoted as broken red boxes).The associated elastic fibers relax and form new cross-links, leading to a new relaxed BM configuration (and hence a different elastic constant, equilibrium length, and thickness).
When n springs of elastic constants K 1 , K 2 , . . ., K n are arranged in series, their equivalent spring constant is if K 1 = K 2 = . . .= K n = K.Therefore, if the length of the segment increases (i.e., we add springs), K decreases; vice versa, K increases when the segment is compressed.Hence, we can state On the other hand, the springs arranged along the thickness of the segments can be considered as parallel springs, and their equivalent constant is if Thus, the equivalent stiffness is directly proportional to the thickness T of the segment: The simplest relationship which takes into account these two properties is: where c is a constant.If K 0 , 0 , T 0 are, respectively, the initial stiffness, length and thickness of the segment, where A 0 = T 0 0 is the initial area of the segment.Summarizing the last two equations, the stiffness of a BM segment is given by The ECM also has a fibrous composition, and so we can use a similar approach.The only difference is the lack of thickness; thus, the stiffness of the ECM only depends on the length of the springs: 2.3.Elasto-plastic model.In the previous sections, for simplicity we assumed that the basement membrane and the extracellular matrix experience elastic deformations.This hypothesis allows us to focus our attention on the role of the stiffness, but neglects key mechanobiology: when a tissue is stretched, the energy may be not elastically stored, but spent in breaking the molecular bonds between its components.Hence there is an internal reorganization which results-at a macroscopic level-in a material rearrangement, that is, a plastic deformation.If the imposed stress is released after some time from the beginning of the experiment, the tissue will not return to its original configuration, because in the meantime the natural configuration has changed.
To describe situations like this, in [25] the deformation gradient F is decomposed in a part that can be elastically recovered F n and a plastic part due to the reorganization of bonds between the elements of the material.The former is related to the constitutive equations for the stress T, in this case in the basement membrane.The latter requires a constitutive equation to describe how the natural configuration evolves.In [25] the following constitutive equation is suggested where [ • ] + stands for the positive part of the argument, f is an invariant measure of the stress, and τ is the value of stress above which unbinding events between single fibers occur at the microscopic scale.
In our discrete description of the membrane elements, the above equation can be adapted by saying that the rest length of the spring evolves if the force in the spring is above a threshold value, In order to take into account the fact that fiber reorganization only occurs under tension, we use where r = k/η gives an indication of the characteristic time needed to relax the internal stress.In our simulations, we assume τ = 0.
2.4.Forces acting upon the cells.The cells interact with the membrane through adhesive and repulsive forces, according to Equations ( 3) and ( 4).They are also subject to interactions with other cells and to a drag force, which are modeled as in [20].

Cell-cell adhesion.
Cell adhesion is essential for many physiological functions (survival, proliferation, differentiation, migration) and pathological conditions (inflammation, metastasis, atherosclerosis) [33].Adhesion molecules on a cell's surface bond with ligands on neighboring cells.These junctions are responsible for maintaining the structural integrity of tissues.Similarly to cell-BM adhesion, the strength of the interaction increases as the cells are closer, because more surface area (and the receptor-ligand pairs) is in direct contact.The force imparted by cell j on cell i is modeled by where c cca is a constant, f ij describes the specific molecular biology of the adhesion; n cca and the potential function ∇ϕ are described in Section 2.5.We use Cell-cell repulsion.Cells mechanically resist compression because of the structure of their cytoskeletons, the incompressibility of cytoplasm and the surface tension of their membranes.Hence, we introduce a repulsive strength which is zero when the cells are just touching, and increases rapidly as they are pressed together, particularly when nuclei are in close proximity: where c ccr is a constant, R i N and R i are, respectively, the nuclear and the equivalent radii (as described earlier); M , n ccr and the potential function ∇ψ are detailed in Section 2.5.

Drag force.
Similarly to the BM, the cells experience a drag force during their motion in the interstitial fluid: Summarizing these interactions, we can express the balance of forces acting on the cell i: . Balance of forces for cell i.

Potential functions.
In our agent-based model, we model mechanical interactions between the cells and the basement membrane using potential functions.In this Section, we introduce the functions ∇ϕ and ∇ψ which describe, respectively, mechanical adhesion and repulsion.For the adhesion potential, let R A be the cell's maximum adhesive interaction distance.For any n ∈ N, define ϕ first by its gradient: where n is said exponent of the potential.Similarly, we define the repulsion potential through its gradient.If R N is the nuclear radius, R is the cell's radius, M ≥ 1 is the cell's maximum repulsive force (i.e. the maximum value of |∇ψ|) and m is a fixed nonnegative integer, define where The potentials ϕ and ψ can be obtained by integrating with respect to |r|, but it is not necessary for our model.Note that ϕ, ψ and their derivatives have compact support, in order to model the finite interaction distance; their modulus is maximum when |r| = 0, but their signs are different because they are opposite forces.
The left panel in Figure 7 shows the balance of adhesive and repulsive forces acting on two identical cells as a function of |r|, the distance between the centers of their nuclei.We can see that there is no interaction if |r| > R A , while the balance is negative (i.e.there is repulsion) for small values of |r|: this prevents overlapping of cell nuclei.When |r| = r 0 , adhesion and repulsion forces balance each other: this is the equilibrium distance between the two cell centers when they do not interact with other cells or with the basement membrane.
Similarly, the right panel in Figure 7 shows the balance of adhesive and repulsive forces acting on a cell interacting with the BM: in this case, |r| is the distance between the center of the nucleus and the membrane.When this is the only interaction force acting on the cell, the equilibrium cell-BM distance is r 0 .

2.6.
Inertialess assumption and model summary.We introduce the inertialess assumption: we suppose that the forces equilibrate quickly, hence |m i vi | 0. This assumption is made both for cells and BM segments; thus, Equations ( 11) and ( 26) can be approximated F = 0 and solved for the velocity.In detail, defining The model we use for the points of the BM is where x k,0 describes the initial discretized positions of the BM.We assume that initially the basement membrane is in an equilibrium configuration, hence the equilibrium lengths of the springs can be evaluated through x k,0 .Similarly, in order to describe cells' motion, we define and use the following model: where x i,0 describes the initial positions of cells' centers.
3. Numerical methods.In this Section, we describe the algorithm implemented in MATLAB.31) and (33) using a finite difference scheme with a fixed time step: (d) Update springs: Given the new BM position, evaluate the new springs' length; the elastic constants of the BM and the ECM are updated according to Equations ( 18) and (19).(e) Update simulation time: Increment t by ∆t.
In our simulations, we chose a 2-D rectangular domain whose length is 194 µm.The ECM region is 50 µm thick, which is approximately half the distance between two breast ducts.Assuming N p = 100, the BM is discretized in N p + 1 segments, which are 1.9404 µm long and, according to Section 1.1, 100 nm thick.The model can readily be extended to 3-D by describing the motion of vertices of an irregular triangular mesh.For what concerns the time grid, the time step is ∆t = 0.02 minutes and the total time is t max = 60 minutes.
As a boundary condition, we suppose that there is no adhesion between the cells and the computational boundary, while repulsion force is necessary to maintain cells in the computational domain.4. Coefficients and parameters.The coefficients used in our model are principally taken from [20]; we only modified the cell-BM interactions, and introduced the elastic constants of the springs.
Before investigating the coefficients, let us introduce the initial configuration of the system.Assuming a density ρ = 0.0031 cells/ µm 2 , the equilibrium distance between two adjacent cells' centers is [20]: Since our aim is to focus on cell-BM mechanics, we arrange the cells to satisfy Equation (36), as shown in Figure 8.
Similarly to [20], we want to calibrate cell-BM interactions so that the equilibrium distance between the cell's center and the membrane is d = s 2 = 9.7023 µm by balancing the adhesion and repulsion forces at the equilibrium distance d: As R N < d < R, we can replace the potential functions and evaluate the ratio c cba /c cbr : In this way, we obtain the constraint c cba = β c cbr , but we still have a degree of freedom.Assuming c cba = 10c cca as in [20], the coefficients of cell-BM interactions are fully defined.An exhaustive list of all the parameters used in our model in given in Appendix A.
For what concerns the stiffness of the springs, in our simulations we compared the behaviour of the membrane for different values of the elastic constants of the BM and the ECM.Since the basement membrane has a density bigger than the extracellular matrix, we assume K ECM,unit < K BM,unit .
5. Simulation results.The initial configuration used in all our simulations is shown in Figure 8; the black circles represent the N p = 100 points which discretize the basement membrane.We simulate 60 minutes in each case.
Figure 8.The initial configuration used in our simulations.

Impact of BM and ECM elasticity (without plasticity).
In order to focus our attention on the elastic cell-BM interactions and on the role of K ECM and K BM , we first neglect plasticity and set r = 0; in the following sections, elastoplastic deformations will be considered.Further, simulations are run for 60 minutes, which is a time scale much smaller than the typical one for mitosis and apoptosis, allowing us to neglect macroscopic growth.For simplicity, in this Section we will denote K ECM,unit with K ECM , and K BM,unit with K BM .
The first result we introduce shows the final configurations when K ECM = 10 −3 is fixed and the stiffness of the membrane varies (see Figure 9).Because of the smallness of the ECM elastic constant, the membrane is quite deformed.The BM segments near the boundary of the cells are compressed, while the farther ones are stretched (Figure 9a); as K BM increases the segments of the BM are slightly less dense on the boundary of the cells, resulting in a lower stretch far from the cells (Figure 9b): as example, when K BM = 2K ECM the maximum length of the BM segments is 4.8191 µm, and when K BM = 20K ECM it is 4.5992 µm; the initial length was 1.9404 µm.The strain of the BM must be monitored since, if excessive, it may lacerate the membrane.In the second set of simulations we chose K ECM = 10 −2 .In Figure 10a the ECM constant is dominant, and the final deformation of the membrane is very slight.Conversely, in Figure 10b, BM stiffness dominates, and the membrane is more deformed near the boundary of the cells.Similarly to the previous case, as K BM increases the segments of the BM are distributed more evenly around the cells.Because of the larger deformations, the segments farther from the cells are more stretched than the previous case: the maximum length of BM segments ranges between 5.7839 µm and 6.3279 µm, and increases as K BM decreases.
Figure 11 shows the final configuration when K ECM = 10 −1 .Once again, for small values of BM stiffness, ECM is dominant; in this case, because of the high value of the ECM constant, the final shape of the membrane is very similar to the initial one (compare Figure 11a to 10a).As K BM increases, the segments of BM are more evenly distributed on the cells' surface (Figure 11b).In Table 1 we summarize the maximum and mean length typical ranges of the ECM springs at the final configuration.
Finally, to fully understand the role of ECM stiffness, we ran a simulation setting K ECM = 0 and K BM = 1.The configuration at time = 60 minutes is shown in Figure  12.The membrane is very deformed: it is wrapping around the boundaries of the cells and all its segments are evenly distributed.Thus, we can state that the final shape of the membrane is due to the dominating parameter between K ECM and K BM , and high values of K ECM prevent the deformation; K BM also determines the distribution of the BM segments on the surface of the cells.
Figure 13 shows the mean and the maximum length of BM and ECM springs.When K ECM = 10 −3 (red curve), BM stiffness is dominant: the segments of membrane are quite evenly distributed on the boundary of the cells for all the values of K BM , and the maximum and mean BM strain are almost constant.ECM strains are slightly increasing because, as K BM increases, the membrane segments farther from the cells are pulled by the adjacent segments with an increasing strength.When K ECM = 10 −2 (blue curve) there is a discontinuity for K BM /K ECM ≈ 7: it is the transition between the ECM and the BM dominance.For smaller values of the ratio, K ECM is dominant and the membrane is slightly deformed: this results in small values of the maximum ECM strain.Conversely, when K BM is dominant the deformations of the membrane are larger, and the ECM strain is higher.The mean ECM strain behaves in the same way.Further, in both the cases, the maximum BM strain is decreasing because the increase of K BM results in a different final disposition of the BM segments: they are less compressed around the cells, and this results in a lower stretch for the segments farther from the cells.Finally, a similar transition also occurs for K ECM = 10 −1 (black curve).The maximum ECM strain is slightly increasing as K BM increases, but the mean strain is almost flat because of the small deformations in the BM shape (see Figure 11).Also, as the membrane stiffness increases, the segments are more evenly distributed on the boundary of the cells, thus the maximum strain of the BM is smaller.
Figure 14 shows the strains for different values of the ECM constant, with a fixed ratio K BM /K ECM .The black and green curves represent simulations in which BM stiffness is dominant: as the ECM constant increases the membrane is affected by smaller deformations, thus both the BM and the ECM strains decrease.The red curve represents a simulation in which ECM stiffness is dominant: as K ECM increases, the BM segments farther from the cells are more stretched, hence the strain of these springs is increasing.
Finally, Figure 15 shows the strains for different values of K BM .For small values of K ECM the membrane stiffness is dominating, but the influence of ECM is getting stronger, thus the mean BM strain and the maximum ECM decrease; when K ECM dominates, there are no significant differences in the shape of the membrane, and the strains are almost flat.

Impact of BM and ECM plasticity.
In this Section we now simulate the full elasto-plastic model.The numerical method is conveniently modified: in the main loop, when the BM and ECM springs are updated (step 2d), their new equilibrium lengths are evaluated through Equation ( 22) and a finite difference scheme: We use the same initial configuration of the previous simulations (it is shown in Figure 8); for each simulation, we use the same value of r, both for BM and ECM.We simulate the case K ECM = 10 −3 .In Figure 16 the maximum and mean strain of the BM and the ECM are shown: each line represents a different value of r.For r ≤ 10 −3 there are no substantial differences between elastic and elasto-plastic deformations, both in the membrane and ECM strains: hence, when r ≤ 10 −3 plasticity is negligible.We note that as the plastic response increases, the maximum BM strain decreases, while the mean BM strain slightly increases: because of plasticity, BM segments are distributed more evenly on the boundary of the cells, and the farther ones are less stretched (as shown in Figure 17).Further, since the membrane is more deformed as r increases, ECM strains are increasing.Figures 16c and  16d show that for a fixed value of the ratio K BM /K ECM , the logarithm of the strain is a linear function of the plastic response, thus the strains vary exponentially.
In Figure 18 we show the results of a set of simulations in which the constitutive equations introduced in Section 2.2 are switched off.While there are no significant differences in the maximum strain, both for BM and ECM (the curves are almost overlapped), for high values of K BM the constitutive equations increase the mean strain: Equations ( 18) and (19) show that the stiffness is inversely proportional to the length of the springs, hence as the length increases, the stiffness decreases, resulting in a larger deformation.Because of this, since the strain is larger as the BM constant increases, the displacement between the red and the blue curve increases too.Finally, this displacement is larger in the mean BM strain, because it is influenced by the thickness of the BM segments: as the strain increases, the thickness decreases because the area of each segment is constant over time.From Equation 18, the stiffness decreases as the thickness decreases.Hence, when BM

Mean relative strain of the ECM springs
Figure 14.Relative strain at time = 60 minutes for different values of K ECM ; each curve represents a different value of the ratio KECM KBM .The legend is the same for all the pictures.is stretched, it becomes thinner, its stiffness decreases and it experiences larger deformations.

5.3.
Model test near proliferating and apoptosing cells.We now briefly demonstrate the model when simulating cells with nonuniform size, as would be expected in growing tissues in vivo.We replace one cell with two daughter cells with half volume to simulate two cells in G 1 phase immediately after a recent mitosis event.We also remove one of the cells to simulate the gap left in the epithelium after an apoptosis event.All other parameters are the same as before, with K ECM = 10 −3 , K BM = 2K ECM , and r = 10 −2 .In Figure 19, we see that the simulation behaves as expected: the basement membrane conforms to the smaller post-mitotic cells, and the membrane is relativley flat in the cell-free fregion.In ongoing work, we are integrating the full agent model [20] with the BM/ECM biomechanics model presented here, which will allow an investigation of the impact of tumor biomechanics on the evolving BM morphology.See further discussion below (Section 6).

Locally-weakened BM causes passive cell protrusion into the stroma.
We close by simulating the impact of a locally-weakened BM, as would be expected in cancer cell invasion [4].To do this, we reduce the BM thickness to 10 nm near the 5th cell, and leave the remainder of the BM 100 nm thick as before.This is a simplified model of localized BM degradation by diffusing matrix metalloproteinases (MMPs, such as MMP9) or cell membrane-bound MMPs (e.g., MT2-MMP), as is observed during cancer cell invasion [34].The reduced BM thickness reduces its stiffness (according to Equation 18), resulting in a smaller elastic force.Consequently, the BM-ECM elastic interaction is stronger relative to the BM elastic force along the weakened section, resulting in a greater membrane deformation towards the stroma.Because the cell-BM adhesive force is still intact, the cell is pulled along with the membrane into the stroma, with a final displacement of 5.12 µm from its starting position after 60 minutes of simulation.See the plot at t = 60 minutes in Figure 20.It is striking that protrusion of tumor cells into the stroma can be facilitated by passive adhesive, repulsive, and elastic force responses to the weakened BM alone, without need for additional motile forces.Including active motile forces would only increase this effect, potentially including a rupture of the membrane and full invasion of the stroma.
6. Discussion and future work.In this work, we developed an agent-based model to describe elasto-plastic mechanical interactions between cells and a basement membrane.Our BM/ECM model is broadly applicable to current lattice-free agent-based cell models, such as [20].We also note that our mechanics model should be compatible with other models that simulate cell morphology, such as the subcellular element model [28,29] and cellular Potts models [12,13,24,30,31,32].In this last case, since cellular Potts models are based on an energetic formulation, one should use a similar framework to handle cell-ECM interactions and basement    membrane deformation and remodeling.The former aspect is already included in standard cellular Potts models.For the latter, one should introduce an elastic energy for the basement membrane possibly in a discretized form, and take into account suitable plastic phenomena related to network remodeling, similar to those presented here in Figure 5.
To eliminate confounding influences such as strains arising from rapidly and heterogeneously-proliferating tissues (e.g., cancer), we simulated a constant number of cells over time, neglecting proliferation and death.This allowed us to more clearly examine the relative impact of elastic and plastic effects in the BM and ECM.While our results apply most directly to slowly-growing tissues, or to short      The 100 nm BM (red segments) is weakened to 10 nm thick (blue segments) near the 5th cell.The altered BM mechanical properties lead to increased membrane deformation.After 60 minutes, the cell is displaced 5.12 µm into the stroma (plotted above).Inclusion of an active motile force would amplify this result, potentially including membrane rupture and full invasion into the stroma.times for quickly-growing tissues, the BM/ECM model is capable of simulating interactions with large and small cells.See Section 5.3.Adding the volume change dynamics of proliferation, apoptosis, and necrosis, as well as cell number changes from proliferation and cell death, will be straightforward, allowing investigation of the dynamics of tumor-basement membrane-ECM biomechanical interactions.We are currently integrating the BM/ECM mechanics model presented here with the full agent-based model of Macklin et al. [20] to investigate the impact of basement membrane mechanics on in situ and invasive carcinoma progression.Indeed, recent work by Kim et al. has demonstrated a significant impact of basement membrane mechanics on DCIS progression [16].Conversely, we are also interested in the impact of stresses generated by proliferating cells and mechanical relaxations generated by cell death [20,22] on the evolving basement membrane morphology; these can also be studied by integrating the BM/ECM biomechanics model introduced in this paper with the full agent model in [20].
We ran simulations to test the key parameters-the elastic constants of the BM and ECM-across several orders of magnitude (Section 5.1).We found that the final shape of the membrane is determined by the relative balance of K BM and K ECM .When K BM dominates, we observe large deformations of the BM into the lumen between cells, in an attempt to conform to the cell morphology.When K ECM is dominant, deformations are much more constrained.K BM also determines the distribution of the BM segments on the surface of the cells.For small K BM , the BM elastic force cannot overcome cell-BM adhesion and BM-ECM elastic force to evenly redistribute BM nodes along the cell-BM contact area; larger values of K BM facilitate this redistribution of BM nodes.In the elasto-plastic model, the rate of relaxation r plays an important role, too: increasing the plastic relaxation permits greater BM deformation.Because there are no external applied forces in this model, the center of mass of the cells and BM nodes is preserved.Hence, as BM nodes are pulled towards the lumen, the cells are pulled slightly into the stroma.As a result, BM-ECM springs are elongated (positive strain) between cells, and slightly compressed (negative strain) nearest to the cells.This has implications for the distribution of BM nodes along its length, particularly when considering the impact of the cell morphology.
The cell morphology plays a critical role in the distribution of BM nodes.For purely elastic deformations (or for shorter time scales relative to the plastic relaxation rate), the combination of cell-BM adhesion and BM elasticity drives the BM to conform to the cell morphology, particularly for weaker values of BM-ECM elasticity.For a circular (or spherical) cell morphology, the cell-BM forces equilibrate at a fixed distance from the cell center, and the BM elastic force redistributes the BM nodes along this distance, with fixed spacing when K ECM = 0 (see Figure 12).For K ECM = 0, a nonuniform force BM-ECM force is introduced across the cell-BM contact area: BM nodes near the edge (where there is net positive ECM strain) are pulled towards the stroma, which tends to "sweep" them back towards the center of the contact area.Nodes near the center of the contact area (where there is net negative ECM strain) are pushed out from the stroma, which tends to sweep them out towards the edge of the contact area.The net result is the clustering of BM nodes between the center and edge of the contact area (see . Hence, we see that the cell shape strongly influences the distribution of BM nodes.However, cells can also experience plasto-viscoelastic responses to strains.The restorative force of the ECM springs should impart a lumen-directed force near the center of the cell-BM contact area where ECM springs are compressed, and a stroma-directed force near the edges (where ECM springs are elongated), causing the cell to flatten.This, in turn, should facilitate a more even redistribution of the BM nodes along the cell-BM contact area.This can be addressed in our modeling framework by increasing the cell adhesive interaction distance (approximates cell deformation [20]) and decreasing the relative strength of the cell-BM interactive forces.Other approaches might include ellipsoidal cell morphology approximations [6,16], or applying our elasto-plastic model to the cell morphology.In this case, the "BM" springs represent the cell membrane, and the "ECM" springs connect the membrane to the center of mass or nuclear envelope.An alternative approach would be to model each cell with multiple adhering agents, similarly to Newman's subcellular element model [28,29]; in this case, we would use the same potential functions for adhesion and repulsion between subcell agents, and only allow boundary elements to adhere to neighboring cells and the basement membrane.We note that such an approach would trivially model cell polarization by only allowing a subset of the boundary subcell agents to adhere.We are currently investigating these approaches.
The non-uniform membrane node distribution also highlights the potential importance of nonlinear, non-Hookean spring effects in the basement membrane.For a dense, tightly cross-linked matrix like a membrane, compression may be energetically unfavorable compared to elongation, and may promote (neglected) bending forces instead.If bending forces are not permitted (e.g.due to the elastic interactions with the ECM), the BM springs would require K BM to increase quickly with eq − .This should result in less clustering of BM nodes along the cell-BM interface, and would also be more consistent with our current plastic reorganization model (plastic response is disallowed during compression).We will consider these nonlinear spring effects in future work.Now that we have identified the range of behaviors predicted by the model, we must work to further constrain and calibrate the parameter values.Experiments should be performed to understand the elastic forces of the membrane and ECM, and their rate of plastic reorganization.Finally, in order to have more accurate predictions, we should introduce proteolytic degradation of the membrane and ECM: In response to various microenvironmental signals, some cell may secrete matrixmetalloproteinases (MMPs), a group of enzymes that promote degradation of collagens and other matrix components during normal and pathologic tissue remodeling.[5,14,15].In case of tumor growth, MMPs may be overexpressed and unresponsive to signals in tumor cells, tumor-associated fibroblasts (TAFs), and tumor-associated macrophages (TAMs), resulting in extensive degradation and modification of the basement membrane [15].As we saw in Section 5.4, this creates a significant link between ECM and BM deformations,stresses generated by rapidly proliferating tumor cells, and MMP-regulating signaling between tumor, stroma, and immune cells.We are currently investigating a multiscale model of MMP transport and proteolytic degradation of the BM and ECM capable of resolving the typical 100 nm BM thickness and ∼10 µm MMP transport length scale [7].The combination of membrane deformation, increasing strains, and proteolytic membrane and matrix degradation can create tears in the BM, allowing tumor cells to invade the stroma.This allows progression from in situ to invasive carcinoma, ultimately resulting in metastatic disease.

Figure 1 .
Figure 1.Example of tissue: integrins on the basal side adhere to ligands in the BM, and E-cadherin on the lateral sides adheres to E-cadherin on adjacent cells.

Figure 2 .
Figure 2. The domain is a single layer of 10 cells (pale blue with dark blue nuclei) in contact with the BM (red line), separated from the ECM (orange region).

Figure 3 .
Figure 3. Left: A sample discretization of the BM.x k and x k+1 identify the k-th rectangle, with length k and thickness T k .Right: Normal and tangential vectors for point x k ; see Equations 1.

Figure 4 .
Figure 4. Balance of forces for the point k of the BM.

Figure 7 .
Figure 7. Balance of adhesive and repulsive force in cell-cell interaction (left panel) and in in cell-BM interaction (right panel).

1 . 2 .
Initialization Routines (a) Global Variables: Define the time grid by choosing a time step ∆t and a final time t max , and fix the coefficients c used to express the interaction forces.Coefficients and parameters are initialized according to the values listed in Appendix A. (b) Initialize cells: Set the position of the centers of the N c cells, and create the variables which describe their geometrical features.(c) Initialize BM: Choosing a number of points N p , the basement membrane is discretized and placed at the required position.The nodes of discretization are equally spaced (the script determines automatically the distance h p ), while the first and the last nodes are placed hp 2 from the computational boundary walls: this is an attempt to reduce the influence of boundaries.(d) Initialize springs: Evaluate the initial length of the springs used to model BM-ECM interactions and the BM elastic forces.The stiffness coefficients for BM and ECM are also set.Main program loop: While t < t max (a) Distances: Evaluate distances between each pair (i, j) of cells, and each cell-BM point pairing (i, k).Euclidean distances are used to evaluate the strength of the interaction, but distances over x and y axes are also necessary, in order to evaluate the two components of the forces.(b) Interactions: Evaluate the interaction forces using the relationships introduced in Sections 2.1 and 2.4.(c) Update BM and cell positions: Solve Equations (

Figure 9 .
Figure 9. Configuration at time = 60 minutes for K ECM = 10 −3 and different values of K BM .

Figure 13 .
Figure 13.Relative strain at time = 60 minutes for different values of K BM ; each curve represents a different value of K ECM .The legend is the same for all the pictures.

Figure 15 .
Figure 15.Relative strain at time = 60 minutes for different values of K ECM ; each curve represents a different value of K BM .The legend is the same for all the pictures.

Figure 16 .
Figure 16.Relative strain at time = 60 minutes for K ECM = 10 −3 and different values of K BM ; each curve represents a different value of r.The legend is the same for pictures (a), (b), (c), (d).

Figure 17 .
Figure 17.Comparison between two final configurations for K ECM = 10 −3 and K BM = 2K ECM .Picture (a) is the result of elastic deformations; in pictures (b) and (c) elasto-plastic deformations are introduced, respectively with r = 10 −3 and r = 10 −2 .

Figure 18 .
Figure 18.Relative strain at time = 60 minutes for K ECM = 10 −3 and different values of K BM .The legend is the same for all the pictures.

Figure 19 .
Figure19.Simulation of the basement membrane following apoptosis (the gap) and a recent mitosis event (the two smaller cells).The BM conforms to the varied cell sizes, as expected.

Figure 20 .
Figure 20.Passive protrusion of cells into the stroma near weakened BM segments: The 100 nm BM (red segments) is weakened to 10 nm thick (blue segments) near the 5th cell.The altered BM mechanical properties lead to increased membrane deformation.After 60 minutes, the cell is displaced 5.12 µm into the stroma (plotted above).Inclusion of an active motile force would amplify this result, potentially including membrane rupture and full invasion into the stroma.

Table 1 .
−1and different values of K BM .Maximum and mean length of the ECM springs at the final configuration: for each entry, first value uses K BM = 2K ECM , second value uses K BM = 20K ECM .

Table 2 .
Parameters value shared by all simulations.