AN AGENT-BASED MODEL FOR ELASTO-PLASTIC MECHANICAL INTERACTIONS BETWEEN CELLS, BASEMENT MEMBRANE AND EXTRACELLULAR MATRIX

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. Simpliftcations.

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

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.

When n springs of elastic constants K₁, K₂,…, K_n are arranged in series, their equivalent spring constant is

if $K_1=K_2=\dots =K_n=\overline{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 $K_1=K_2=\dots =K_n=\overline{K}$. 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₀, ℓ₀, T₀ are, respectively, the initial stiffness, length and thickness of the segment,

where A₀ = T₀ℓ₀ 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 $\tau$ 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 $\tau =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].

2.5. 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\in \mathbb{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 $\left|\nabla \psi \right|$) and m is a fixed nonnegative integer, define

where

The potentials φ and ψ can be obtained by integrating with respect to $\left|\text{r}\right|$, 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 $\left|\text{r}\right|=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 $\left|\text{r}\right|=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, $\left|\text{r}\right|$ 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₀.

2.6. Inertialess assumption and model summary.

We introduce the inertialess assumption: we suppose that the forces equilibrate quickly, hence $\left|m_i{\dot{v}}_i\right|\simeq 0$. This assumption is made both for cells and BM segments; thus, Equations (11) and (26) can be approximated $\sum \text{F}=\text{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.

5.1. 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, elasto-plastic 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⁻³ 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⁻². 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⁻¹. 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⁻³ (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⁻² (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⁻¹ (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.

5.2. 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⁻³. 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⁻³ there are no substantial differences between elastic and elasto-plastic deformations, both in the membrane and ECM strains: hence, when r ≤ 10⁻³ 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 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₁ 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⁻³, K_BM = 2K_ECM, and r = 10⁻². 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).

5.4. 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.