Junctional and cytoplasmatic contributions in wound healing

Wound healing is characterised by the re-epitheliation of a tissue through the activation of contractile forces concentrated mainly at the wound edge. While the formation of an actin purse string has been identified as one of the main mechanisms, far less is known about the effects of the viscoelastic properties of the surrounding cells, and the different contribution of the junctional and cytoplasmic contractilities. In this paper we simulate the wound healing process, resorting to a hybrid vertex model that includes cell boundary and cytoplasmatic contractilities explicitly, together with a differentiated viscoelastic rheology based on an adaptive rest-length. From experimental measurements of the recoil and closure phases of wounds in the Drosophila wing disc epithelium, we fit tissue viscoelastic properties. We then analyse in terms of closure rate and energy requirements the contributions of junctional and cytoplasmatic contractilities. Our results suggest that reduction of junctional stiffness rather than cytoplasmatic stiffness has a more pronounced effect on shortening closure times, and that intercalation rate has a minor effect on the stored energy, but contributes significantly to shortening the healing process, mostly in the later stages. Author summary We simulate the wound healing process of epithelia in the absence of substrate. By analysing the recoil process we are able to fit the viscoelastic properties of the monolayer, and study the influence of contractility at junctions and at the interior polymer network. We numerically simulate the whole wound opening and closure process, and inspect which mechanism has a more pronounced effect in terms of energy barrier and wound closure rate. We conclude that while junctional stiffness seems to be more effective than bulk stiffness at speeding up the closure, the increase of intercalation process is the mechanism with the lowest energy cost.

Wound healing is a fundamental process for maintaining integrity and functionality 2 in epithelia, tissues, and organisms [1]. The most accepted mechanisms for wound 3 closure are Rac-dependent crawling using cell protrusions [2], and the formation of a 4 contractile supra-cellular actomyosin cable at the wound edge (purse-string 5 mechanism) [3]. However, in embryonic and larval tissues, either the elastic substrate is 6 not present or no crawling has been observed. In these cases, the first mechanism is 7 absent and closure relies mostly on the actomyosin cable and cell-cell interactions [4], 8 which is the focus of our study in this article. 9 Extensive analyses have reported the key mechanical factors affecting wound closure. 10 The effects of myosin heterogeneity [5,6], purse-string assembly rate [6], or tissue 11 fluidity [4] have been studied. However, far less is known about the contribution of 12 junctional and cell medial contractility and stiffness during closure, and the energy 13 barrier that they may impose during the intercalation process. We here simulate 14 wounding of the Drosophila imaginal disc epithelium [4], and numerically analyse the 15 influence of junctional versus medial cortex stiffness and contractility in a vertex model. 16 Their differential roles during tissue remodelling have been already reported [7], and 17 thus we expect that their different properties will influence cell polarisation and 18 intercalation that takes place during wound healing. 19 The modelling of wound healing on substrates using vertex models has been 20 presented for instance in [6,[8][9][10][11], while homogenised continuum models can be found 21 in [12][13][14], which allow the inclusion of myosin and calcium dynamics [13,15]. Due to 22 the importance of junctional mechanics, and the absence of a force contribution between 23 cells and their underlying substrate in our reference ex vivo system, we will employ a 24 customised vertex model with no crawling. However, both medial and junctional 25 cortices are differentiated in conjunction with a viscoelastic rheological model. 26 The tailored vertex model and the comparison with the experimental rate of tissue 27 recoil will enable us to calibrate the material viscoelastic properties. The employed 28 rheological model resorts to a Maxwell-like viscous model where the rest-length is able 29 to dynamically adapt [16]. Similar rest-length based models have been recently adopted 30 when modelling embryonic tissues [17,18] or stress relaxation in suspended 31 monolayers [19]. We will apply these ideas to the wound healing process. 32 We will assume uniform constant viscoelastic material properties, as it have been 33 measured experimentally [5], and uniform ring contractility at the wound edge. 34 Although is has been experimentally observed that the supra-cellular actin cable is not 35 uniformly continuous [5,10], we will coarsen its net effect at the wound front.

36
Our numerical model incorporates some differences with respect to our previous 37 hybrid model [9], mainly on the treatment of the wound edge and wound interior. We 38 first highlight these differences, and then discuss our results and compare them to 39 relevant experimental outputs. We place special emphasis on the role of the 40 intercalation and remodelling process, as recently studied in [4]. We additionally inspect 41 the role of the intercalation in terms of contributions from the cell junctions and medial 42 forces across cell boundaries, and show that they have unequal effects on the evolution 43 of the stored energy and the closure time. The proposed computational approach is based on a discrete description of the tissue 47 forming a flat monolayer. The cell boundaries are described through a set of connected 48 vertices y I forming the vertex network, as is customary in vertex models [20]. In 49 addition, we also include the cell centres network, defined by a set of nodes x i , 50 i ∈ {1, 2, . . . , N nodes }, with N nodes the number of cells in the tissue, and a triangulation 51 that defines the connectivity between cells, i.e. cells i and j are connected if a bar 52 element ij is present in the triangulation.

53
In order to obtain the connectivities between cells, we use a Graded Delaunay triangulation (GDT). It corresponds to a modification of the Delaunay algorithm, where the aspect ratio of triangles is not necessarily optimal. Instead, connectivity between two adjacent triangles is swapped only if the aspect ratio of one of the triangles is larger than a parameter δ. More specifically, given two triangles (1) AspectRatio(T I ) is measured as the ratio between the radius of the circumcircle and 54 the incircle of triangle T I . Condition in (1) is checked for I ∈ {1, 2, . . . , N tri }, with N tri 55 the number of triangles in the tissue. We point out that as a consequence of condition 56 (1), the model will have a higher number of elongated cells as we increase parameter δ, 57 and that for δ = 0 standard Delaunay triangulation is used. For this reason, we call this 58 parameter intercalation restriction.

59
The topology and geometry of the cells centres can be defined by a matrix X = [x 1 , x 2 , . . . , x N nodes ] and a matrix T = [T 1 ; T 2 ; . . . ; T Ntri ] defining the connectivity of the cells. The positions of cell centres and vertices in the two networks are not independent. Vertices are in general located at the barycentres of each triangle. In other words, given a triangle T I , the coordinates of the associated vertex y I are given by, with N i three interpolation weights, which we set to N I i = 1 3 , for i = 1, 2, 3 and all 60 triangles I = 1, . . . , N tri . While the cell centre network eases the description of the 61 connectivity changes, the vertex network allows us to include the cortex contractility, as 62 will be explained later. We note that similar coupled interpolations have been used on 63 other vertex models [21,22]. vertices are defined, preserving the initial area of cells, while mechanical equilibrium is 71 described in the next section.

72
As a result, the whole tissue geometry and connectivity is given by matrices X, T connectivity T is recomputed at each new time step t n+1 from the resulting coordinates 76 resorting to the GDT. Figure 1 illustrates the process of obtaining Since the external boundary of the cells at the tissue edge is not defined directly through tessellation of the triangular nodal network, the cell connectivity at the tissue edge cannot be computed by resorting to the GDT. Instead, we employ another criterion which is based on the ratio between the length of the cell side at the tissue edge, l, and the cell perimeter, P . An intercalation is applied when this ratio is smaller than a certain threshold. The restriction on intercalation on the external cells then reads, where the numerical values in 3 have been chosen so that the intercalations at the tissue 80 edge resemble those in the bulk of the tissue, while being controlled by the same 81 parameter δ. Similar to other vertex models [6,11,20], our equilibrium equations are constructed by minimising a total potential function, which in our case is given by, The terms W D (x) and W V (y(x), y b ) account respectively for the contribution of the bar elements in the cell centre and vertex networks. We also add the term W A (y(x), y b ) which penalises area changes. The specific form of these contributions will be detailed below. As yet, the equations of mechanical equilibrium can be derived by minimising the total elastic energy W with respect to the system principal kinematic variables (degrees of freedom): nodal positions in X and the position of exterior vertices in Y b . This minimisation yields the following system of non-linear equations, The derivatives given above can be obtained through standard mathematical Newton-Raphson process, which requires full linearisation of the equations. The explicit 87 expression of the equations and the iterative scheme can be found in [9].

88
Nodal and vertex mechanics

89
Triangulation T of the tissue domain results in N D pairs ij of connected nodes, defining a bar element each, and N V pairs IJ of vertices. Each nodal bar element represents the forces between two cells, while vertex bar elements represent forces at cell junctions. Forces are deduced from the nodal and vertex energy function W D and W V , which read, with the elemental functions given by Material parameters k D and k V are respectively the intra-cellular and cell boundary 90  As the number and the size of the cells is considered approximately constant within the tissue, we consider a cell volume invariance in the model, which is applied here by imposing a two-dimensional area constraint through the energy term, where λ A is a penalisation coefficient and A i 0 and A i are the initial and the current 96 areas of cell i, respectively. The value of λ A is chosen such that relative area variations 97 do not exceed 10%, as as experimentally reported in 3D measurements [5]. 98 Tissue Ablation

99
Wounding in the Drosophila wing disc is carried out by laser ablation of a number of 100 cells within a region located at the interior of the tissue [4]. In the proposed 101 computational approach, we simulate this by removing a number of cell centre nodes 102 inside a circle with radius r and centred at the interior of the nodal network, which 103 represents the area ablated in experiments. This will leave an empty area in the interior 104 of the nodal network, as well as a new tissue edge (wound boundary), which needs to be 105 constructed by new non-coupled boundary vertices y W b .
106 Figure 3 depicts the wounding process in the vertex model. Due to the tissue 107 elasticity, a recoil of the tissue and associated stretching of the boundary at the wound 108 edge takes place.

Rheological model 110
The strain measures in (7) are written in terms of the difference between the measurable current length l and the rest-length L, which is used as an internal variable, and is not necessarily constant. In order to mimic the viscoelastic response of cells and cell junctions [23][24][25], we resort to a dynamic law of the rest-length evolutioṅ where γ is the remodelling rate, and ε is the elastic strain used in (7). It has been 111 previously shown that such a rheological model reproduces Maxwell viscoelastic 112 behaviour [16], and that it can be employed to simulate tissue fluidisation [19,26] or cell 113 cortex response in embryogenesis [17,27].

114
Since cells exhibit inherent contractility exerted at the cortex cross-linked structure [28], we modify the evolution law in (9) aṡ February 17, 2020 6/17 with ε c a contractility parameter. Note that according to this law, the rest-length will 115 remain constant whenever the elastic strain ε equals the value of the contractility ε c , 116 which in each simulation takes a constant value.

117
The differential equation in (10) is discretised resorting to a Backward Euler scheme. 118 The unknown rest-lengths L n+1 are thus written in terms of the current lengths l n+1 119 and inserted in the non-linear equations in (5), which are solved at each time-step t n+1 , 120 allowing the new nodal positions x i n+1 and boundary vertices y m b to be found.

121
Due to cell reorganisation, new bar elements may be created at each time-step. This 122 requires the redefinition of the rest-lengths for some of the elements. We resort to an 123 Equilibrium Preserving Mapping, which computes new elemental rest-lengths that 124 minimise equilibrium errors at vertices and nodes. Further details of this mapping may 125 be found in s [9]. The hybrid vertex model has been implemented in Matlab (R2018, 126 The Mathworks, Inc.) and is available for its use (see online Supplemental Information). 127 We use different contractilities for the nodal and vertex networks, denoted Our numerical model depends on eight material and numerical parameters: the 138 material stiffnesses k D and k V , the remodelling rates γ D and γ V , the GDT intercalation 139 restriction, δ, and the contractilities ε c , ε c and ε c w .

140
We use the time evolution of the wounded area A(t) in order to fit these parameters. 141 Indeed, the recoil process allows us to set some reference values of material constants 142 (k, γ, ε) D and (k, γ, ε) V , while the closure rate is mostly determined by the values of ε c w 143 and δ, in conjunction also with (γ D , γ V ). In Appendix S2 we explain how each factor 144 affects the resulting curve A(t) and how the fitting is accomplished. Table 1 gives the 145 obtained reference values of the parameters, and in Figure 5 we plot the evolution of 146 A(t) for the experimental and numerical results using the reference values in Table 1.  It is worth pointing out that the presence of intercalation events at the wound edge 158 gives rise to a disruption of the cell orientation towards the wound edge [29], allowing 159 cells shape to relax towards a more isotropic distribution, despite imposing a purse 160 string mechanism. Furthermore, although we impose a strain driven closure (rest-length 161 at wound edge is such that the elastic strains converge towards the target value ε c w ), the 162 total stress and strain are not homogeneous, as also experimentally measured [5]. 163 With the aim of gaining further insight into the mechanisms that drive and prevent 164 wound healing and that may increase closure rate, we also plot the time evolution of the 165 elastic energy W E in Figures 7 and 8  In both plots the faster the wound closes, the higher the final stored energy, 172 although energy differences are relatively small. Also, we point out that energy levels at 173 closure are always lower than energy before ablation, but higher than at maximum 174 recoil. The fact that far fewer junctions are present during ablation and intercalation, 175 may explain this energy drop with respect to initial state before wounding. However, the application of ε c w implies an energy supply, which may be biologically explained by 177 actin flow, or higher consumption of ATP by myosin motors at the wound edge.  Interestingly, the intercalation restriction has no effect on the energy evolution. Indeed, 184 all the curves in Figure 10 collapse almost on a single energy-area curve, suggesting that 185 the presence of intercalation, despite minor fluctuations, has no justification in terms of 186 energy terms, but may be able to even reduce the closure time by approximately %75. 187 The vertex model presented here has the ability to distinguish between junctional 188 and non-junctional cortex contributions. For this reason, we have analysed the relative 189 effect of nodal and vertex stiffness parameters, k D and k V . Figure (11) indicates that   although increasing values of both stiffnesses lead to longer closure times, the increase 191 of vertex stiffness has a more pronounced effect on the total time. When inspecting the 192 energy profiles in Figure 12, we also observe that too large values of k V may prevent full 193 closure due to the wound contraction being insufficient to surmount the required energy 194 barrier. Also, for sufficiently low values of k V , the energy level after closure 195 approximates the one after the recoil phase, although a slight energy difference remains 196 in all cases.  We also inspect the number of intercalations during the closure process, N remodel .
198 Figure 13 shows the total number in the (k D , k V ) and (δ, ε c w ) diagrams. We note that 199 remodelling events may be occurring in the bulk tissue rather than specifically at the 200 wound edge, and consequently may not necessarily prevent the formation of rosette-like 201 wounds. For large values of δ, most of the intercalations take place at the later stages, 202 while for lower values, they start in earlier stages a proceed at a more constant rate. It is 203 also worth noting that by increasing δ at higher values of ε c w , the closure is nevertheless 204 obtained in absence of intercalations. In both cases though, N remodel has larger The dependence of the closure process on the intercalation rate and tissue fluidity 212 has been recently revealed [4] by varying forces at the cell boundaries. We show here 213 that this dependence also holds true for non-junctional cortically generated forces, with 214 similar but minor effects. Importantly, it has been shown that pulling forces triggered 215 by medial Myosin II flow cannot be discounted [30] in morphogenesis. Our results 216 suggest that such non-junctional cortically generated forces may similarly have an 217 impact on epithelial wound closure.

219
We have presented a vertex computational model that is able to distinguish 220 viscoelastic properties at cell junctions and cytoplasmatic regions. The model 221 parameters have been fitted from experimental wound recoil and closure measurements. 222 After inspection of the energy evolution as a function of stiffnesses in the vertex and 223 nodal network, (k D , k V ), wound contractility ε c w , and cell intercalation restriction δ, we 224 have identified different mechanisms that can increase closure rate. From our results, an 225 increase in intercalation rate or a reduction in vertex stiffness is more effective than 226 reducing cytoplasmic stiffness or increasing purse cable contractility. The latter induces 227 an important increase in the tissue stored energy confirming the actomyosin purse string 228 is actually dispensable for closure, as also suggested in dorsal closure of Drosophila 229 embryo [31]. 230 Our analyses show that the stored energy after closure is slightly higher than just 231 after wounding. This energy gap is higher for increasing actomyosin cable 232 contractiltities, but in all cases lower than the energy of the tissue before wounding. It 233 is yet experimentally unclear though whether the tissue recovers its initial stress and 234 energetic state. Further experiments with rewounding should allow us to quantify 235 whether reepithelialisation rates are altered in recently healed tissues. 236 Our results confirm that increasing cellular rearrangements can ease wound 237 closure [4]. We further studied the prominent role of junctions instead of other cell The radius of the arc, defining the initial position of the vertices is defined from the 256 initial area of the cell. The actual position of the vertices is found through mechanical 257 balance, as described in section Mechanical Equilibrium.

259
In this section, we briefly explain the fitting of the material parameters in order to 260 best mimic the wound area evolution during the course of healing in experimental 261 tissues.

262
Before proceeding to fit material parameters in order to simulate experimental 263 results, we studied the role of each parameter separately.

264
Wound healing in the Drosophila wing disc has been reported to consist of two 265 phases. The first phase corresponds to a nearly immediate fast recoil of the wound edge, 266 while the second phase consists of the gradual shrinkage of the wound area ending in the 267 full closure of the wound. These two phases can be identified in Figure 5.

268
Since we use the active model to simulate the rheology in both the nodal and vertex 269 networks, representing mechanics of cell bulk and boundary matter respectively, we 270 define the parameters stiffness k, remodelling γ and contractility ε c using subscripts D 271 and V for either of these networks, respectively. Figure

278
Note that the intercalation δ is not considered during the first phase, since due to 279 the short time scale of the area evolution in this phase, no cell intercalation is observed. 280 ε c w was neither considered during the first phase, since the actomyosin concentration 281 responsible for contractility gain at the wound edge is assumed to occur instantly with a 282 delay respecting the wounding instant. ε c w is denoted in units of vertex background 283 contractility ε c V . Note also that k D and k V have opposite effects on the peak of the area 284 before closure starts.