Role of contact inhibition of locomotion and junctional mechanics in epithelial collective responses to injury

In this model, the apical layer of a cell $i$ (of $N$ cells that form the monolayer) consists of a region of sites in the lattice with spin $i$ where $i\left(x,y \right)$ gives the lattice position. In the initialization of the CPM, we choose a resolution $R$ and subdivide the domain into $N$ squares each of side $R$ pixels. Here $R$ is the square root of the preferred area $R=\sqrt{{{a}_{{\rm p}}}}.$ The tissue is modeled with periodic boundary conditions and contiguity is enforced, thereby preventing cell partitioning. The energy equation used in our simulations is the same that we used previously (Noppe et al 2015) with an added apical-basal crosstalk term (Coburn et al 2016):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{E}_{i}}=KL_{i}^{2}-J{{L}_{i}}+\lambda {{({{a}_{i}}-{{a}_{{\rm p}}})}^{2}}+s{{\left| \delta \boldsymbol{r}(t)~ \right|}^{2}}\nonumber \end{align} \tag{ 1 }$$

where ${{E}_{i}}$ is the energy of the cell $i$, $K$ the strength of the junction contractility, ${{L}_{i}}$ the perimeter of the cell $i$, $J$ the strength of cell–cell adhesion, $\lambda$ the volume elasticity, ${{a}_{i}}$ the area of the cell $i$ and $s$ the spring constant use to model the coupling of apical and basal surface of the cell. The lateral displacement $\delta {\bf r}(t)$ is defined as the $xy$ projection of the vector that connects the apical (${\bf r}_{i}^{{\rm a}}$) and basal (${\bf r}_{i}^{{\rm b}}$) centroids of each cell (figure 2(a)). The first term in equation (1) accounts for the contractility of cell–cell junctions while the second one for adhesion between cells. The third term accounts for the volume conservation of cells and constraints the fluctuations of the apical area at constant height (Coburn et al 2016). The fourth term is responsible for the apical/basal crosstalk. As is described in Coburn et al (2016), the cell cytoskeleton serves to limit the skew of the columnar cell shape or the lateral displacement $(\left| \delta {\bf r} \right|)$. Thus, to account for the crosstalk between both adhesion systems we assume that the cytoskeleton functions like a spring (with a spring constant $s$) that controls the value of $|\delta {{{\bf r}}_{i}}|=\left| {\bf r}_{i}^{{\rm a}}-{\bf r}_{i}^{{\rm b}} \right|$, always opposing its increase (figure 2(a)).

Zoom In Zoom Out Reset image size

As we discussed previously (Noppe et al 2015), further analysis of the Hamiltonian in equation (1) could be performed by rearranging it as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{E}_{i}}=K{{\left({{L}_{i}}-\frac{J}{2K} \right)}^{2}}+\lambda {{({{a}_{i}}-{{a}_{{\rm p}}})}^{2}}+s{{\left| \delta {\bf r}(t)~ \right|}^{2}}-\frac{{{J}^{2}}}{4{{K}^{2}}}.\nonumber \end{align} \tag{ 2 }$$

Equation (2) shows that the $K{{\left({{L}_{i}}-\frac{J}{2K} \right)}^{2}}$ term effectively provides a relaxation of cell–cell junctions to a preferred junction length ${{L}_{i}}=\frac{J}{2K}$ so that the ratio $J/2K$ determines the preferred boundary length of cells for a monolayer packed at a density $1/{{a}_{{\rm p}}}$. In this presentation of the energy, it becomes clear that when the ratio $J/2K$ is below the minimum permissible boundary length, which corresponds to the packing of regular hexagons, the system solidifies (transits into the hard regime (Farhadifar et al 2007, Magno et al 2015, Noppe et al 2015, Coburn et al 2016)). Note, however, that in simulations we use the first form (equation (1)) of the energy function.

We account for the cell–substrate interaction by introducing cellular protrusions that appear on the basal surface of the cells (figure 2(b)). Cells now have independent basal and apical centroids and these protrusions spread from the basal centroid. These basal protrusions retract in the event of contact with the protrusions of neighboring cells and their shape is used to determine the traction a cell will gain from the substrate (figure 2(b), Coburn et al (2013, 2016)). Thus, at the beginning of a simulation, a circular protrusion distribution about the cell basal centroid is assigned to each cell in cylindrical co-ordinates as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{0}}\left(\theta \right)={{A}_{1}}\nonumber \end{align} \tag{ 3 }$$

where ${{A}_{1}}$ is the initial protrusion radius. In the event of an overlap—of any size—of the cellular protrusions of neighboring cells, the length of the offending portions of the protrusion contours are retracted (in equal proportions for both cells) towards the basal centroid by reducing their radial length by 1%. This is repeated until the overlap is finally removed in a process based on CIL (figure 2(b), Roycroft and Mayor et al (2015) and Coburn et al (2013, 2016)).

In addition to this mechanism, we incorporated cellular traction forces based on protrusion shape and a mechanism for the regrowth of lost protrusions. In this model, cellular protrusions impart a net force F_i on the cell in the direction of their growth, whose magnitude is proportional to their length (Caballero et al 2014) according to:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{{\bf F}}_{i}}(t)=\int_{0}^{2\pi }{{{P}_{i}}(\theta,t){{{\bf n}}_{\theta }}{\rm d}\theta }\nonumber \end{align} \tag{ 4 }$$

where ${{{\bf n}}_{\theta }}$ is the unit vector in direction θ. In figure 2(b), one can see that after the overlap is removed the protrusion contour develops an asymmetry and as a result we see non-zero traction vectors for both cells appearing. The cell centroid position is then updated using:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}{\bf r}_{i}^{{\rm b}}(t)}{{\rm d}t}={{h}_{0}}\left[ {{{\bf F}}_{i}}(t)-cs\delta {\bf r}(t) \right]\nonumber \end{align} \tag{ 5 }$$

where s is again the spring constant that accounts for the crosstalk between the apical and basal areas and h₀ is a parameter related to the ability of cells to generate force on the substrate that allows them to migrate. In equation (5) we also introduce the constant c which is a factor that scales the spring constant used in the CIL model. This new constant is needed to match the times scales obtained from the Monte Carlo CPM algorithm to the ones from the CIL model.

In each simulation time step the cell shape relaxes from its current shape to the target shape (equation (3)) as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial {{P}_{i}}\left(\theta,t \right)}{\partial t}=-\alpha \left[ {{P}_{i}}\left(\theta,t \right)-{{P}_{0}}\left(\theta \right) \right]\nonumber \end{align} \tag{ 6 }$$

where $\alpha$ determines the rate of regrowth and ${{P}_{i}}(\theta,t)$ is the instantaneous cell shape of the ith cell.

The simulation box size is chosen to make the area of the domain equal $N{{a}_{{\rm p}}}$ so that cells will exactly fill the box thus eliminating internal cell pressure. Apical cell configurations are sampled using a Monte Carlo algorithm in which, a Monte Carlo attempt (MCA) consists of randomly selecting a pixel at the boundary of cells and changing its spin value (cell ID) to that of its neighbor. The probability (Prob) of accepting this change is set by the energy change of the entire system $\Delta E$ caused by this spin change according to the Metropolis procedure (Graner and Glazier 1992),

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm Prob}\left(i\left(x,y \right)\to i\left({{x}^{\prime}},{{y}^{\prime}} \right) \right)=\left\{\begin{array}{@{}cccccccccccccccccccc@{}} {{{\rm e}}^{-\Delta E/\eta }},\quad \Delta E>0 \nonumber \\ 1,\quad \quad \quad \Delta E\leqslant 0 \nonumber \end{array} \right.\nonumber \end{align} \tag{ 7 }$$

where $\newcommand{\e}{{\rm e}} \eta$ is a temperature-like parameter that determines the noise in the system and which was set the same ($\newcommand{\e}{{\rm e}} \eta = ~ 3\times {{10}^{3}}$) for all simulations. This value was determined by performing a set of simulations of cell tissues in the hard regime in which the value of the temperature was raised until the system could escape from its initial local minimum energy configuration and when fluctuations of boundary length and of the cell area could be observed (data not shown). In our simulations, we also define a Monte Carlo Cycle (MCC) to be $n_v\,{\rm MCA}$ where $n_v$ is the number of boundary pixels in the system. Although in principle, selecting $n_v$ attempts for a MCC suggests that a system with longer cell boundaries will lead to more attempted MCAs and possibly slower dynamics, we found overall this does not influence the temporal scales of collective responses in our simulations since, for example, cells in the soft regime that have longer perimeters have faster collective responses (see, for example, figure 4).

In addition, to determine the boundary between the hard and soft regimes (Farhadifar et al 2007, Magno et al 2015, Noppe et al 2015, Coburn et al 2016) we performed simulations of cells in the absence of cell–substrate adhesion in the $(J,~K)$ phase space. For each simulation, we determined the average ratio $\langle \frac{{{L}_{{\rm p}}}}{{{L}_{{\rm min}}}}\rangle$, where ${{L}_{{\rm p}}}$ is the boundary length of the cells and ${{L}_{{\rm min}}}$ is the minimum possible boundary length for cells with the given cell-density which corresponds to the packing of regular hexagons. In practice, the criterion $\langle \frac{{{L}_{{\rm p}}}}{{{L}_{{\rm min}}}}\rangle$  >  1.05 was used to classify the system to be in the soft regime. Supplementary figure 1(b) shows a surface plot of $\langle \frac{{{L}_{{\rm p}}}}{{{L}_{{\rm min}}}}\rangle$ in the $(J,~K)$ phase space indicating the boundary of the hard and soft regime. For all subsequent simulations, we fix $J=5875$ and vary $K$ to obtain a set of parameters that cross the boundary between the hard and soft regime. Finally, to tune the parameter s in equations (1) and (4) and, h₀, in equation (3), to the experimental data we carried out calibration experiments by live cell 3D-imaging of small epithelial islands (~20 cells, labeled as above) using confocal microscopy (supplementary figure 2(a)). These cells were imaged and the positions of the apical and basal centroids of the bulk cells were recorded. We then calculated the horizontal projection of the distance between apical and basal centroids, the offset distance $\left| \delta {{{{{\bf r}}}}_{i}} \right|$ for each cell, and normalized this distance with the square root of average area of bulk cells. These results where then collected into a histogram (supplementary figure 2(b)) and we found that the apical and basal centroids of the bulk cells have on average an offset with a median value of ca. 10% of the mean cell diameter. We used this observation to tune the values of the parameters s and h₀, so the simulations of cells in the hard regime showed a similar average value for $\left| \delta {{{{{\bf r}}}}_{i}} \right|$. Thus, we used $s=1.0~\times {{10}^{3}}$, $c=5\times {{10}^{-5}},\alpha =2\times {{10}^{-4}}~$ and ${{h}_{0}}=5~\times {{10}^{-3}}$ in all simulations unless otherwise stated. The remaining value of the parameters used in the simulations were: $N=144,R=40,{{A}_{1}}~=~28,K~=~30~$ and $\lambda =1$ as summarized in supplementary table 1.

Several models have been developed previously to analyze the local rearrangement of epithelial cells in response to an injury. In particular, computational simulations of apical adherens junctions using vertex or CPMs have been useful to describe the mechanical behavior of cells that exert low level of traction forces on the substrate (Rozbicki et al 2015) and the possible tissue response scenarios in the case of injury (Kuipers et al 2014, Noppe et al 2015). Thus, we first investigated in silico the response of the tissue to an injury for monolayers in the CPM (Apical) model. In our simulations, monolayers were first allowed to equilibrate (i.e. achieve the minimum of the total energy) and then a local injury was introduced in a group of ten cells. In experiments, irradiation of cells by a laser causes apoptotic cell death and its posterior extrusion of the monolayer; described in detail in Kuipers et al (2014) and Michael et al (2016). Under these circumstances, dying cells tend to no longer conserve their volume (Saias et al 2015), maintain their cell–cell interactions (Andrade and Rosenblatt 2011) as well as exert forces and adhere to the substrate (Kocgozlu et al 2016), and thus no longer inhibit the formation of protrusions in the neighboring cells. To imitate cell death or injury in our simulations, we effectively removed them from the monolayer by changing the cell energy parameters ($K,J,\lambda,{{a}_{{\rm p}}},s$) to zero, a process comparable to the removal of dead cells from the epithelial layer that is observed experimentally (Rosenblatt et al 2001, Kuipers et al 2014, Lubkov and Bar-Sagi 2014). As we described before and showed in supplementary figure 2, in the absence of injury this type of modeling leads to steady-state epithelial monolayers in either a hard or a soft regime, in which cells exhibit high or low junctional tension, respectively (Magno et al 2015, Noppe et al 2015, Coburn et al 2016). Therefore, we investigated the response to injury in these two regimes.

Figure 3 shows results from simulations of the injury response of epithelial cells in the soft and hard regimes. For this series, we kept the adhesion strength fixed ($J=5875$) and varied only the contractility. We found that for high values of the junctional contractility term ($K\geqslant 50$, hard regime) the injury area starts growing right after the injury. Therefore, this recoil is the product of a local mechanical relaxation caused by the loss of tension on the side of damaged cells (figures 3(a) and (b)). After approximately 100 000 MCC all these curves reach a plateau and show only moderate growth in injury size. The plateau area scales with the junctional contractility parameter ($K$) consistent with the idea that the tissue recoil after injury is related to a mechanical relaxation (note however that the time at which the open area reaches equilibrium increases for higher K values, which reflect a 'high viscosity' of the tissue for these values of high contractility). In contrast, for low values of junctional contractility ($K=30$, soft regime), we observed an immediate decrease—yet not complete healing—of the area occupied by the dying cells after injury (figures 3(a) and (b)). These results suggest that in the absence of cell–substrate adhesion (and CIL) the relative contributions of the contractility ($K$) and adhesion parameters ($J$) determines whether the injury area will 'open' (until it reaches the plateau area) or partially 'close', which agrees with previous computational analysis of wound healing (Kuipers et al 2014, Noppe et al 2015).

We then looked for the presence of collective response in the model by measuring aspect ratios and cell orientation in either soft (K  =  30) or hard regimes (K  ⩾  40). We found in these simulations that during the phase of expansion (hard regime) or contraction (soft regime) of the injured area, cells exhibit a change in their aspect ratio. In the case of soft cells (we take K  =  30 and 40 for reference), cells immediately adjacent (row 1) to the wound site show an increase in their aspect ratio, which then becomes negative (more pronounced when K  =  40, 50 and 60). This is replicated in cells with the highest values of junctional contractility (K  =  70 and 80), where there is an increase in aspect ratio for cells located in the first row but this is less than for cells in the soft regime (K  =  30, 40). Again, cells in the hard regime and located in rows two and higher exhibit a reduction in aspect ratio (figure 3(c)). In addition, we found that the cells next to the injury (row  =  1) reorient their principal axis in the direction of the injury as shown in figure 1(g). This change is propagated up to four cell diameters in the case of K  =  80 and decreases in intensity for lower values of K becoming positive for K  =  30. Altogether, this data shows that the first row of cells re-orient almost independently of their capacity to contract their junctions—i.e. possibly just by reacting to the 'free' space—but higher junctional contractility between healthy cells when compared to junctions between healthy and damaged cells contributes to the reorientation of this first row of cells caused by the injury. Yet, although hard cells can reorient collectively, they cannot easily elongate/ reshape (i.e. increase their aspect ratio, 3c), which is needed for proper injury repair (figure 3(a)).

Although the model that only incorporates cell–cell adhesion and junctional contractility predicts partial healing of the epithelial monolayer, it is also clear that it has two major limitations when compared to experimental results: (1) it has a limited capacity to predict collective morphological changes in the hard regime in regions that surround the injury site, such as the elongation of cells in the direction of the injury; (2) it predicts partial healing of the tissue in a range of parameters that do not agree with experimental observations, i.e. this partial healing occurs in the soft regime only, when cells do not exhibit junctional tension.

The results from the previous section motivated us to investigate how the cells' traction on the substrate and CIL (Stramer and Mayor 2016) contribute to epithelial collective responses in the absence of cell–cell adhesion. Therefore, we performed a simulation of injury in monolayers that lacked cell–cell adhesions and junctional tension (figure 4) and analyzed the responses of the neighboring cells to injury. We found that loss of adhesion to the substrate of the dying cells stimulates the cells next to the injury site to extend their protrusions into the injured area (figure 4(a)). This allowed a gain of net axial orientation of neighboring cells in the radial direction to the injury site, which permitted these cells to gain traction, migrate and heal the tissue in silico (figure 4(b)). As expected, we found that the wound closes faster for cells with a higher motility parameter (h₀, figure 4(b)).

One interesting property we found for this type of in silico tissue, compared to the one with only cell–cell adhesion and no adhesion to the substrate (figure 3), is that the CIL interaction between the cells next to the injury site guides them to migrate into the free space. As they migrate, free space opens behind these cells at the basal level, which induces the axial orientation of cells in the next row, thus generating a simple mode of collective cell migration that permeates many rows (ca. 3–4 cell diameters) away from the site of injury (figure 4(c)). This observation is in line with measurements of the changes in aspect ratio of cells and the orientation of cells in the direction of the injury (figures 4(d) and (e), Brugues et al (2014)). Thus, and very simplistically, the CIL mechanism allows cells to collectively respond to the presence of injury. Since within this model wounds closed in all cases, we chose to take statistics of collective responses in the system as average values observed during the time window that spans 0% to 70% wound closure (figure 4) and between 85% and 95% of closure (supplementary figure 3). We noticed that in this model, rapid regrowth of protrusions leads to a reduction in collective response when closure is almost complete (supplementary figure 3).

To test how CIL and adhesion to the substrate may contribute to the collective responses to injury of epithelial monolayers in the hard regime we then combined the above two models: (1) the CPM, to describe cell–cell adhesion and junctional tension and (2) the CIL model, to describe cell substrate anchoring, cell migration and the interaction of cellular protrusions via CIL. In particular, we were interested in analyzing the reaction of cells that are in the hard regime (figure 5(a)), in which they exhibit significant amounts of junctional tension (figure 1(b)), and whose response to injury was not reproduced by the previous reduced models (figures 3 and 4). Figure 5(b) gives the area occupied by dying cells normalized to its pre-injury value versus time for a range of values of the junctional contractility parameter $K$. Strikingly, we found that even for injuries made on monolayers in the hard regime, the apical area occupied by dying cells shows an initial expansion followed by a slow shrinkage, which matches precisely observations from our experiments (figure 1). These results were obtained by using a high value of the apical/basal mechanical crosstalk (i.e. intracellular stiffness) parameter s (s  =  1000) that couples the basal and apical area of the in silico cells in this model. As before, this initial expansion scales with the junctional contractility (K parameter, figure 5(b)), showing that it occurs as a result of the loss of tissue tension due to the removal of the dying cells. Moreover, after this local relaxation of the tissue, the basal protrusions that form in the cells bordering the injured area stop tissue recoil and initiate the healing process by starting to pull surrounding cells into the injury area. Notably, we also found that increasing the parameter controlling the cells velocity (h₀) and/or their capacity to regrow their protrusions in the basal region ($\alpha$), increases the speed of wound closure (supplementary figure 4). Thus, cell–substrate adhesion and CIL in this model favors epithelial repair even for monolayers of cells with high levels of junctional tension.

We then characterized the collective responses to injury within this model in the hard regime (K  =  40, since it better matches the biomechanics of non-injured monolayers), with various levels of intracellular stiffness parameter s. We found that just before the healing (90% wound closure) cells in contact with the injury site extend cellular protrusions and orientate their basal area in the direction of the injury as follows from the aspect ratio and angle with respect to the direction of injury (figures 5(d) and (f), Basal), and these changes were more pronounced for higher values of intracellular stiffness. On the contrary, in the apical region we observed a different response in which increasing intracellular stiffness prevents cell shape changes in the apical region (figures 5(c) and (e), Apical). Moreover, increasing s, also leads to a reduction in the value of $\Delta \left({\rm comp}~\delta r \right)$ meaning that higher stiffness prevents the establishment of the off-set between apical and basal regions, an index of cryptic lamellipodia formation within the monolayer (figure 5(g)).

These results suggest that while coupling apical and basal layers in the model led to an accurate description of the wound area versus time response and the CIL model accounted for morphometric changes in the basal region, this was not sufficient for a proper description of the morphometric changes that we observed on the apical region of cells as the dying cells are being extruded from the monolayer. These could be better explained by either considering that in response to injury the cells surrounding the injury site soften their junctions and/or increase contractility in the interface of dying/healthy cells (i.e. purse string mechanism).

MCF-7 cells were from ATCC and cultured in DMEM; supplemented with 10% fetal bovine serum (FBS), 1% non-essential amino acids, 1% L-glutamine, 100 U ml⁻¹ penicillin and 100 U ml⁻¹ streptomycin. Cells were infected with lentivirus expressing mCherry–K-Ras^C14 (mCherry-MT) and Hist2b-GFP (NLS-GFP) or MRLC-GFP as it was described in Leerberg et al (2014), Wu et al (2014)and Michael et al (2016). Cells were then isolated by fluorescent activated cell sorted and subsequently maintained in DMEM+10% FBS plus antibiotics as described above.

Laser microirradiation was performed as described in Michael et al (2016). Confluent cells stably expressing a plasma membrane targeted mCherry construct (mCherry–K-Ras^C14, (Wu et al 2014). A 354  ×  354 µm region was imaged at 90 s intervals for about 5 h with a total of 11 z-slices (1 µm thick). To induce cell injury, a circular 85 µm diameter (ca. ten cells) centered in the field of view was irradiated at the vertical center of the monolayer (z-slice position six) after one frame of starting the acquisition. Injury was carried out using 35% transmission of the 790 nm laser for 35 iterations.

The use of laser ablation technique to assess junctional tension has been described in detail previously (Liang et al 2016). Briefly, cells stably expressing E-cadherin-GFP in an E-cadherin shRNA knockdown background (Smutny et al 2011, Priya and Gomez 2013) were used to identify the apical region of cell–cell contacts. These experiments were carried out at 37 °C on a Zeiss LSM710 system (63×, 1.4NA Plan Apo objective) using a 488 nm laser for time lapse imaging (GFP and DIC imaging) and a MaiTai (Coherent) laser set at 28% transmission and 790 nm for ablation. Time lapse imaging for about 3 min (20 frames) of a 90 µm  ×  90 µm region was taken at 2 s intervals and ablation was done after the second frame of acquisition on a 2 µm diameter circular region on the apical cell–cell junctions. Values shown average recoil curves for 15 ablated contacts.

Injury area over time

Collective cell rearrangements

Analysis of collective cell rearrangements in silico

The codes used for numerical simulations are available upon request.