Development of a New 3D Hybrid Model for Epithelia Morphogenesis

1. Introduction

Many morphogenetic events in epithelia can be successfully described with two-dimensional models. Some examples include tissue intercalation (Munjal et al., 2015), jamming transitions (Bi et al., 2016), or collective cell migration (Sunyer et al., 2016). However, monolayers are also subjected to observable out of plane deformations like tissue folding (Tozluoğlu et al., 2019), invagination (Bielmeier et al., 2016), extrusion (Deforet et al., 2014), or delamination (Eisenhoffer et al., 2012). In these cases, two-dimensional models with folding capabilities (Misra et al., 2016), or purely three-dimensional models, usually implemented in the continuum context, seem necessary to capture the underlying contractile mechanisms.

There are also other problems that, despite being studied extensively in two-dimensions, contain three-dimensional contributions that have not been included in their modeling and description. Some examples are neural crest zippering (Hashimoto et al., 2015) or wound healing in monolayers (Antunes et al., 2013; Brugués et al., 2014). The former in fact originates from a precedent tissue folding, while experimental observations of the latter indicate that wound closure has some variations along the cell apicobasal axis (Zulueta-Coarasa et al., 2014). Understanding forces in these processes requires models that are able to reproduce both junctional mechanics and deformations that are different at the apical and basal sides of the monolayer. The present paper introduces a vertex model that includes these ingredients.

The discrete nature of tissues makes vertex modeling an ideal approach that has been successfully employed to simulate cell dynamics of monolayers and study, for example, T1 transitions (Bi et al., 2015), phase transformations (Bi et al., 2016), and wound healing (Brugués et al., 2014; Staddon et al., 2018; Tetley et al., 2019). For further reference on vertex models see for instance the recent review (Alt et al., 2017).

Three-dimensional versions are more scarce, but have been also recently developed to study curved monolayers (Gómez-Gálvez et al., 2018), cyst formation in monolayers (Bielmeier et al., 2016), folding in epithelial shells (Misra et al., 2016), or general morphogenesis (Okuda et al., 2013, 2015).

We here extend a previous version of our 2D hybrid model to three dimensions (Mosaffa et al., 2018), with the ability to uncouple intercalation on the apical and basal sides of the monolayer. Differential intercalation and the geometrical definition of the cell poses special computational challenges. The analysis of shape transitions between the two sides of a monolayer have been analyzed for curved monolayers (Gómez-Gálvez et al., 2018). Differential apical and basal intercalation has also been seen in the Drosophila salivary gland–termed “interleaving” (Sánchez-Corrales et al., 2018). A similar idea is employed here in flat geometries by defining an intermediate vertex that facilitates neighbors changes.

In our model, we resort to a hybrid approach, where the cell-centers of the apical and basal layer are the main degrees of freedom. As is customary in vertex models, each cell-center is surrounded by a set of vertices, where mechanical balance is also imposed between all connected bar elements (Barton et al., 2017). Vertex positions are constrained by the triangulation of the cell-centers. By inserting these constraints into the equilibrium equations, we manage to reduce the computational size of the model, which is solely described by cell-center positions, and thus simplifies the topological definition of the monolayer. We note though that despite this reduction, our final equations still include vertex equilibrium.

Illustrative results are shown by applying the model for the analysis of wound healing, with specific contractility evolution at the wound edge. Each bar element of the model adopts a viscoelastic rheology, which is based on a dynamic change of the rest-length (Muñoz and Albo, 2013; Staddon et al., 2018), and allows calibrating the short term recoil process.

2. Methodology

2.1. Three-Dimensional Vertex Model

2.1.1. Monolayer Geometry

The cells in the monolayer are initially defined by the centers (nodes) of their surfaces at the apical (top) and basal (bottom) sides, respectively denoted by $x_A^i$ and $x_B^i$, i = 1, …, N_cells. In our simulations, these cell centers are equal at the apical and basal sides, and correspond to the measured experimental locations of the cell center of mass.

After applying a 2D Delaunay triangulation on each side, the cell apical, and basal boundaries are then constructed by joining the vertices $y^I,I=1,\dots ,N_y$ on each side. These vertices are located at the barycenters of the triangles that surround each cell center. The two 2D layers are then joined with vertical and diagonal segments that join apical and basal nodes, and also basal and apical vertices, forming prism-like polyhedra. As a result, each cell is formed by an apical and basal center (node) and two sets of vertices defining the apical and basal boundaries. Figure 1 shows the construction process and final polyhedra.

FIGURE 1

Figure 1. (A) Scheme of a cell with vertex bar elements (red) and nodal bar elements (black). (B) Adding apical vertex (red) and nodal (black) network. (C) Construction of cells with vertical bar elements. Diagonal elements in vertex network have been omitted for clarity.

The initial position of the apical and basal cell-centers (nodes) are taken from two dimensional experimental images, where cell center positions are measured. The cell boundaries are located at the barycenters of the resulting triangles computed from a Delaunay triangulation of the cell centers. Consequently, the initial cell-center locations and shape of the cell areas at the apical and basal surfaces are equal. However, they are allowed to change their connectivity independently in subsequent time-steps, forming polyhedra that may have different polygonal shapes and number of sides at each apical and basal surface. The apico-basal transition between two polygons with different number of segments is facilitated by the definition of intermediate vertices, located between the apical and basal surfaces, as shown in Figure 2.

FIGURE 2

Figure 2. (A) Scheme of a cell with vertex bar elements (red) and nodal bar elements (black). (B) Definition of intermediate vertices in order to ease different connectivity at basal surface (B1) and apical surface (B3).

More specifically, the positions of the vertices at the apical and basal surfaces are computed from the interpolation of the corresponding cell-centers by using the following constraint equation,

$\begin{array}{ll}{y}^I={\displaystyle \sum _{i=1}^3}{N}^i({\xi }^I)x_I^{i}s,I=1,\dots ,N_{tri}& (1)\end{array}$

where Nⁱ(ξ^I) are interpolation functions associated to each one of the nodes $x_I^1,x_I^2,x_I^3$ forming triangle I, and ξⁱ a parametric coordinate. Hereafter we set Nⁱ(ξ^I) = 1/3, so that vertices are located at the triangle barycenters. As a consequence, the model is a hybrid version between a cell-centered and a purely vertex model. Vertices are used for defining the cell boundary (cortex), but cell-centers are kept as degrees of freedom (DOF) describing the whole cell kinematics. Similar approaches can be found in Barton et al. (2017) and Mosaffa et al. (2018).

The resulting geometrical construction is formed by two coupled networks: the set of segments ij joining nodes nodal network, and the segments IJ joining the vertices vertex network. The motivation of the hybrid approach is threefold: (i) it reduces the number of DOF, and thus the size of the resulting system of equations to solve at each time-step, (ii) it allows us to model mechanical interactions between cell-centers (nodal network of triangles) and joint mechanics (vertex network of polygonal cell boundaries), and (iii) it also provides a mechanical coupling between the two networks.

In order to define the vertices at the boundary of the patch, we add a set of boundary nodes which are not the center of any cell. We note also that intermediate vertices are free to move as independent DOF, so that they are not interpolated according to Equation (1). Furthermore, we will also relax the constraint in Equation (1) for those vertices that are at the edge of the wound, and let all those free vertices, denoted by $y_w^I$, to be additional DOFs. As a result, the kinematics of the monolayer is fully defined by the cell centers positions xⁱ and the coordinates of the relaxed vertices $y_w^I$. We note that this relaxation is introduced in order to avoid zig-zag effects at the wound edge. The reader is referred to Mosaffa et al. (2018) for further analysis on this effect.

2.1.2. Mechanical Equilibrium

In each one of the nodal and vertex network, we distinguish apical, basal and lateral segments. The total energy of the system W(x, y_w) is defined by the sum of nodal and vertex elastic contributions, and a volume penalization as,

$\begin{array}{ll}W(x,y_w)=W_N(x)+W_V(y(x),y_w)+W_{Vol}(y(x),y_w)& (2)\end{array}$

Each term is defined by,

$\begin{array}{l}{W}_N(x)=\frac{k_N}{2}{\displaystyle \sum _{ij}(}{l}^{ij}-L^{ij}{)}^2\\ W_V\left(y(x),y_w\right)=\frac{k_V}{2}{\displaystyle \sum _{IJ}{(l^{IJ}-L^{IJ})}^2}\\ W_{Vol}\left(y(x),y_w\right)=\frac{{\lambda }_{Vol}}{2}{\displaystyle \sum _{i=1}^{N_{cells}}{\left(\frac{V^i-V_0^i}{V{i}_0}\right)}^2}\end{array}$

The material parameters k_N and k_V measure the stiffness of the nodal and vertex segments, respectively. The measures l^ij = ||xⁱ − x^j|| and l^IJ = ||y^I − y^J|| correspond to the current observable lengths of each nodal and vertex segment. The rest-lengths L^ij and L^IJ are internal variables, not necessarily constant. Their evolution will furnish viscous properties to the monolayer, and will be defined in section 2.1.3. The penalization parameter λ_Vol is set to a value such that relative volume difference $|V^i-V_0^i|/V_0^i$ at each cell i is kept between 5 and 10%, as experimentally measured (Gelbart et al., 2012).

Mechanical equilibrium is achieved by minimizing the total energy with respect to the positions of the nodes and the relaxed vertices, i.e.,

$\begin{array}{ll}\frac{\partial W_N}{\partial x^i}+\frac{W_V}{\partial y^I}\frac{\partial y^I}{\partial x^i}+\frac{\partial W_{Vol}}{\partial y^I}\frac{\partial y^I}{\partial x^i}=0,i=1,\dots ,N_{nodes}& (3)\end{array}$

$\begin{array}{ll}\frac{W_V}{\partial y_w^I}+\frac{\partial W_{Vol}}{\partial y_w^I}=0,I=1,\dots ,N_{relax}& (4)\end{array}$

The terms $\frac{\partial y^I}{\partial x^i}$ are computed from the constraint equation in (1) as, $\frac{\partial y^I}{\partial x^i}=N^i({\xi }^I)\text{I}=1/3\text{I}$, with $\text{I}=\frac{\partial x^i}{\partial x^i}$ the three-dimensional identity matrix. Due to the non-linearity of the equations in (3), the solution is found using a Newton-Raphson strategy enhanced with line-search strategies.

After each converged step at time t_n, with connectivity C_n and nodal and vertex positions $x_n^i$ and $y_w^I$, we compute new nodal and vertex positions, $x_{n+1}^i$ and $y_{w,n+1}^I$, respectively, and a new connectivity C_n+1 using a modified Delaunay triangulation of the apical and basal surfaces. New triangles are formed if their aspect ratio at time t_n+1, denoted by r_n+1, is improved according to the relation

$r_{n+1}<(1+to{l}_r)r_n$

with tol_r a non-negative numerical parameter. When tol_r = 0, the standard Delaunay algorithm is recovered, while for tol_r > 0, suboptimal stretched triangles and cells are permitted.

2.1.3. Rheological Model

Each vertex segment in the model has the ability to respond according to a viscoelastic rheological law. This law is implemented by resorting to a variable rest-length L^ij or L^IJ, which evolves with the following equation:

$\begin{array}{ll}\frac{1}{L}\frac{dL}{dt}=\gamma \left(\frac{l-L}{L}-{\epsilon }^c\right)& (5)\end{array}$

This law reflects the fact that as far as the strain measure (l − L)/L is different from a contractility ε^c, the rest-length will evolve, and that L will remain unchanged when the strain reaches the value ε^c. It has been proved that for ε^c = 0, such evolution law gives a similar response to a Maxwell viscous model (Muñoz and Albo, 2013). The intrinsic contractility ε^c has been included in order to mimic the contractile state of cells (Khalilgharibi et al., 2019; Wyatt et al., 2020).

In addition, we also consider the local actin concentration at the wound edge. This is implemented through an additional contractility ${\widehat{{\rm Y}}}^c$, which increases the tension at the elastic branch of the vertex segment as

${\sigma }_V=k_{V0}({\epsilon }^e+{{\rm Y}}^c)$

We will use time-varying values of ${\widehat{{\rm Y}}}^c$ at the apical and lateral sides of the cells, as it is explained in section 2.2. As a result, an additional tension equal to $k_{V0}{\widehat{{\rm Y}}}^c$ is being applied on those vertex segments.

Figure 3 depicts the rheological model employed at the nodal and vertex segment. The purely elastic behavior is obtained by the minimization of energies W_N and W_V using, respectively, stiffnesses k_N and k_V0, while the viscous response results from the implementation of the evolution law in (5) in an additional energy term W_V with stiffness k_V. For simplicity and to avoid having to fit too many parameters, we have set γ = 0 for the nodal elements, which yields a purely elastic behavior, as depicted in Figure 3. The calibration of the material parameters will be explained in section 3.1.

FIGURE 3

Figure 3. Rheological model for nodal (left) and vertex (right) segments. Contractility parameters ${{\rm Y}}_A^c$ and ${{\rm Y}}_L^c$ are only used at the wound edge.

2.1.4. Numerical Solution

The set of non-linear equations in (3) depends non-linearly on the nodal positions xⁱ and on the relaxed vertices $y_w^I$. In addition, these equations, which may be expressed as,

$\begin{array}{ll}g(x,y_w)=0& (6)\end{array}$

also include the rest-lengths L of each bar element. These rest-lengths obey the differential equation in (5), which is discretized in time using a θ − weighted scheme,

$\begin{array}{ll}{L}_{n+1}-Ln=\Delta t{L}_{n+\theta }\gamma \left(\frac{l_{n+\theta }-L_{n+\theta }}{L_{n+\theta }}-{\epsilon }^c\right)& (7)\end{array}$

with (•)_n+θ = (1 − θ)(•)_n + θ(•)_n+1. We used the value θ = 0.5, which yields a second-order accurate and unconditionally stable scheme in linear systems. The relation in (7) allows obtaining an expression of L_n+1 as a function of l_n+1 and other values at time t_n. This expression is inserted in the system of equations in (6) at each time t_n+1, and solved in an implicit manner with a fully Newton-Raphon iterative process. We set the convergence tolerance tol = 1E−10 and impose the convergence condition ||δx|| < tol and ||g(x, y_w)|| < tol. We have used the time-step size Δt₁ = 0.6 min during the recoil process (up to t = 6 min), and Δt₂ = 1 min during the closure process. In some simulations, the latter time-step was halved in some of the increment in order to achieve convergence. The total simulation when using a patch of 205 cells took around 150 increments, and a run time of approximately 20 min in Matlab R2018a in a Windows machine with Intel(R) Core(TM) i7-6700 CPU @ 3.4 GHz, and 16 GB RAM memory, working with 2 processors.

2.2. Wounding and Contractility

The Drosophila larval wing imaginal disc is a pseudostratified epithelium containing highly columnar cells. This tissue is therefore an ideal experimental system for investigating wound healing cell behaviors in 3-dimensions. We wounded wing disc epithelia by ablating multiple tricellular junctions at the level of adherens junctions using a pulsed TiSa laser (Tetley et al., 2019). We then imaged the wound healing response in wing discs expressing a GFP tagged form of non-muscle Myosin II using 3D time-lapse confocal microscopy with intervals of 3 min between successive time points. A representative sequence of images is shown in Figure 5.

We simulate the in silico wounding of the vertex model by degrading the stiffness of the ablated cells down to 1% of their initial value, and removing the volume penalization term in the total energy for these cells, that is, reducing λ_Vol to 0 on those ablated cells. This softening, together with the subsequent removal of these degraded cells and the progressive intercalation of the cells at the wound edge allows us to simulate progressive wound closure.

We explicitly implement the evolution of an actomyosin purse-string at the wound edge by applying a decreasing trend to the apical vertex segments, while applying a constant value at the lateral sides. In the model, the contractilities are explicitly given by,

$\begin{array}{l}{\widehat{{\rm Y}}}_A^c=\{\begin{array}{ll}{{\rm Y}}_A^c\left(1-\frac{t-t_w}{400}\right),& t>t_w\\ 0,& \text{otherwise}\end{array}\\ {\widehat{{\rm Y}}}_L^c=\{\begin{array}{ll}{{\rm Y}}_L^c,& t>t_w\\ 0,& \text{otherwise}\end{array}\end{array}$

where t_w ≈ 6 min is the time at which purse-string contractility is activated, with t = 0 the time for tissue ablation. Our simulations last in average around 150 min, so that the factor 1/400 in ${\widehat{{\rm Y}}}_A^c$ aims at reducing the purse-string contractility with a similar trend of the concentration of Myosin II measured and showed in Figure 4. Model parameters ${{\rm Y}}_A^c$ and ${{\rm Y}}_L^c$ will be calibrated in section 3.1.

FIGURE 4

Figure 4. Purse-string Myosin II intensity and wound area evolution during the first hour of wound closure, averaged over 5 wing disc wounds. Myosin II intensity in the purse-string gradually reduces over time, as the wound closes (Error bars = S.D.).

FIGURE 5

Figure 5. Drosophila wing disc wound evolution. Apical (a,b,c) and lateral views (a’,b’,c’) of a Drosophila wing disc expressing a GFP-tagged form of non-muscle Myosin II at 0 min (a,a’), 6 min (maximum recoil, b,b’) and 45 min (partial closure, c,c’) after wounding. The margin of killed cells (a) and the wound (b,c) are marked by a green dotted line. The apical surface of the tissue is marked by a magenta dotted line (a’,b’,c’). Scale bar = 5μm.

2.3. Experimental Measurements

Experimental quantifications were averaged across wounds in five separate wing discs. We quantified the evolution of wing disc apical wound area by manually tracing the periphery of the wound for the first 72 min after wounding. The wound periphery was particularly clear while wound healing progressed, due to the formation of an apical actomyosin purse-string (Figure 5). To quantify the evolution of wound apical indentation depth, we first generated orthogonal image views and fitted a line between the highest points of the apical surface either side of the wound edge. We then calculated the distance along the apicobasal axis between this line and the position of the wound periphery (most clearly marked by the actomyosin purse-string, Figure 5). The relative height was computed from the depth measurements by assuming that the monolayer had an average height of 35 μm.

We also quantified the number of T1 transitions by analysing the gradual reduction in the number of cells at the wound edge. The evolution of the latter will be compared with our simulations in the next section. Figure 6 demonstrates that T1 transitions also occur spatially, along the apicobasal axis in a single timeframe, as well as temporally in the plane of the epithelium. This justifies our inclusion of intermediate vertices (Figure 2).

FIGURE 6

Figure 6. Wing disc cell outlines (cells a–h, marked by CAAX-GFP) at two different apicobasal positions, demonstrating that neighbor exchanges occur along the apicobasal axis at a single time point. Persistent cell-cell junctions are in blue, lost junctions in red and gained junctions in yellow.

3. Results

3.1. Model Calibration

We use the recoil process for calibrating tissue contractility ε^c, vertex stiffness k_V0, and the remodeling rate γ, which measures the viscous response. In order to avoid stiffness redundancy between vertex and nodal networks, we fix k_N = 0.5. Figure 7 shows the sensitivity of the recoil to these material parameters.

FIGURE 7

Figure 7. Apical wounded area as a function of time for the first 6 min following ablation as simulated using a 65 cell computational tissue. Different values of (A) background contractility ε^c, (B) vertex stiffness k_V, and (C) and remodeling rates γ.

Tissue contractility ε^c increases the asymptotic line tension in cells, and as such increases the final apical and basal areas after recoil. Larger values of recoil for MyoII activation and higher tension have been reported experimentally (Tetley et al., 2019), supporting this effect. Material stiffness has an inverse trend, reducing recoil for larger values of k_V0, which increases the relative energy cost of length changes. In our formulation, mechanical equilibrium is reached through energy minimization, which is proportional to material stiffness and line strains. Consequently, when stiffness increases, line stretching is in general reduced.

The remodeling rate mimics the viscous and fluid response of the cells. For higher values of γ, the characteristic time and the viscosity of the fluid is reduced, in agreement with the model (Muñoz and Albo, 2013). Consequently, the area after recoil increases for a fixed time. The final values of the tissue are given in Table 1. Other values of nodal and vertex stiffness, k_N and k_V have been manually fitted so that the final wound area remains stable at the experimental values, and that no element is under compression. Although there is some redundancy on their values (multiple combinations giving similar wound area), we have chosen the values k_N = 0.3 and k_V = 1.0.

TABLE 1

Table 1. Values of parameters fitted in recoil phase.

3.2. Wound Healing Simulations

We tested squared patches with different size, from 80 to 205 cells, and generated for each case cell positions similar to those measured experimentally. We fixed the positions of the cell-centers at the patch boundary. For the tested sizes of ablation (from 5 to 11 cells) and patch dimensions, the assumption of zero displacements on those external cells agreed also with the observed deformations. We measured the experimental displacements of the boundary cell centers, and the mean of their norm was in all cases below 1% of the side of the patch.

In order to test the effects of the wound edge contractility after ablation, we measured the time evolution of the relative projected area and the relative height at the wound edge in the in silico model and in vivo. Figure 8 shows snapshots of the full simulated tissue. Two videos showing apical and basal view of the simulation can be found in the Supplemental Material (Supplementary Videos 1, 2). Figure 9 shows cross-sections through the wounded region, where unequal closure at the apical and basal sides can be observed. Figure 10 shows the standard deviation of the experimental and simulated wound area evolutions, when using the same cell center positions and similar areas for each one of the patches tested. The deviations from the mean trend are similar, but the mean values of the numerical simulations close slightly earlier than the in vivo wounds. More sophisticated contractility profiles were needed in order to delay the simulated closure.

FIGURE 8

Figure 8. Computational epithelial tissue with 205 cells and 8 ablated cells, shown in side view for applied apical contractility of ${{\rm Y}}_A^c=2.3$. (A) Apical view of epithelial computational tissue shows the maximum recoil (t = 10 min). (B) Wound edge cell before intercalation t = 30 min. (C) Wound edge cell after intercalation t = 32 min. (D) Final time step of wound closure process t = 52 min. The color map illustrates the relative volume change. Videos of the simulation showing the apical and basal view can be found in the Supplemental Material.

FIGURE 9

Figure 9. Computational epithelial tissue with 205 cells and 8 ablated cells, shown in cross-section view for applied apical contractility of ${{\rm Y}}_A^c=2.3$. (A) Apical of epithelial computational tissue shows the maximum recoil t = 10 min. (B) Wound edge cell before intercalation t = 30 min. (C) Wound edge cell after intercalation t = 52 min. (D) Final time step of wound closure process. The color map illustrates the relative volume change.

FIGURE 10

Figure 10. Mean values and standard deviation of the experimental and numerical area evolution.

The parameter study of apical contraction ${{\rm Y}}_A^c$ in Figure 11 shows that below a threshold value of ≈ 2.2, contractility is insufficient to close the wound, for the tested material parameters. However, the changes on the purse-string tension have a minimal effect on the relative height evolution (see Figure 11B). We remark that the material parameters also have an effect on the closure time and profile: higher values of stiffness and tissue contractility ε^c delay or may prevent closure, while higher viscosity (lower value of γ) may also delay the closure process. These material properties have been calibrated in order to match the recoil, but also the rate of closure.

FIGURE 11

Figure 11. Area (A) and height (B) evolution as a function of apical purse-string contractility ${{\rm Y}}_A^c$.

Recently, it has been shown that changes in tissue height occur as a result of increased apicobasal contractility (Monier et al., 2015; Sui et al., 2018). We therefore also tested the effect of increased lateral contractility on wound closure. For increasing values of ${{\rm Y}}_L^c$, the height diminishes, as expected, and the closure of the area is in turn also accelerated (see the plots in Figure 12). Despite the fact that lateral edges are inclined due to the higher area reduction in the apical side, and oppose closure, their global effect is to contribute to closure due to the volume constraint.

FIGURE 12

Figure 12. Area (A) and height (B) evolution as a function of lateral purse-string contractility ${{\rm Y}}_L^c$.

We also analyse the evolution of the transitions, or equivalently, the number of cells at the wound edge. Figure 13A shows the evolution of the number of cells at the wound edge. While the experimental evolution is progressive, our simulations exhibit a more sudden concentration of the transitions. This may be due to the geometrical control of the transitions in the vertex discretization. We are currently investigating more accurate cell descriptions in order to obtain less drastic T1 transitions.

FIGURE 13

Figure 13. (A) Evolution of the mean number of junctions at apical wound perimeter in vivo and in the simulations. (B) Evolution of volume in the wounded region as a function penalization factor λ_Vol. The cases with λ₁0 and λ₂5 did not converge.

In order to measure the effect of the volume constraint, we checked the evolution of the apical area for different values of λ_Vol. Figure 13B shows that in fact, to strict or too relaxed volume preservation may impede wound closure. In our simulations, when λ_Vol = 25 or λ_Vol = 10, the wound does not close. In the first case this is due to the need to increase the cell size to some extent in order to recover the same sized patch with fewer cells (we do not simulate cell proliferation). In the second case, when λ_Vol is too low, the recoil is too large and cell adapt to the purse-string by changing size instead of closing the wound through intercalation. In our simulations we used λ_Vol = 20, which garantees that the mean deviation of each cells remains below 10%.

4. Conclusions and Discussion

Recent experimental analyses have shown that wound closure is not only driven by tissue tension and contractility, but also by the rate of intercalation, so-called tissue fluidity (Tetley et al., 2019). The present model aims at extending this analysis to three dimensions, by including lateral contractility, and allowing the simulation of different intercalation at the apical and basal side.

We used the recoil process to calibrate the material parameters that characterize tissue viscoelastic properties. The rate of wound expansion just after ablation is a useful measure of tissue viscosity, which we simulated through a variable rest-length L(t) that adapts according to a remodeling rate γ. This parameter, and tissue contractility ε^c have been experimentally measured in vitro (Wyatt et al., 2020) for suspended monolayers. We showed that ε^c takes the value ≈ 0.3, which is close to our fitted value from the in vivo measurements of the area evolution. We here also showed that the recoil can be employed to evaluate this material property and that it indeed determined the final expansion.

Based on experimental observations of the evolution of myosin concentration, we applied non constant trends on the apical and lateral surfaces of the wound. We encountered a minimum value of apical tension at the wound edge, below which no closure takes place. When ${{\rm Y}}_A^c$ is lower than approximately 2.3, the tissue in unable to surmount line tension between cells, given by ε^c in the model. It appears thus that the ratio between purse-string tension (given in the model by ${{\rm Y}}_A^c$) and cell line tension (in the model represented by ε^c) modulates the speed of closure, and that for too low values, closure may not succeed. Experimentally, too high values of line tension have been also shown to slow down or even prevent closure (Tetley et al., 2019).

We additionally observed that although lateral tension contributes in general to wound closure, depending on the tissue thickness, the net contribution of lateral purse-string contractility, regulated by ${{\rm Y}}_L^c$, the duration of the healing process is shortened. Further inspection of the whole shape and cross-section of the wound along the apicobasal axis and accurate measurements of the height profile are necessary to corroborate this fact. This is a challenging task, given the high aspect ratio (thickness/diameter) of the cells and their high light scattering. We also note that purse-string apical tension has very minor effects on the height evolution, but lateral contractility does modify the values of tissue thickness and importantly, also area evolution (see Figure 12).

Our numerical results indicated that apical purse-string tension, when applied together with volume preservation induced a reduction of the height, due to the expansion of the tissue. We point out though that in our model, lateral contractility is applied on the whole height of the tissue. This may not be so in the real tissue, where lateral myosin may not be homogeneous along the thickness. Further discretization of the monolayer along the apicobasal axis seems necessary to simulate the specific localization of lateral myosin on the apical side. We expect that more accurate experimental measurements and model enhancements will allow us to quantify cell mechanical contribution and regulation more closely during wound closure, which as shown, are more complex than a 2D analysis may reveal.

The evolution of relative height in vivo and in vitro plotted in Figures 11B, 12B also revealed that the initial height of the monolayer was not recovered, even when the wound was fully closed. In our model, this fact can be explained by the reduced number of cells in the patch after ablation. The reduction in the total volume reduction is compensated by a height reduction, since the patch area is fixed. In the experiments, whether the patch recovers the initial height after a sufficiently long period is still under study. This analysis, and the evolution of the material properties after successive re-wounding is left for further investigations.

Data Availability Statement

The data generated and analyzed for this study can be provided upon request to the corresponding authors.

Author Contributions

YM conceived the project and designed the experiments. JM conceived the in silico model. FI, MD, and JM developed the numerical simulations. RT carried out experiments and performed experimental quantification. All authors participated in the discussions of the experiments and the modeling. JM, FI, and YM drafted the manuscript, which was completed and revised by all the authors.

Funding

JM and MD have been financially supported by the Spanish Ministry of Science, Innovation and Universities (MICINN) with grant DPI2016-74929-R and by the local government Generalitat de Catalunya with grant 2017 SGR 1278. RT was funded by a Medical Research Council Skills Development Fellowship (MR/N014529/1). YM was funded by a Medical Research Council Fellowship MR/L009056/1, a Lister Institute Research Prize, and EMBO Young Investigator Programme. FI was funded by a Marie Curie ITN PolarNet. This work was also supported by MRC funding to the MRC LMCB University Unit at UCL, award code MC_U12266B.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2020.00405/full#supplementary-material

Supplementary Video 1. Apical view of the wound healing simulation using fitted parameters.

Supplementary Video 2. Basal view of the wound healing simulation using fitted parameters.

References

Alt, S., Ganguly, P., and Salbreux, G. (2017). Vertex models: from cell mechanics to tissue morphogenesis. Philos. Trans. R. Soc. Lond. B 372:20150520. doi: 10.1098/rstb.2015.0520

PubMed Abstract | CrossRef Full Text | Google Scholar

Antunes, M., Pereira, T., Cordeiro, J. V., Almeida, L., and Jacinto, A. (2013). Coordinated waves of actomyosin flow and apical cell constriction immediately after wounding. J. Cell Biol. 202, 365–379. doi: 10.1083/jcb.201211039

PubMed Abstract | CrossRef Full Text | Google Scholar

Barton, D., Henkes, S., Weijer, C., and Sknepnek, R. (2017). Active vertex model for cell-resolution description of epithelial tissue mechanics. PLoS Comput. Biol. 13:e1005569. doi: 10.1371/journal.pcbi.1005569

PubMed Abstract | CrossRef Full Text | Google Scholar

Bi, D., Lopez, J., Schwarz, J., and Manning, M. (2015). A density-independent rigidity transition in biological tissues. Nat. Phys. 11, 1074–1079. doi: 10.1038/nphys3471

CrossRef Full Text | Google Scholar

Bi, D., Yang, X., Marchetti, M. C., and Manning, M. L. (2016). Motility-driven glass and jamming transitions in biological tissues. Phys. Rev. X 6:021011. doi: 10.1103/PhysRevX.6.021011

PubMed Abstract | CrossRef Full Text | Google Scholar

Bielmeier, C., Alt, S., Weichselberger, V., Fortezza, M., Harz, H., Jülicher, F., et al. (2016). Interface contractility between differently fated cells drives cell elimination and cyst formation. Curr. Biol. 25, 563–574. doi: 10.1016/j.cub.2015.12.063

CrossRef Full Text | Google Scholar

Brugués, A., Anon, E., Conte, V., Veldhuis, J., Gupta, M., Collombelli, J., et al. (2014). Forces driving epithelial wound healing. Nat. Phys. 10, 683–690. doi: 10.1038/nphys3040

PubMed Abstract | CrossRef Full Text | Google Scholar

Deforet, M., Hakim, V., Yevick, H., Duclos, G., and Silberzan, P. (2014). Emergence of collective modes and tri-dimensional structures from epithelial confinement. Nat. Commun. 5:3747. doi: 10.1038/ncomms4747

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisenhoffer, G., Loftus, P., Yoshigi, M., Otsuna, H., Chien, C. B., Morcos, P. A., et al. (2012). Crowding induces live cell extrusion to maintain homeostatic cell numbers in epithelia. Nature 484, 546–549. doi: 10.1038/nature10999

PubMed Abstract | CrossRef Full Text | Google Scholar

Gelbart, M., He, B., Martin, A., Thiberge, S., Wieschaus, E., and Kaschube, M. (2012). Volume conservation principle involved in cell lengthening and nucleus movement during tissue morphogenesis. Proc. Natl. Acad. Sci. U.S.A. 109, 19298–19303. doi: 10.1073/pnas.1205258109

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez-Gálvez, P., Vicente-Munuera, P., Tagua, A., Forja, C., Castro, A., Letrán, M., et al. (2018). Scutoids are a geometrical solution to three-dimensional packing of epithelia. Nat. Commun. 9:2960. doi: 10.1038/s41467-018-05376-1

CrossRef Full Text | Google Scholar

Hashimoto, H., Robin, F., Sherrard, K., and Munro, E. (2015). Sequential contraction and exchange of apical junctions drives zippering and neural tube closure in a simple chordate. Dev. Cell 32, 241–255. doi: 10.1016/j.devcel.2014.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Khalilgharibi, N., Fouchard, J., Asadipour, N., Barrientos, R., Duda, M., Bonfanti, A., et al. (2019). Stress relaxation in epithelial monolayers is controlled by the actomyosin cortex. Nat. Phys. 15, 839–847. doi: 10.1038/s41567-019-0516-6

CrossRef Full Text | Google Scholar

Misra, M., Audoly, B., Kevrekidis, I., and Shvartsman, S. (2016). Shape transformations of epithelial shells. Biophys. J. 110, 1670–1678. doi: 10.1016/j.bpj.2016.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Monier, B., Gettings, M., Gay, G., Mangeat, T., Schott, S., Guarner, A., et al. (2015). Apico-basal forces exerted by apoptotic cells drive epithelium folding. Nature 7538, 245–248. doi: 10.1038/nature14152

CrossRef Full Text | Google Scholar

Mosaffa, P., Rodríguez-Ferran, A., and Muñoz, J. (2018). Hybrid cell-centred/vertex model for multicellular systems with equilibrium-preserving remodelling. Int. J. Num. Methods Biomed. Eng. 34, 1–24. doi: 10.1002/cnm.2928

PubMed Abstract | CrossRef Full Text | Google Scholar

Munjal, A., Philippe, J., Munro, E., and Lecuit, T. (2015). A self-organized biomechanical network drives shape changes during tissue morphogenesis. Nature 524, 351–355. doi: 10.1038/nature14603

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz, J., and Albo, S. (2013). Physiology-based model of cell viscoelasticity. Phys. Rev. E 88:012708. doi: 10.1103/PhysRevE.88.012708

PubMed Abstract | CrossRef Full Text | Google Scholar

Okuda, S., Inoue, Y., Eiraku, M., Adachi, T., and Sasai, Y. (2015). Vertex dynamics simulations of viscosity-dependent deformation during tissue morphogenesis. Biomech. Model. Mechanobiol. 14, 413–425. doi: 10.1007/s10237-014-0613-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Okuda, S., Inoue, Y., Eiraku, M., Sasai, Y., and Adachi, T. (2013). Modeling cell proliferation for simulating three-dimensional tissue morphogenesis based on a reversible network reconnection framework. Biomech. Model. Mechanobiol. 12, 987–996. doi: 10.1007/s10237-012-0458-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Sánchez-Corrales, E., Blanchard, G., and Röper, K. (2018). Radially patterned cell behaviours during tube budding from an epithelium. eLife 7:e35717. doi: 10.7554/eLife.35717

PubMed Abstract | CrossRef Full Text | Google Scholar

Staddon, M., Bi, D., Tabatabai, A. P., Ajeti, V., Murrell, M., and Banerjee, S. (2018). Cooperation of dual modes of cell motility promotes epithelial stress relaxation to accelerate wound healing. PLoS Comput Biol. 14:e1006502. doi: 10.1371/journal.pcbi.1006502

PubMed Abstract | CrossRef Full Text | Google Scholar

Sui, L., Alt, S., Weigert, M., Dye, N., Eaton, S., Jug, F., et al. (2018). Differential lateral and basal tension drive folding of drosophila wing discs through two distinct mechanisms. Nat. Commun. 4320, 1–13. doi: 10.1038/s41467-018-06497-3

CrossRef Full Text | Google Scholar

Sunyer, R., Conte, V., Escribano, J., Elosegui-Artola, A., Labernadie, A., Valon, L., et al. (2016). Collective cell durotaxis emerges from long-range intercellular force transmission. Science 353, 1157–1161. doi: 10.1126/science.aaf7119

PubMed Abstract | CrossRef Full Text | Google Scholar

Tetley, R., Staddon, M., Banerjee, S., and Mao, Y. (2019). Tissue fluidity promotes epithelial wound healing. Nat. Phys. 15, 1195–1203. doi: 10.1038/s41567-019-0618-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Tozluoğlu, M., Duda, M., Kirkland, N., Barrientos, R., Burden, J., Muñoz, J., et al. (2019). Planar differential growth rates initiate precise fold positions in complex epithelia. Dev. Cell 51, 299–312. doi: 10.1016/j.devcel.2019.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Wyatt, T. J., Fouchard, J., Lisica, A., Khalilgharibi, N., Baum, B., Recho, P., et al. (2020). Actomyosin controls planarity and folding of epithelia in response to compression. Nat. Mat. 19, 109–117. doi: 10.1101/422196

PubMed Abstract | CrossRef Full Text | Google Scholar

Zulueta-Coarasa, T., Tamada, M., Lee, E., and Fernandez-Gonzalez, R. (2014). Automated multidimensional image analysis reveals a role for Abl in embryonic wound repair. Development 141, 2901–2911. doi: 10.1242/dev.106898

PubMed Abstract | CrossRef Full Text | Google Scholar