Wound opening in a thin incompressible viscoelastic tissue
WWound opening in a thin incompressible viscoelastic tissue
G. M. Carvalho,
1, 2, 3, ∗ N. A. M. Ara´ujo,
1, 2, † and P. Patr´ıcio
1, 4, ‡ Centro de F´ısica Te´orica e Computacional, Faculdade de Ciˆencias,Universidade de Lisboa, 1749-016 Lisboa, Portugal Departamento de F´ısica, Faculdade de Ciˆencias,Universidade de Lisboa, 1749-016 Lisboa, Portugal Instituto Federal de Educa¸c˜ao, Ciˆencia e Tecnologia Catarinense,89283-064 S˜ao Bento do Sul, Santa Catarina, Brasil Instituto Superior de Engenharia de Lisboa, Instituto Polit´ecnico de Lisboa, 1959-007 Lisboa, Portugal
We develop a model to investigate analytically and numerically the mechanics of wound openingmade in a viscoelastic, isotropic, homogeneous, and incompressive thin tissue. This process occursjust immediately after the wound infliction. Before any active biological action has taken place, thetissue relaxes, and the wound opens mostly due to the initial homeostatic tension of the tissue, itselastic and viscous properties, and the existing friction between the tissue and its substrate. We findthat for a circular wound the regimes of deformation are defined by a single adimensional parameter λ , which characterizes the relative importance of viscosity over friction. I. INTRODUCTION
Cell and subcellular dynamics, in vivo and in vitro,in response to external stimuli or during the embryonicdevelopment of plants and animals, have been the subjectof intense research activity during the last decades [1,2]. Among the most popular topics is wound healing,which is a process of tissue regeneration for wound closure.This implies collective cellular migration and formationof a contractile cable that connects the cells along thewound edge [3–5]. Modeling the biophysical mechanismsassociated with wound healing is a non trivial challenge [6–9].Recent advances on experimental techniques allow toaccess the mechanical properties of the tissues and followtheir dynamics in real time. These have opened thepossibility of studying the process of wound infliction andhealing in a systematic and quantitative way [10–14].Along with the development of experimental techniques,there is a need to develop theoretical and numerical mod-els to shed light on the biochemical and biophysical pro-cesses involved in tissue regeneration. There are severalways to model the movement of cells in a tissue [15–20]. The existing models are classified as continuous,particle-based or hybrid models. Continuous models areusually considered to access large length- and time-scales.Particle-based models are appropriate when some levelof detail of the particle-particle interaction is of rele-vance [20, 21]. Hybrid models combine properties of bothtypes [22]. Usually the wound healing models are macro-scopic and continuous, being used to investigate the globalbehavior of cells in the tissue [23–25].Inspired by recent experimental results [26], we proposea theoretical model to investigate analytically and numer-ically the deformation of epithelial tissues of Drosophila ∗ [email protected] † [email protected] ‡ [email protected] larvae after wound infliction. We focus our study in thefirst moments of the deformation, the first tens of seconds,before any active biological process comes into play. Atthis stage, the wound opens under the influence of theinitial homeostatic tension of the tissue, its elastic andviscous properties, and the friction between the epithelialtissue and its surroundings. We use a 3D Kelvin-Voigtcontinuous model, which combines both the elastic andviscous properties of the tissue, and allows the tissue tobe initially stretched, resisting to an existing homeostatictension. The tissue dynamics is given by the Newton’slaws in the overdamped regime. By choosing appropri-ate length and time scales, we find different deformationregimes, which depend on a unique adimensional param-eter λ , that characterizes the relative importance of theviscosity over friction.To our knowledge, this adimensional parameter was firstintroduced in the context of cell mechanics in Ref. [27].In this work, the authors severed in vivo the adherensjunctions around a disc-shaped domain of Drosophilapupa dorsal thorax epithelium, comprising typically ahundred cells. They compared the observed deformationof the disk, as it shrunk and relaxed, with the resultsobtained using a 1D Kelvin-Voigt model to find that therelative importance of viscosity over friction increasedwith pupa’s age (see Fig. 5 of Ref. [27]). In the contextof a Kelvin-Voigt model, this parameter also appears inRef. [28], where different rheological models are reviewed.In the following section, we give the mathematicaldetails of our continuum Kelvin-Voigt model. In the limitof an incompressible thin tissue, we obtain a 2D generalequation of motion (Eq. 10). In section III, we solve thisequation to obtain the dynamics of the tissue after theinfliction of a circular wound. We draw some conclusionsin the last section. a r X i v : . [ c ond - m a t . s o f t ] F e b II. MODEL
We model the tissue as an isotropic and homogeneousmaterial, whose mechanical deformation follows the 3DKelvin-Voigt model. The total stress tensor may be writ-ten as [29]: σ = 2 G ε + λ e Tr( ε ) I + 2 η γ + λ v Tr( γ ) I , (1)where ε is the strain tensor, γ is the strain-rate tensor, I is the identity tensor, G and λ e are the Lam´e elastic con-stants, and η and λ v are the dynamic and bulk viscositycoefficients. The first and second (third and fourth) termsof Eq. (1) correspond to the elastic (viscous) part of thestress tensor, and describe the elastic (viscous) responseof each volume element to forces applied tangential andnormal to the different surfaces of the element, respec-tively. In the Kelvin-Voigt model, elastic and viscousstress terms add up.We use the equilibrium position of the mate-rial as the vector coordinate reference X . Its de-formed position x ( X , t ) defines the displacement vec-tor u ( X , t ) = x ( X , t ) − X . For small displacements, thestrain tensor and the rate of strain tensor retain only thelinear terms in u : ε = 12 (cid:0) ∇ u + ( ∇ u ) T (cid:1) , γ = 12 (cid:0) ∇ v + ( ∇ v ) T (cid:1) , (2)where the differential operator ∇ is defined with respectto the coordinate reference X , and v = ∂ u /∂t = ˙ u .The dynamics of the tissue is described by Newton’slaw of motion: ρ D v Dt = ρ g + ∇ · σ , (3)where ρ is the tissue density and g the acceleration ofgravity. The total time derivative of the velocity is definedas D v /Dt = ∂ v /∂t + v · ∇ v . If the deformation of thetissue is of the order of the size of individual cells, wemay neglect the inertial terms and the equation of motionbecomes: ∇ · σ = 0 . (4)Let us suppose now that the tissue is also incompressible.In this case, we have:Tr( ε ) = Tr( γ ) = 0 . (5)By introducing the pressure field p , a Lagrange multiplierwhich ensures this condition, the stress tensor may bewritten in a simplified form: σ = 2 G ε + 2 η γ − p I . (6)A thin cellular tissue may be represented as a 2D surfacein the plane x − y . In this limit, we assume that the normalforces applied to the lower and upper sides of the tissueare much smaller than the longitudinal forces in the bulk. Since the tissue is thin, the normal forces inside the tissueare also negligible and so, σ zz = 0 ⇔ p = 2 Gε zz + 2 ηγ zz , (7)everywhere inside the tissue [30]. Taking into account theincompressibility condition, we obtain: p = − G Tr (cid:0) ε t (cid:1) − η Tr (cid:0) γ t (cid:1) , (8)where ε t and γ t correspond to the 2D strain and strain-rate tensors, defined in the plane x − y of the tissue.Therefore, in the 2D approximation of an incompress-ible isotropic Kelvin-Voigt tissue, the in-plane 2D stresstensor is given by: σ t = 2 G (cid:0) ε t + Tr (cid:0) ε t (cid:1) I t (cid:1) + 2 η (cid:0) γ t + Tr (cid:0) γ t (cid:1) I t (cid:1) , (9)where I t is the identity tensor in 2D.The friction between the tissue and the substrate isa tangential contact force, exerted on the bottom. Dueto the small thickness of the tissue, this force is spreadthrough the interior and it may be described, for thesake of simplicity, as an in-plane force per unit volume f t = − ζ v t . The 2D equation of motion becomes: ∇ t · σ t − ζ v t = 0 . (10) III. RESULTSA. Tissue under a uniform stress
We consider first a circular tissue under a uniformdistribution of forces, as shown in Fig. 1b). At rest,the radius of the circular tissue is R ∞ (see Fig. 1a).In cylindrical coordinates, the coordinate reference is: X t = ( r, θ ), 0 ≤ r ≤ R ∞ and 0 ≤ θ < π . The boundarycondition at the periphery of the tissue is σ rr ( r = R ∞ , θ, t ) = σ. (11)The radial symmetry of the tissue and applied forcessuggests a deformation in the form: u t ( r, θ, t ) = u ( r, t ) e r . (12)The in-plane strain tensor becomes: ε rr = ∂u∂r , ε θθ = ur , ε rθ = 0 . (13)The in-plane components of the stress tensor are: σ rr = 2 G (cid:18) ∂u∂r + ur (cid:19) + 2 η (cid:18) ∂ ˙ u∂r + ˙ ur (cid:19) , (14) σ θθ = 2 G (cid:18) ∂u∂r + 2 ur (cid:19) + 2 η (cid:18) ∂ ˙ u∂r + 2 ˙ ur (cid:19) , (15) σ rθ = 0 . (16) FIG. 1. a) A thin circular tissue of radius R ∞ at rest. b) Tissue under a uniform radial distribution of forces, where σ is the force per unit area, applied at the lateral border of thetissue, leaving it stretched. c) The tissue under tension, justimmediately after a circular wound, of radius R , has beeninflicted. There are no forces applied to the internal lateralborder of the tissue. d) Final deformation of the tissue, afterrelaxation. The hole increases to reach a new larger radius R . The gradation of colors reflects the intensity of the radialtension of the tissue. By symmetry, the polar component of the equation ofmotion is immediately satisfied. On the other hand, theradial component of the equation of motion (10) is: ∂σ rr ∂r + 1 r ( σ rr − σ θθ ) − ζ ˙ u = 0 , (17)which, after some algebra, becomes simply4 (cid:0) GD u + ηD ˙ u (cid:1) − ζ ˙ u = 0 , (18)where D u = ∂∂r (cid:18) r ∂∂r ( ru ) (cid:19) . (19)At equilibrium, the deformation of the tissue must obeythe equation D u = 0, which has the general solution: u ( r ) = ar + br , (20)and the coefficients a and b may be calculated from theboundary conditions. Since u ( r = 0) = 0 at the centerof the tissue, we have b = 0. The stress is in this caseconstant at every point of the tissue σ rr = σ θθ = 6 Ga = σ .Thus, we obtain: u ( r ) = σr G . (21)
B. Circular wound
We now consider that a circular wound, of radius R ,is made at the center of the tissue under tension (seeFig. 1c)). Some of the tension is released, and the holeincreases to reach a larger radius R (see Fig. 1d)). Thedeformation maintains its radial symmetry and the equa-tion of motion is the same (Eq. (18)).In the general case, this equation is solved numerically.However, in the limit of no friction: (cid:0) GD u + ηD ˙ u (cid:1) = 0 , (22)we may find analytical solutions. In this limit, we have D u = Ae − t/τ , where τ = η/G is a relaxation time and A an integration constant. A = 0, since initially, justimmediately after the wound, we have D u = 0. So, wehave the general solution: u ( r, t ) = a ( t ) r + b ( t ) r . (23)The coordinate reference of the wound radius r = R is given by: R = R + u ( R , t = 0) = R (cid:16) σ G (cid:17) . (24)We consider that no forces are applied to the tissue atthe wound border, so σ rr ( r = R , t ) = 0, or:2 G (cid:18) a ( t ) − b ( t ) R (cid:19) + 2 η (cid:32) a ( t ) − ˙ b ( t ) R (cid:33) = 0 . (25)At the periphery of the tissue, σ rr ( r = R ∞ , t ) = σ :2 G (cid:18) a ( t ) − b ( t ) R ∞ (cid:19) + 2 η (cid:32) a ( t ) − ˙ b ( t ) R ∞ (cid:33) = σ. (26)If R (cid:28) R ∞ , the final set of equations for a and b maythen be written in the form: a + τ ˙ a = σ G (cid:18) R ∞ R ∞ − R (cid:19) ≈ σ G , (27) b + τ ˙ b = σR G (cid:18) R ∞ R ∞ − R (cid:19) ≈ σR G . (28)Using the initial conditions: a ( t = 0) = σ G , b ( t = 0) = 0 . (29)we obtain: a ( t ) ≈ σ G , (30) b ( t ) ≈ σR G (cid:16) − e − t/τ (cid:17) . (31)In particular, we may now determine the time evolutionof the wound opening, R ( t ) = R + u ( R , t ). In the limitof small displacements, to linear order in σ/G , we have: R ( t ) ≈ R (cid:104) σ G (cid:16) − e − t/τ (cid:17)(cid:105) . (32)In the presence of friction, the deformation of the tissueunder tension is different. Nevertheless, their equilibriuminitial and final states is the same, as the velocities of theseconfigurations are zero. To obtain the time dependence ofthe deformation, we must solve this equation numerically.If we choose L = R and T = τ as unit length and timescales, the adimensional equation of motion is then: D u + D ˙ u − ˙ uλ = 0 , (33)where λ = L η /R , and L η = 2 (cid:112) η/ζ is the viscous length.The parameter λ sets the possible regimes of deforma-tion. When λ → ∞ , we recover the limit where frictionis negligible, discussed above.Numerically, it is convenient to use T = τ /λ as thetime unit. This choice yields the following equation ofmotion: D u + λ D ˙ u − ˙ u = 0 . (34)If λ →
0, we may use an explicit method to integratethis equation. It will be stable if we choose a sufficientlysmall time step. However, as λ increases, the methodrapidly gets unstable. To solve this equation for everychoice of λ , we use the implicit method described in theappendix, with second order precision in space and time.The boundary conditions are expressed through the stresstensor. The stress tensor is adimensionalized by dividingit by the elastic modulus G . Using this latter choice forunit space and time scales, we have: σ rr = Du + λ D ˙ u, (35)where D = 2 (cid:18) ∂∂r + 1 r (cid:19) . (36)So, the explicit adimensionalized boundary conditions are σ rr ( r = 1) = 0 and σ rr ( r = R ∞ /R ) = σ .Figures 2, 3 and 4 show the results for the displacement u ( r, t ) and the velocity v ( r, t ) = ˙ u ( r, t ), for three differentvalues of λ = 0 . , ,
10 (corresponding respectively to thefriction, intermediate and viscous regimes), R ∞ = 10 and σ = 1.The panel a) of each figure shows the displacement as afunction of the radius, for different instants t . The dashedred and black lines represent respectively the displacementbefore wound infliction ( t = 0), and at the final relaxedopened state of the wound ( t = ∞ ). These representationsshow already the differences between the three differentdeformation regimes. However these differences stand outin the other plots, as we discuss in what follows.The panels b) show the velocity as a function of theradius, also for different instants t it is clear that thevelocity is initially (at t = 0) higher closer the wound.The velocity decreases rapidly as we move away from thewound. In fact, the characteristic length of this decay isof the order of the viscous length L η = 2 (cid:112) η/ζ = λR . In D i s p l ace m e n t Radius0.000.450.901.351.80 1 4 7 10 0.000.501.001.502.00 1 4 7 10b) V e l o c it y Radiust = 0.0 t = 0.1 t = 0.5 t = 2.0 t = 5.0 t = 50.00.000.501.001.502.00 1 4 7 100.000.250.500.751.00 0 30 60 90 120c) N o r m . D i s p l ace m e n t Timer = 1.0 r = 2.5 r = 5.0 r = 7.5 r = 10.00.000.250.500.751.00 0 30 60 90 120 0.000.250.500.751.00 0 30 60 90 120d) N o r m . V e l o c it y Time0.000.250.500.751.00 0 30 60 90 1200.00.51.0 0 4 8 12Time0.00.51.0 0 4 8 12
FIG. 2. Tissue deformation after wound infliction, for λ = 0 . R ∞ = 10 and σ = 1 (friction regime). a) Displacement u ( r, t ) vs radius r , for different instants t of the deformation.The dashed red and black lines represent respectively thedisplacement before wound infliction ( t = 0), and at the fi-nal relaxed opened state of the wound ( t = ∞ ); b) Velocity v ( r, t ) = ˙ u ( r, t ) vs radius r , for different instants t of the de-formation; c) Normalized displacement (( u − u i ) / ( u f − u i ),with u i the initial displacement and u f the final displacement) vs time, in different positions r of the tissue. d) Normalizedvelocity ( v/v max , with v max the maximum velocity attainedat a particular position) vs time, in different positions r of thetissue. D i s p l ace m e n t Radius0.000.450.901.351.80 1 4 7 10 0.000.110.220.330.44 1 4 7 10b) V e l o c it y Radiust = 0.0 t = 0.1 t = 0.5 t = 2.0 t = 5.0 t = 50.00.000.110.220.330.44 1 4 7 100.000.250.500.751.00 0 30 60 90 120c) N o r m . D i s p l ace m e n t Timer = 1.0 r = 2.5 r = 5.0 r = 7.5 r = 10.00.000.250.500.751.00 0 30 60 90 120 0.000.250.500.751.00 0 30 60 90 120d) N o r m . V e l o c it y Time0.000.250.500.751.00 0 30 60 90 1200.00.51.0 0 4 8 12Time0.00.51.0 0 4 8 12
FIG. 3. The same as in Fig. 2, but for λ = 1 (intermediateregime). the friction regime (Fig. 2 b)), the decay length is verysmall ( L η = 0 . L η = 10), of the same size asthe tissue itself. The inital velocity decay length is thena signature of the regime, and may be of experimentalinterest to access the relative importance of the viscosityover friction.Panels c) show the time dependence of the normal-ized displacement ( u − u i ) / ( u f − u i ), with u i the initial D i s p l ace m e n t Radius0.000.450.901.351.80 1 4 7 10 0.0000.0030.0050.0070.010 1 4 7 10b) V e l o c it y Radiust = 0.0 t = 20.0 t = 50.0 t = 100.0t = 200.0t = 400.00.0000.0030.0050.0070.010 1 4 7 100.000.250.500.751.00 0 125 250 375 500c) N o r m . D i s p l ace m e n t Timer = 1.0 r = 2.5 r = 5.0 r = 7.5 r = 10.00.000.250.500.751.00 0 125 250 375 500 0.000.250.500.751.00 0 125 250 375 500d) N o r m . V e l o c it y Time0.000.250.500.751.00 0 125 250 375 500
FIG. 4. The same as in Fig. 2, but with λ = 10 (viscousregime). displacement and u f the final displacement, for differentpositions r . With these figures, it is possible to understandthat in friction regimes (small λ ), the wound opening oc-curs faster near the wound ( r = 1), whereas away fromthe wound ( r > . λ ), the tissue globally relaxes with a characteristic timewhich depends solely on the ratio between its viscous andelastic properties.Panels d) show, for different positions r , the time de-pendence of the normalized velocity v/v max , with v max the maximum velocity attained at a particular position.From these figures, it is clear that in friction regimes,the positions away from the wound do not move immedi-ately. Instead, their speeds increase up to a certain value,and only afterwards start to decrease until they achievetheir relaxed states. The results suggest that there isa propagation wave affecting the maximum normalizedvelocity for each position. This effect is however hardlynoticed, because the velocities far from the wound arealready very small. In the viscous regimes, the velocityprofile is almost the same at each point of the tissue: thetissue feels almost instantaneously the wound inflictioneverywhere. IV. FINAL REMARKS
We developed a theoretical model to investigate ana-lytically and numerically the mechanical deformation ofan epithelial tissue after a circular wound has been in-flicted. The tissue was described as a continuous isotropic,homogeneous and incompressible thin material, obeyingthe 3D Kelvin-Voigt model. This model takes into ac-count the elastic and viscous properties of the tissue, andallows the tissue to be stretched at equilibrium, underthe effect of a homeostatic pressure. Friction between the tissue and its surroundings was also considered. Wedetermined the passive mechanical response of the tissue,after wound infliction, without considering any active bio-logical effects which will try to close the wound, and healthe tissue. This behavior is consistent with the passivephysical behaviour observed for the first tens of secondsfor the Drosophila larvae studied in Ref. [26].By choosing appropriate length and time scales, wefound different deformation regimes, depending on aunique adimensional parameter λ , which characterizesthe relative importance of the viscosity η over friction ζ . Although the final relaxed state is the same, for allvalues of the viscosity or friction, the dynamics of thedeformation presents distinct features. In friction regimes,for small λ , the deformation is initially concentrated atthe border of the wound, whereas in viscous regimes, forlarge λ , the deformation evolves globally. In fact, if onlyviscosity is present everywhere, and only depends frictionis negligeable, the normalized displacement time evolutionis exactly the same, everywhere, and only depends on arelaxation time simply defined by the ratio between theviscous and elastic properties of the tissue.The experimental characterization of these differentregimes may be accessed through the initial velocity spaceprofile, v ( r, t = 0). The initial velocity field typicallydecays from the wound boundary in a typical lengthgiven by L η = λR = 2 (cid:112) η/ζ , which relates viscosity andfriction. Appendix: Numerical integration
To integrate the equations numerically, we discretizespace and time: r i = 1 + ( i − r, ( i = 1 , ..., N ) , (A.1) t n = n ∆ t, ( n = 0 , , ... ) , (A.2)with ∆ r = ( R ∞ − / ( N − u ( r i , t n ) = u ni ,for simplicity. The discretized adimensional equation ofmotion (Eq. (34)) is approximated by: D u n +1 i + λ D u n +1 i − D u n − i t − u n +1 i − u n − i t = 0(A.3)with D u i = u i +1 − u i + u i − ∆ r + 1 r i u i +1 − u i − r − r i u i . (A.4)If we use matrix notation, we may write the N − i = 2 , ..., N − (cid:0) − δ ij + (cid:0) λ + 2∆ t (cid:1) D ij (cid:1) u n +1 j = (cid:0) − δ ij + λ D ij (cid:1) u n − j , (A.5)where we used Einstein’s convention for the sum of re-peated indexes, and D i,i − = 1∆ r − r i ∆ r (A.6) D i,i = − r − r i (A.7) D i,i +1 = 1∆ r + 12 r i ∆ r (A.8)and D ij = 0 otherwise.The other two equations are given from the boundaryconditions. For r = 1, we have σ rr = 0. The discretizedadimensional boundary condition becomes: Du n +1 i + λ Du n +1 i − Du n − i t = 0 (A.9)or in matricial notation: (cid:0) λ + 2∆ t (cid:1) D j u n +1 j = λ D j u n − j , (A.10)with D = − r + 2 r , (A.11) D = 4∆ r , (A.12)and D j = 0 for all other values of j . For r N = R ∞ , σ rr = σ . The discretized adimensional boundary condi-tion becomes: (cid:0) λ + 2∆ t (cid:1) D Nj u n +1 j = λ D Nj u n − j + 2∆ tσ, (A.13) with D N,N − = − r , (A.14) D N,N = 4∆ r + 2 r N , (A.15)and D Nj = 0 for all other values of j .In sum, we have a matricial equation of the kind: A ij u n +1 j = B ij u n − j + C i , (A.16)where the matrix A ij is A j = (cid:0) λ + 2∆ t (cid:1) D j , (A.17) A ij = − δ ij + (cid:0) λ + 2∆ t (cid:1) D ij , ( i = 2 , . . . , N − , (A.18) A Nj = (cid:0) λ + 2∆ t (cid:1) D Nj , (A.19) the matrix B ij is B j = λ D j , (A.20) B ij = − δ ij + λ D ij , ( i = 2 , . . . , N − , (A.21) B Nj = λ D Nj , (A.22)and the vector C i is given by C N = 2∆ tσ and C i = 0otherwise. [1] D. Weihs, A. Gefen, and F. J. Vermolen, Interface Focus , 20160038 (2016).[2] S. N. Jorgensen and J. R. Sanders, Med. Biol. Eng. Com-put. , 1297 (2016).[3] V. Ajeti, A. P. Tabatabai, A. J. Fleszar, M. F. Staddon,D. S. Seara, C. Suarez, M. S. Yousafzai, D. Bi, D. R.Kovar, S. Banerjee, et al., Nat. Phys. , 696 (2019).[4] A. Jacinto, A. Martinez-Arias, and P. Martin, Nat. CellBiol. , E117 (2001).[5] P. Martin and J. Lewis, Nature , 179 (1992).[6] R. J. Tetley, M. F. Staddon, D. Heller, A. Hoppe, S. Baner-jee, and Y. Mao, Nat. Phys. , 1195 (2019).[7] A. Brugu´es, E. Anon, V. Conte, J. H. Veldhuis, M. Gupta,J. Colombelli, J. J. Mu˜noz, G. W. Brodland, B. Ladoux,and X. Trepat, Nat. Phys. , 683 (2014).[8] E. Javierre, F. Vermolen, C. Vuik, and S. Van der Zwaag,J. Math. Biol. , 605 (2009).[9] R. T. Tranquillo and J. Murray, J. Surg. Res. , 233(1993).[10] D. G. Sami, H. H. Heiba, and A. Abdellatif, Wound Med. , 8 (2019).[11] S. A. Eming, T. A. Wynn, and P. Martin, Science ,1026 (2017).[12] B. A. Purnell and P. J. Hines, Science , 1020 (2017). [13] J. Reinke and H. Sorg, Eur. Surg. Res. , 35 (2012).[14] F. Huber, J. Schnauß, S. R¨onicke, P. Rauch, K. M¨uller,C. F¨utterer, and J. K¨as, Adv. Phys. , 1 (2013).[15] J. J. Lee, L. Talman, S. M. Peirce, and J. W. Holmes,Biomech. Model. Mechanobiol. , 1297 (2019).[16] L. Rold´an, J. J. Mu˜noz, and P. S´aez, Comput. MethodsAppl. Mech. Eng. , 28 (2019).[17] A. Guerra, J. Belinha, and R. N. Jorge, J. Theor. Biol. , 1 (2018).[18] B. A. Camley and W.-J. Rappel, J. Phys. D: Appl. Phys. , 113002 (2017).[19] F. J. Vermolen, in Encyclopedia of Cell Biology (AcademicPress, 2016), vol. 4, p. 117.[20] D. Tartarini and E. Mele, Front. Bioeng. Biotechnol. ,206 (2016).[21] R. O’Dea, H. Byrne, and S. Waters, in ComputationalModeling in Tissue Engineering (Springer, 2012), vol. 10,p. 229.[22] B. D. Cumming, D. McElwain, and Z. Upton, J. R. Soc.Interface , 19 (2010).[23] J. C. Arciero, Q. Mi, M. F. Branca, D. J. Hackam, andD. Swigon, Biophys. J. , 535 (2011).[24] J. C. Arciero, Q. Mi, M. Branca, D. Hackam, andD. Swigon, Wound Repair Regen , 256 (2013). [25] L. Geris et al., Computational modeling in tissue engi-neering (Springer, 2013).[26] L. Carvalho, P. Patricio, S. Ponte, C. P. Heisenberg,L. Almeida, A. S. Nunes, N. A. M. Ara´ujo, and A. Jacinto,J. Cell Biol. , 4267 (2018).[27] I. Bonnet, P. Marcq, F. Bosveld, L. Fetler, Y. Bella¨ıche, and F. Graner, J. R. Soc. Interface , 2614 (2012).[28] S. Tlili, C. Gay, F. Graner, P. Marcq, F. Molino, andP. Saramito, Eur. Phys. J. E , 1 (2015).[29] E. H. Dill, Continuum mechanics: elasticity, plasticity,viscoelasticity (CRC press, 2007).[30] L. Landau and E. Lifchitz,