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

Epithelial tissues form physically integrated barriers against the external environment protecting organs from infection and invasion. Within each tissue, epithelial cells respond to different challenges that can potentially compromise tissue integrity. In particular, cells collectively respond to injuries by reorganizing their cell–cell junctions and migrating directionally towards the sites of damage. Notwithstanding, the mechanisms that drive collective responses in epithelial aggregates remain poorly understood. In this work, we develop a minimal mechanistic model that is able to capture the essential features of epithelial collective responses to injuries. We show that a model that integrates the mechanics of cells at the cell–cell and cell–substrate interfaces as well as contact inhibition of locomotion (CIL) correctly predicts two key properties of epithelial response to injury as: (1) local relaxation of the tissue and (2) collective reorganization involving the extension of cryptic lamellipodia that extend, on average, up to 3 cell diameters from the site of injury and morphometric changes in the basal regions. Our model also suggests that active responses (like the actomyosin purse string and softening of cell–cell junctions) are needed to drive morphometric changes in the apical region. Therefore, our results highlight the importance of the crosstalk between junctional biomechanics, cell substrate adhesion, and CIL, as well as active responses, in guiding the collective rearrangements that are required to preserve the epithelial barrier in response to injury.


L Coburn et al
, which revealed that neighboring cells up to 5 cell diameters from the site of injury collectively contribute to epithelial healing (Antunes et al 2013, Lubkov andBar-Sagi 2014). However, cells play different roles within this neighborhood, depending on their location. In particular, the cells that immediately border the injured site can form an intra-cellular 'purse string' of acto-myosin that encircles the wound and which, upon contraction, leads to a closure of the wound site (or extrusion of dying cells, Martin and Lewis 1992, Bement et al 1993, Abreu-Blanco et al 2012, Antunes et al 2013, Fernandez-Gonzalez and Zallen 2013. However, the 'purse string' mechanism is believed to be most important for small wounds only. As the area to be repaired becomes larger, surrounding cells shift their behavior and extend more lamellipodia to migrate into the wound area, a process named 'lamellipodial crawling' (Abreu-Blanco et al 2012). Similarly, a transition from 'purse string' to 'lamellipodial crawling' type of behavior has been found during epithelial cell extrusion at different monolayer densities (Kocgozlu et al 2016) suggesting that, in general, changes in epithelial mechanics (larger wounds or crowded epithelia) determine the type of behavior that cells use to preserve the barrier integrity.
Irrespective of the precise repair mechanism, epithelial cells surrounding the wound must collectively remodel their cell-cell junctions and change their morphology as well as reorient themselves in the direction of the injury in order to minimize friction during their migration towards the damaged area (Tambe et al 2011). Although different mechanisms have been proposed to underlie this junctional remodeling (Antunes et al 2013, Hunter et al 2015 the biomechanical principles that govern epithelial collective responses are still not fully understood. To address this question, we performed a quantitative morphometric analysis of collective responses to injury and attempted to reproduce these changes using different minimal models of epithelial monolayers in which cells adhere to one another and/or to the substrate as well as exhibit contact inhibition of locomotion (CIL, Coburn et al 2016). By matching the results of different models with experimental data we provide quantitative and mechanistic insights into the forces that drive epithelial collective behavior in response to injury.

Experimental results: epithelial collective cell responses to laser induced-injury
Collective cell behavior is an intrinsic property exhibited by epithelial cells in different contexts such as morphogenesis, wound healing and cancer invasion (Friedl and Mayor 2017). To test the ability of our model to reproduce the collective behavior, we first performed quantitative experiments and then used their results to evaluate the predictions of the different modeling approaches. Several experimental techniques have been previously employed to analyze collective cell behavior in response to injury. Of these, laser micro-irradiation is readily implemented and allows changes in cell morphology and movement to be monitored with high spatiotemporal resolution (Abreu-Blanco et al 2012, Antunes et al 2013). Also, laser micro-irradiation causes apoptotic cell death (Kuipers et al 2014, which removes the cell's capacity to conserve its volume (Saias et al 2015), maintain its cell-cell interactions (Andrade and Rosenblatt 2011) as well as exert forces on and adhere to the substrate (Kocgozlu et al 2016) As such, these cells no longer inhibit the protrusions (i.e. lamellipodia formation) of neighboring cells. Thus, we grew to confluence, epithelial MCF-7 cells stably co-expressing a membrane targeted mCherry (MT-mCherry) and nuclearlocalized-GFP (NLS-GFP) and then introduced an injury by laser micro-irradiation (supplementary movie 1 (stacks.iop.org/PhysBio/15/024001/mmedia)) and characterized quantitatively the dynamic changes of epithelial cell shape and collective behavior in response to injury. We found that immediately after an injury there was, on average, a small immediate increase in the apical area occupied by the ca. 10 cells that were damaged (figures 1(a) and (b)), being evident some variations between individual experiments (supplementary figure 1(a)). This initial expansion of the injured area is probably due to a local mechanical relaxation of the tissue since MCF-7 cells exhibit a significant amount of junctional tension (figure 1(c)). After this local relaxation, the injured cells (Andrade andRosenblatt 2011, Michael et al 2016) started to extrude and at the same time neighboring cells started to extend lamellipodia and occupy the substrate previously occupied by the damaged cells before micro-irradiation (figure 1(a), supplementary movie 1). In addition, we found that as the healing process progressed, epithelial cells immediately adjacent to, and few cell diameters beyond, the site of damage became more elongated with their major axis orientated in the direction of the injury (figure 1(a), t = 2 h).
To quantitatively assess these changes we then measured different shape descriptors and cell orientations in the direction of the injury as a function of their distance (in cell rows) to the injury site (figure 1(d)). First we measured the aspect ratio of cells in the apical region where the zonula adherens of epithelial MCF-7 cells is localized  and in the basal region corresponding to the lowest part of the cell-cell interface. Note that due to technical limitations we did not measure the cell area at the actual cell-substrate interface, as it is difficult to segment cell boundaries at this location based on the fluorescence of the MT-mCherry marker. Our results show that at time points close to healing (90% closure) the apical and basal aspect ratio of cells situated between 1 and 3 cell diameters from the injury site increased with respect to their aspect ratio at the time before injury (figure 1(e)). Epithelial collective responses to injury. (A) Confluent MCF-7 cells expressing a plasma membrane targeted mCherry (MT-Cherry) were locally injured by micro-irradiation (dashed yellow line) with a two-photon laser and the morphological responses were analyzed by 3D (z-stack) time lapse imaging. Panels show merged projections of the plane of cells in contact with the substrate (basal, pseudo-colored in magenta) and the plane of cells that contain the cell-cell junctions (apical, pseudo-colored green) at times before injury and at 0 min, 25 min and 2 h after injury. Scale bar 50 µm. (B) Quantifications of the area within the monolayer occupied by injured cells as a function of time. Plot shows the mean ± S.E.M of 3 independent experiments. (C) Average recoil of the Zonula Adherens after ablation. Experiments were performed on cells expressing E-cadherin-GFP as described in materials and methods. Data show the average recoil for 16 analyzed junctions and its S.E.M. (D) Definition of cell positions with respect to the injury centre used for quantitative analysis. (E) Changes in aspect ratio at the apical and basal planes of the cells (Δ = Final-initial, see also materials and methods) in response to injury. Changes were measured as a function of the distance of cells from the site of injury. (F) and (G) Scheme (F) and quantification (G) of the changes in relative orientation of cells (measured either at its apical or basal plane, see the materials and methods section) as a function of their distance from the site of injury. (H) and (I). Scheme (H) and quantitation (I) of the component (comp.) of the vector (δr) that defines the offset between apical (r a ) and basal (r b ) centroids in the direction of the injury. In (E), (G) and (I), p values of statistical comparison of linear dependence versus no change across the rows is reported. Results shown correspond to average values obtained from three independent experiments/ movies (details of number of cells quantitated in each row are included in the materials and methods section).
We then analyzed the relative orientation of cells with respect to the site of injury at the moments before healing (90% closure). We measured this as the angle between the major axis of an ellipse that we fit to the cell (either in the apical or basal area) and the line that defines the cell's position with respect to the injury site (figure 1(f)). We found that for the apical and basal region, this angle becomes smaller (i.e. cells are more aligned, ∆γ = γ after − γ before < 0) as the healing process proceeds. This change (absolute value) is larger for cells closer to the site of injury site and decreases to approximately zero for cells located at 4 cells diameters from the site of injury (figure 1(g)), which resembles the phenomenon of collective alignment of cell's principal axis during collective epithelial migration (Zaritsky et al 2015). Finally, we also quantitatively measured the skewness of cells in the direction of the injury. In our measurements, we refer to skewness as the relative displacement between the basal and apical centroids and this is an index of how 'tilted' the cells become as they migrate (figures 1(h) and (i) Coburn et al (2016)). This parameter, ∆ (compδr) = compδr after − compδr before , compares at some extent to the formation of cryptic lamellipodia that have been observed in collectively migrating cells (Farooqui andFenteany 2005, Trepat et al 2009). In our description (see materials and methods), a 'positive' skewness parameter means that cells extend their basal area more than the apical area in the direction of the injury whereas a 'negative' means cells extend their basal area in the direction opposite to the injury site (figure 1(h)), (ii). Our measurements of skewness show that after injury cells, on average, extend their basal area in the direction of the injury and this is observed for cells up to 3 cells diameters from the injury site (figure 1(i)). Overall, this quantitative data shows that the presence of a micro-injury causes a local relaxation in the tissue and triggers collective responses that involve the elongation of cells and orientation of their major axes in the direction of injury as well as an average extension of cryptic lamellipodia in the direction of the injury in cells located up to 3 cell diameters from the site of injury.

Computational model of epithelial cells
To investigate the driving forces of collective cellular response to injury we develop a lattice-based mechanistic model of epithelial cells based on their capacity to interact with one another and with the substrate (Coburn et al 2016). Cell-cell adhesion (apical region) is modeled using a cellular Potts model (CPM) algorithm (Graner and Glazier 1992, Kabla 2012, Noppe et al 2015, Magno et al 2015, Albert and Schwarz 2016 where a cell is made up of a number R 2 of pixels (lattice sites) that are allowed to change the index (cell ID or spin) to that of the neighboring cell, thus moving the location of the cell-cell boundary. The pixel attribution switching is performed according to a probabilistic Metropolis rule, which leads to minimization of the tissue energy function (described in detail below). The energy function includes terms related to cell-cell adhesion, apical ring contractility, and cell volume as well as cell motility (Kabla 2012). For cell-substrate adhesion we used our previously described protrusiondriven cell propulsion model that also incorporates CIL (Coburn et al 2013). Both cell-cell adhesion and cell-substrate adhesion were then coupled by an elastic spring that represents the mechanical crosstalk between these two adhesion systems and as such is a measure of the intracellular stiffness (Coburn et al 2016). This leads to a modeling approach that, although not strictly 3D, allows the analysis in simplified terms of the main properties that define collective epithelial cell behavior in response to injury.

Cellular Potts model of the apical region of cells: CPM (Apical) model
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 (x, y) 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 = √ a 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): 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, λ 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 δr(t) is defined as the xy projection of the vector that connects the apical (r a i ) and basal (r b i ) 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 (|δr|). 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 |δr i | = r a i − r b i , always opposing its increase (figure 2(a)).
As we discussed previously (Noppe et al 2015), further analysis of the Hamiltonian in equation (1) could be performed by rearranging it as: (2) Equation (2) shows that the K L i − J 2K 2 term effectively provides a relaxation of cell-cell junctions to a preferred junction length L i = J 2K so that the ratio J/2K determines the preferred boundary length of cells for a monolayer packed at a density 1/a 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.

Cell substrate adhesion with CIL model: CIL (Basal) model
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 (2013Coburn et al ( , 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 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 (2013Coburn et al ( , 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: where n θ 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: where s is again the spring constant that accounts for the crosstalk between the apical and basal areas and h 0 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: where α determines the rate of regrowth and P i (θ, t) is the instantaneous cell shape of the ith cell. To account for crosstalk between cell-cell and cell-substrate adhesion systems we assume that the cytoskeleton functions like a spring that controls the displacement between apical r a i and basal r b i centroids (B) Implementation of CIL in the model. In the first panel the protrusion of two cells are initially symmetrically distributed thus neither cell experiences a net force from traction on the substrate. There is also an overlap of the protrusion contours. In the second panel, the CIL interaction retracts protrusions in the radial direction until the overlap is removed. This results in an asymmetry in the protrusion distribution for both cells and as such both cells now gain a net force of traction from the substrate.

Simulation details and parametrization
The simulation box size is chosen to make the area of the domain equal Na 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 ∆E caused by this spin change according to the Metropolis procedure (Graner and Glazier 1992), where η is a temperature-like parameter that determines the noise in the system and which was set the same (η = 3 × 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 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 Lp Lmin , where L p is the boundary length of the cells and L 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 Lp Lmin > 1.05 was used to classify the system to be in the soft regime. Supplementary figure 1(b) shows a surface plot of Lp Lmin 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 0 , 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 |δr i | 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 0 , so the simulations of cells in the hard regime showed a similar average value for |δr i |. Thus, we used s = 1.0 × 10 3 , c = 5 × 10 −5 , α = 2 × 10 −4 and h 0 = 5 × 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 λ = 1 as summarized in supplementary table 1.

Numerical simulation results
To analyze the mechanisms responsible for different features of epithelial collective reorganization we performed numerical simulations of epithelial injury in the following three scenarios (summarized in table 1) of our model of epithelial cells: i. Using only the CPM (Apical) model to simulate monolayers of cells with only cell-cell adhesion and junctional contractility; no CIL, no random motion and no coupling to the basal surface ( figure 3).
Apical region: Phys. Biol. 15 (2018) 024001 ii. Using only the CIL (Basal) model to simulate monolayers of cells that exhibit adhesion to the substrate and interact with each other through CIL. The CPM model and therefore the coupling terms between this and CIL model are not included in these simulations (figure 4). iii. Using the CPM + CIL (Apical + Basal) model. This is a combination of the above reduced models, in which cell-cell adhesion is mechanically coupled to cell-substrate adhesion through an elastic spring of constant s ( figure 5).

Injury response in the CPM (apical) model
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, λ, a 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 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 adja- cent (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 contrac-tility (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 cellcell 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 exper imental observations, i.e. this partial healing occurs in the soft regime only, when cells do not exhibit junctional tension.

Injury response in the CIL (basal) model
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 cellcell 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 0 , figure 4(b)).
One interesting property we found for this type of in silico tissue, compared to the one with only cellcell 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).

Injury response in the CPM + CIL (apical/basal) model
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 0 ) and/or their capacity to regrow their protrusions in the basal region (α), 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 noninjured 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 ∆ (comp δr) 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).

Discussion
We study a simple model in which we combine CIL and junctional mechanics to analyze collective morphological cell responses to injury in epithelial tissues. In particular, simulations of monolayers in the hard regime have the capacity to exhibit salient properties of epithelial tissues that have been observed experimentally: (i) cells exhibit junctional tension, (ii) after injury there is a local relaxation of the tissue, (iii) the tissue is able to heal and (iv) healing involves morphological collective responses of epithelial cells in the neighborhood of the site of injury.
The modeling demonstrates that cell-substrate adhesion and CIL play crucial roles in the healing process. First, we observed that an initial expansion of the area of injury observed in simulations (and in experiments) is a consequence of the epithelial cells' ability to generate junctional tension. Loss of junctional tension caused by the injury leads to an elastic mechanical relaxation of the cells surrounding the injured area. However, the extent to which this local relaxation is propagated across the tissue is limited by the cell adhesion to the substrate, similar to what is observed in epithelial cell islands (Ng et al 2014, Coburn et al 2016. Second, cell adhesion to the substrate and CIL allow for the extension of protrusions and the migration of cells into the injury area. During this process, cells at the edge of the injured area migrate first and thus leave gaps behind them, which favors the re-orientation, asymmetry and migration of cells in the rows behind so that a collective response is developed. Notably, these collective rearrangements are not observed in simulations in which cells are only allowed to interact with each other and generate junctional contractility (figure 3) or in cells that exhibit high levels of intracellular stiffness ( figure 5(g)).
In addition, the model shows that extension of protrusions both in vivo (figure 1, supplementary movie 1) and in silico (figure 5) occurs later and on a time scale that is slower than the instantaneous elastic relaxation, being this similar for a wide range of parameter values (figure 5). This illustrates the robustness of our model for reproducing these spatiotemporal features of the epithelial response to injury.
We also note that our model provides an explanation for the length scales, at which collective behavior arises in response to injury, in particular on the basal region and with the offset distance between apical and basal cell centroids (Farooqui and Fenteany 2005, Antunes et al 2013, Lubkov and Bar-Sagi 2014. In particular, mechanical coupling between cell-cell adhesion and cell-substrate adhesion and the presence of CIL serves as a mech anism that leads to a stress gradient that propagates to cells behind the sites of injury, altering their shape and polarizing them collectively in the direction of the injury. A similar stress gradient has been reported for collective behavior in expanding epithelial islands (Banerjee et al 2015, Zimmermann et al 2016, thus suggesting that the biomechanical properties of cells profoundly affect the spatial scales of collective responses. Several active response mechanisms that drive collective behavior of cells in response to injury have been proposed in the literature, and we found these could be critical for morphometric changes on the cells' apical region. For example, neighbouring cells may soften their bonds to facilitate cell rearrangements that occur during neural tube closure and cell extrusion (Escuin et al 2015, Hashimoto et al 2015; observations that agreed with our numerical simulations of 'soft' cells (figure 5) and reduction of MRLC-GFP on some of cell-cell junctions in cells that surround the site of injury (supplementary movie 2). In addition, directional flow of actomyosin towards the injury-live cell interface to form an actomyosin 'purse string', that we also observed in our experimental conditions (supplementary movie 2), could be favored by the orientation of protrusions and changes in the cell aspect ratios (Antunes et al 2013), a notion supported by the work of the Ladoux Lab where it is shown that 'purse strings' are preferentially stabilized in the regions of negative curvature (Ravasio et al 2015). Also, we do not exclude the possibility that motility may be enhanced through mechanotransduction in the tissue next to the injury site. Cells next to the injury site gain greater traction at the substrate due to a passive resistance to migration from the cells in the next rows (Coburn et al 2016). This could in principle enhance cells' traction on the substrate through force transduction as it has been observed in other experimental systems (Weber et al 2012, Mertz et al 2013. Although all these phenomena have been shown to contribute, to some extent, to the collective response of neighboring cells to injuries, there is no clear consensus on what limits the spatiotemporal scales of this process. Our results show that the presence of CIL and the ability of cells to adhere to the substrate allow cells to extend their basal area relatively more than the apical area thus polarizing them towards the site of injury as well as limiting the space over which mechanical relaxation of the tissue occurs. This is interesting as it was recently observed that polarization of cells in the basal area limits the amount of available actin to generate actomyosin networks to sustain contractility in the apical region of cells (Lomakin et al 2015). Thus, we hypothesize that these mechanisms allow cells to transmit information about the presence and location of the injury across the tissue that can be used by other cells to generate controlled and local active responses such as proliferation (Aragona et al 2013).

Materials and methods
Cell Culture and Transfections 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 −1 penicillin and 100 U ml −1 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.

Two-photon laser induced epithelial injury
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 , . 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.

Laser ablation experiments to measure junctional tension
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 andGomez 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.

Image analysis
Injury area over time The area of injury and its variation were determined in Image J by drawing a region of interest (ROI) corresponding to the damaged area in the plane of the monolayer and following its changes over time. Data were normalized to the area occupied by cells that died upon injury before ablation. Data shown correspond to the average of three independent movies.

Collective cell rearrangements
First, cells were indexed with respect to their position (in row number) to the site of injury. Then, ROIs corresponding to the boundaries of each cell in its most apical and most basal planes were drawn using the drawing tools in Image J and added to the ROI manager. In addition, an additional ROI that describes the region of injury (or region occupied by dying cells) was drawn. After ROIs were drawn the following options were chosen in the set measurements menu in Image J (i) centroid, (ii) shape descriptors and (iii) fit ellipse. Thus, after 'measure' in Image J, we extracted the following parameters from each ROI: (a) X and Y coordinates of the centroid of each ROI, (b) the angle between the major axis of an ellipse that fits the cell boundary and the X axis of the image and (3) the aspect ratio (ratio between the major axis and minor axis of the ellipse that fits the cell boundary). Using this information, we then calculated the average value of the following quantities for each cell row from the site of injury: i. Change of aspect ratio before and after injury: we measured for each cell how much their aspect ratio changes before and after injury (t = 2 h). These measurements were done for both the basal and apical region of cells. Average values were then calculated for cells in the same row within a movie. ii. Change of cells' orientation with respect to the site injury before and after injury. For these measurements, we first calculated the vector position of a cell with respect to the site of injury using the information of the cells centroid and the centroid of the area of injury. From this information, we calculated the orientation (or slope, m 1 ) of the line that connects the site of the injury and the analyzed cell. Similarly, by fitting a cell boundary with an ellipse we calculated the orientation (or slope, m 2 ) of the major axis of a cell with respect of an image using the 'angle' values obtained from image J. Using the values of m 1 and m 2 we then calculated the acute angle (γ) between cells orientation and the direction to the injury using tanγ = m 1 − m 2 1 + m 1 m 2 and how it changes before and after (2 h) injury (i.e. ∆γ = γ after − γ before ). Average values were then calculated for cells in the same row within a movie. iii. Changes in cell skewness before and after injury. A measure on how much a cell is tilted within the monolayer is given by the 2D projection in the basal plane of the vector that connect the basal and apical centroids (δr, Coburn et al (2016)) and a measure on how much this is aligned in the direction of the injury is given by the component of this vector in the direction that connects the basal centroid of the cell and the centroid of the area of injury (compδr). Thus, we measured the magnitude of compδr using the dot product equation between these two vectors and dividing it by the distance between the injury site (centroid) and the position of the basal centroid of the analyzed cell. Average changes (∆) in compδr before and after 2 h of injury were calculated for each cell row.

Analysis of collective cell rearrangements in silico
In numerical simulations, changes (Δ) in different cell morphometric parameters (i.e. aspect ratio, cell orientation (γ) at the apical and/or basal regions and compδr were determined as follows. For a given simulation the aspect ratio and the orientation of the long axis was found for each cell on each MCC. We then found the change in aspect ratio and orientation of cells from the current MCC as compared to their values just before injury. These values from cells where then collected into centric rings, or rows, from the injury site (we extended up to five rows from the wound site) and averaged in each row. We then had values for the change in aspect ratio and orientation of the long axis in each row for all MCC. For a given simulation, we then chose a window to average these values over. In the CPM model, as wounds never closed, we chose this averaging window to be the period after wound size had equilibrated. In the CIL model this window, was the first 70% of wound closure (i.e. from when the injury was first made to when the injury was 30% open, figure 4) or during the times the 'wounds' are between 85-95% closed or, in the case when the wound does not close, from 85%-100% of equilibrated wound area (supplementary figure 3). In the CPM/CIL model we chose this window to correspond to the times the 'wounds' are between 85-95% closed or, in the case when the wound does not close, from 85%-100% of equilibrated wound area.
For each simulation, we now had average change in aspect ratio and orientation in each row over the chosen window. Simulations where then repeated multiple times (96 times in the CPM (Apical) model, 96 times in the CIL (Basal) model and 24 times in the CPM+CIL (Apical + Basal) model and the mean and SEM (plotted in figures 3-5) was calculated from there.

Statistical analysis
Linear regression was performed on different measured quantities and a statistical test was performed against the no change (Δ quantities = 0) null hypothesis of across all the rows using the nonlinear regression function using GraphPad Prism 7. For this analysis, each experimental replicate was used for non-linear regression whereas the mean values for each row were used when data from simulations were analyzed.

Code availability
The codes used for numerical simulations are available upon request.
Cancer Biology Imaging Facility, established with the generous support of the Australian Cancer Research Foundation.