Stability of a two-dimensional biomorphoelastic model for post-burn contraction

We consider the stability analysis of a two-dimensional model for post-burn contraction. The model is based on morphoelasticity for permanent deformations and combined with a chemical-biological model that incorporates cellular densities, collagen density, and the concentration of chemoattractants. We formulate stability conditions depending on the decay rate of signaling molecules for both the continuous partial differential equations-based problem and the (semi-)discrete representation. We analyze the difference and convergence between the resulting spatial eigenvalues from the continuous and semi-discrete problems.


Introduction
In recent decades, healthcare has made significant progress so that today patients can survive even severe burns. Nevertheless, burns are still life-threatening and difficult to assess and treat (Lang et al. 2019). These injuries will still significantly impact the quality of human life. Additional factors after a burn can include shock, infection, and prolonged stress and are of a physical, mental and social nature. Besides slow wound healing, the prevention of hypertrophic scars and contractures, which always cause a reduction in patient mobility, are significant challenges in burn treatment (Wang et al. 2018).
Almost all full-thickness burns lead to scarring. Deep burns leave no vital skin tissue, including the collagen network, cells (fibroblast), and vasculature. If the resulting scars entail contraction, these scars are subject to change in dimensions. A critical physiological aspect is the occurrence of burn contractures, which are contractions that cause reduced mobility of joints. Mostly, those scars require reconstructive surgery.
Burn healing comprises three partly overlapping phases: inflammation (reactive), proliferation (reparative), and maturation (remodeling). First, immune cells clear contaminants and pathogens within several hours after injury. Secreted growth factors stimulate cells to migrate to the wound from the intact peripheral dermis and subcutaneous tissue to the wound. This migration is a hallmark of proliferation; the cells multiply in the injured area and replace the fibrin network by regenerating collagen. During this phase, cells form a temporary spongy extracellular matrix (ECM) (granulation tissue) which is replaced by a solid matrix at a much later stage (remodeling). Collagen type III fills the granulation tissue as a provisional matrix early, and collagen type III is replaced by the embryonic collagen type I during remodeling. In this final stage, which can last for years, the scar matures and forms a balanced structure.
Post-burn contraction is one sub-process during the proliferative phase. Under the influence of growth factors, fibroblasts can differentiate into myofibroblasts (Desmoulière et al. 1993). Myofibroblasts produce a large amount of collagen to which the cells attach and pull forces. Furthermore, like fibroblasts, these cells stimulate both the production of the components of the new collagen-rich ECM and the release of matrix metalloproteinases (MMPs). Usually, myofibroblasts disappear by apoptosis when the wound closes (Desmoulière et al. 1995). However, if myofibroblasts persist in a closed wound, they continue to pull, produce excess collagen, and develop a hypertrophic scar (Tomasek et al. 2002). Thus, the biomechanical interaction of fibroblasts, myofibroblasts, growth factors, and collagen is essential in post-burn contraction and hypertrophy.
Various scientific disciplines study the prevention of contractures, including biology, the medical sciences, and mathematics. For a couple of decades, mathematical models have simulated wound healing. These models predict the behavior of experimental and clinical wounds and contractions (Tranquillo and Murray 1992;Olsen et al. 1995; Barocas and Tranquillo 1997;Dallon et al. 1999;McDougall et al. 2006;Koppenol 2017;Menon et al. 2017). We can divide most of these models into continuum hypothesis-based, discrete cell-based, and hybrid models (Koppenol 2017). Here, the discrete cell-based model type considers individual cell interactions compared to the continuum models considering cellular densities in a tissue. Hybrid modeling is classically defined as coupling a continuous approach with a discrete one to model a complex phenomenon that cannot be described in a standard homogeneous way mainly because of its inherent multiscale nature (Stéphanou and Volpert 2015). One subcategory of the continuum (partial differential equations-based) hypothesis-based models is the mechano-(bio)chemical model. This formalism and the hybrid model provided the basis for the biomorphoelastic model we use in this study, which Koppenol and Vermolen developed (Koppenol and Vermolen 2017). This biomorphoelastic model can simulate permanent post-burn contraction and reproduce the trends observed in real-life data.
Following up on our previous study (Egberts et al. 2021b), we analyze the stability of the biomorphoelastic model for post-burn contraction around equilibria in a two-dimensional setting. The model's equations are nonlinear and multivariable, making multiple steady states possible. Some steady states might be unstable for specific parameter values, and the model might reach other steady states for linearly unstable parameter values. Therefore, this study aims to avoid the parametric dependence of stable and unstable solutions and to assess the local behavior around the steady states. We analyze the system of equations and provide stability conditions using Gershgorin's theorem, where we cannot compute the eigenvalues exactly. As in Egberts et al. (2021b), we use a linear stability analysis with Fourier series, where the transformations represent perturbations around equilibria.
Further, we give the eigenvalues for a specific case. Biologically this means that the equilibria of the effective strain are determined. We distinguish between the continuous problem, which represents the solution to the system of partial differential equations, and the semi-discrete problem, which represents the solution of a semi-discrete solution method. In the latter case, the spatial finite difference method is carried out, whereas the time remains continuous. We show that the continuous system's stability implies the semi-discrete system's stability. Besides stability conditions, we pay attention to the effects of system instability on actual post-burn contraction. We also discuss specific parts of the model that can be adjusted to reproduce additional biological observations and the application of machine learning to make the application of the model for personalized healthcare possible. The results reported in this paper complement the stability analysis of the model's one-dimensional counterpart.
The structure of this paper is as follows. Section 2 presents the mathematical model, and Sect. 3 presents the stability analysis. Then, Sect. 4 presents the validation of (in)stability. Finally, Sect. 5 presents the discussion and conclusion. total deformation gradient into a deformation gradient by growth or shrinkage, and deformation by mechanical forces (Hall 2008). These processes may occur simultaneously or consecutively. The advantage of this model is that it can simulate permanent deformation resulting from the contraction process in burn wound healing. The most important variable in this model is the skin displacement (u), i.e., the variable that allows us to determine the area of the wound and, in later stages, the scar. This variable is estimated using a set of other variables that we can divide into biochemical and mechanical. The biochemical variables are the fibroblast density (N ), the myofibroblast density (M), the signaling molecules concentration (c), and the collagen density (ρ), and for the mechanics, we have the displacement velocity (v) and the effective strain (ε). These variables are each modeled by a partial differential equation (PDE). Hence, we consider continuous distributions and concentrations in contrast to agentbased modeling considering individual cells. In this section, we show and explain the equations that involve the material time derivative D(·) Dt and passive convection z(∇ · v), z ∈ {N , M, c, ρ, v}. These concepts are introduced because the domain of computation is subject to displacement.
After an injury, the cascade of processes that characterize wound healing is initiated by releasing growth factors and cytokines (by immune cells) that we consider a collective of signaling molecules. These signaling molecules influence cell proliferation (as an activator-inhibitor), myofibroblast differentiation, chemotaxis, and the synthesis and decay of collagen. For the cells, we consider migration towards the gradient of the signaling molecules (Postlethwaite et al. 1987) by a minimal model for chemotaxis (Hillen and Painter 2008), and cell density-dependent Fickian diffusion for the inclusion of random walk. The equations for the cell densities of fibroblasts and myofibroblasts also contain a logistic-like cell proliferation term. The equations are given by: Here, D F represents the (myo)fibroblast diffusion coefficient, which accounts for a random walk, and χ F is the chemotactic parameter, r F is the cell division rate, r max F is the maximum factor of cell division rate enhancement because of the presence of the signaling molecules, a I c is the concentration of the signaling molecules that cause half-maximum enhancement of the cell division rate, κ F (N + M) represents the reduction in the cell division rate because of crowding (Vande Berg et al. 1989), q is a fixed exponent, k F is the signaling molecule-dependent cell differentiation rate constant of fibroblasts into myofibroblasts, and δ N , δ M represent the apoptosis rates of the fibroblasts and myofibroblasts, respectively. An essential difference between these equations is that myofibroblasts are assumed to only proliferate in the presence of signaling molecules.
The following equations describe the evolution of the signaling molecules and collagen. The signaling molecules only migrate because of (Fickian) diffusion, and collagen is not subject to active migration. In both equations, (myo) fibroblasts are responsible for the secretion, and MMPs are responsible for the breakdown: Here D c is the Fickian diffusion coefficient of the signaling molecules, k c is the maximum net secretion rate of the signaling molecules, η I is the ratio of myofibroblasts to fibroblasts in the maximum secretion rate of the signaling molecules and collagen, a I I c is the concentration of the signaling molecules that causes the half-maximum net secretion rate of the signaling molecules, δ c is the proteolytic breakdown rate parameter of the signaling molecules, η I I is the ratio of myofibroblasts to fibroblasts in the secretion rate of the MMPs and 1 + a I I I c c represents the inhibition of the secretion of the MMPs. Further, k ρ is the collagen secretion rate, k max ρ is the maximum factor of secretion rate enhancement because of the presence of the signaling molecules, a I V c is the concentration of the signaling molecules that cause the half-maximum enhancement of the secretion rate of collagen and δ ρ is the degradation rate of collagen. A generic MMP affects the reaction kinetics of the signaling molecules and collagen, and we assume it is always at a local equilibrium concentration. This assumption avoids even more complexity and additional unknown and undocumented parameter values.
We capture the mechanics of the model by two PDEs for the displacement velocity (v) and the effective strain (ε). In the equation for displacement velocity, the Cauchy stress tensor σ relates to the effective strain and displacement velocity gradients by a visco-elastic constitutive relation. The body force f is generated by a pulling force on the ECM by myofibroblasts, which we assume to be proportional to the product of the cell density of the myofibroblasts and a function of the concentration of collagen 1 : This equation represents the balance of momentum. Here ρ t represents the total mass density of the dermal tissues, ξ is the generated stress per unit cell density and the inverse of the unit collagen concentration, and R is a constant. Despite many studies neglecting inertial effects, we have kept the inertia terms (the time derivative) closer to the underlying physics. From a mechanical point of view, we assume the tissue to be isotropic and homogeneous, except for a dependency of the stiffness on the local collagen density (Ramtani 2004). The visco-elastic relation for the Cauchy stress tensor is: where μ 1 , μ 2 are the shear and bulk viscosity, respectively, E √ ρ represents Young's modulus (stiffness), and ν is the Poisson's ratio. Despite possibly obtaining large deformations in the tissue, we use linear elasticity to avoid the need to include additional input parameters. We incorporate permanent deformation because of microstructural changes in the tissue via morphoelasticity. This tensor-based approach is also commonly used in the growth of tissues (such as tumors). For the constitutive law, we use Hooke's law which assumes a linear relationship between stress and strain. Following (Koppenol and Vermolen 2017), we assume the effective Eulerian strain to be modelled by Here, ζ is the rate of morphoelastic change (i.e., the rate at which the effective strain changes actively over time). In particular, the tensor for contraction depends on the product of the concentration of the MMPs, the concentration of the signaling molecules, and it is inversely proportional to the collagen density (note that the collagen density drops out because of the linear dependence of the equilibrium MMP concentration on the collagen density). Further, the linear strain evolution equation (7) provides an appropriate description provided that all components of ε are small. Therefore, as well as eqn. (6), if ε → O(1), then the model requires another constitutive law.

Further assumptions and boundary conditions
For the two-dimensional setting, we locate the x y-plane parallel to the surface of the skin and Variations over skin depth are disregarded; hence the computations are done on an arbitrary skin depth. The domain, and the initial conditions, are symmetrical; hence the solution inherits this property. Therefore, we can perform calculations on a reduced domain to benefit from a computational workload.
We define the computational domain's boundary by x,t = o x,t . Here o represents the outer nonsymmetrical boundaries, h represents the horizontal symmetrical boundary where y = 0, and v represents the vertical symmetrical boundary where x = 0. For the constituents of the dermal layer, the following boundary conditions hold for all time t and for all x ∈ p x,t : where p ∈ {h, v} and n is the outward pointing normal vector. We use similar conditions for the mechanical part of the model, that is, for all time t and for all x ∈ p x,t : v · n = 0 and (σ · n) · τ = 0, where τ is the tangential vector. It is not allowed to specify boundary conditions for ρ and ε because of overdetermination since the equations for ρ and ε are ordinary differential equations with respect to time t. We vary the initial conditions around N , 0, 0, ρ, 0 for N (x, 0), M(x, 0), c(x, 0), ρ(x, 0), and v 1 (x, 0), respectively, and we set v 2 (x, 0) = − v 1 (x, 0) to ensure symmetry. Further, we use ε(x, 0) = 0.

Linear stability of the model
We analyze the stability of the two-dimensional biomorphoelastic model for postburn contraction to understand the a priori behavior of the solution. We cannot derive the (exact) solution to coupled nonlinear partial differential equations system. In line with our previous stability analysis in R 1 , we first analyze the linear stability of the continuous problem. Further, we also analyze the stability of the numerical approximation. We consider the following linearised equations around equilibria (N , M, c, ρ, v 1 , v 2 , ε 11 , ε 12 , ε 22 ) = (N , 0, 0, ρ, 0, 0, ε 11 , ε 12 , ε 22 ), where N , ρ ∈ R ≥0 and ε 11 , ε 12 , ε 22 ∈ R: for the chemical part of the model, where we used that k ρ = δ ρ ρ 2 must hold in equilibrium, and for v 1 (the equation for v 2 is similar, where x, y andv 1 ,v 2 are interchanged), and for the effective strainsε 11 ,ε 12 (the equation forε 22 is similar as forε 11 , where x, y andv 1 ,v 2 are interchanged). In Eqs. (13)-(15),N ,M,ĉ,ρ,v 1 ,v 2 ,ε 11 ,ε 12 , andε 22 are variations around the equilibria. Hence, N (x, t) = N +N (x, t), etc. Note that we only consider the equation for ε 12 and not ε 21 . We demonstrate that if the strain tensor ε is initially symmetric, then it remains symmetric at all later times (Egberts et al. 2021a). , Note that we used sym(∇v) T = sym(∇v) and skw(∇v) T = −skw(∇v). Subtraction gives From the above equation, it is clear that (ε − ε T ) = 0 represents an equilibrium, and hence symmetry of ε represents an equilibrium. Hence, we conclude that initial symmetry implies no symmetry changes for later times. Furthermore, we also prove that ε − ε T is the only solution if ε − ε T = 0 at t = 0.
Performing the matrix scalar product A : B := i, j A i j B i j on the above equation by ε − ε T gives upon setting w = ε − ε T and M = ∇v: Hence we obtain Define ||w|| 2 := w : w, then it follows that Integrating over t from t = 0 and using w = 0 at t = 0, gives 0 ≤ ||w|| 2 = −ζ With ζ, N , η I I , M, c, a I I I c ≥ 0, this implies that ||w|| = 0 on t > 0 if ||w|| = 0 on t = 0. Hence w = 0 for t > 0, which represents symmetry, is the only possibility if w = 0 on t = 0. (7) depends on a linear relationship between stress and strain, hence ε − ε T = 0 is a solution. However, this solution is not guaranteed to be unique; therefore, initial symmetry may change over time because of computing and rounding errors, for instance.

Remark 1 Equation
Theorem 1 motivates why we only need to consider ε 12 as a cross term assuming initial symmetry. Further, we demonstrate that small perturbations around symmetry of ε remain small, which is a characteristic of stability.

Theorem 2 Let Eq. (7) hold on an open Lipschitz domain
for t > 0. Let ε be symmetric for t ≥ 0, then stability of symmetry is warranted if and only if K ≥ 0.
a system of ordinary differential equations. Writing out, then D Dt Hence In matrix-vector form, let ω = [ω 11 , ω 12 , ω 21 , ω 22 ] , then we get For K = 0 (in equilibrium, c = 0 in Eq. (17)), this matrix is skew-symmetric (that is B T = −B), and hence the eigenvalues are zero or purely imaginary. This implies that ω = 0 ⇔ ε = ε T is a null-stable equilibrium. Hence, small perturbations around the symmetry of ε will remain small, which implies stability. For the case that K > 0, it follows that the real part of the eigenvalues is given by K , which also gives A-stability (perturbations from symmetry vanish as t −→ ∞). For K < 0, which corresponds to expansion (instead of contraction in our model), the negative real part of the eigenvalues results in instability of symmetry. Although the current case is not similar to diffusional growth, it is known that diffusional growth in combination with surface processes suffers from Mullins-Sekerka instabilities (Caroli et al. 1986), which exhibits growth of perturbations on spherical surfaces.

Remark 2
Of course v is non-constant. The only thing that happens is that v = v(t) impacts the angular frequency around the equilibrium ε = ε T .

Stability of the continuous problem
We write the variations around the equilibria in terms of a complex Fourier series, for z ∈ {N ,M,ĉ,ρ,v 1 ,v 2 ,ε 11 ,ε 12 ,ε 22 }, where | | denotes the measure of and i represents the unit imaginary number. Substitution of the variations (29) into the linearised equations (13)-(15), multiplication by e −2iπlx e −2iπ py , integration over = (0, 1) 2 (| |=1) and double orthonormality over giveṡ for the chemical part of the model, and for the displacement velocity, anḋ for the effective strain.
Interchanging the third and first equation of (30), these equations together with Eqs. (31) and (32) are in the form y + Ay = 0 (y the vector of the time-dependent components) with We determine the eigenvalues of A by solving |A − λI | = 0 for λ, where I represents the identity matrix. First, we perform Gaussian elimination to see that we can use the first four diagonal values as pivots. Hence, the first four eigenvalues are the first four diagonal entries. The system is linearly stable if and only if the real part of all the eigenvalues is non-negative, hence stability for the chemical part of the model is guaranteed if: These four requirements show that stability for the chemical part of the model is equal to the stability constraints in R 1 (Egberts et al. 2021b). That is, for δ c ≥ k c a I I c ρ and qδ N ≤ κ F r F N 1+q (l = p = 0). We note we need δ N > 0 and hence κ F N < 1. The second and fourth eigenvalues meet the stability condition independent of the chosen parameter values, given that these values are positive. Hence, if the conditions are met for l = p = 0, they hold for all l, p ∈ Z, corresponding to wavelike perturbations. Further, we end up with a 5 × 5-submatrix containing the mechanical part of the model. For this 5 × 5-matrix, we see that the last three columns contain possibly nonzero values at the first two row positions; hence, these columns are linearly dependent. From this, it immediately follows that λ = 0 is an eigenvalue.
Applying Gershgorin's Theorem leads to eigenvalues that can be located anywhere in a union of circles centered around the origin. Hence, Gershgorin's Theorem does not exclude eigenvalues with a negative real part (which could reside in the left half of the complex plane). For this reason, linear stability for the effective strain is guaranteed if |α|(ε 11 + ε 12 + ε 22 − 1) − |β|ε 12 = 0, |α|(ε 11 − 1 for all α, β ∈ Z. Hence, linear stability for wavelike perturbations around equilibria (α, β ∈ Z =0 ) is certainly obtained for ε 11 = 1 2 , ε 22 = 1 2 and ε 12 = 0. This implies that A i, j = 0 for i ∈ {7, 8, 9}, j ∈ {5, 6}. We note that with these equilibria, we have ε → O(1), hence describing a physical situation in which the model can no longer be applied. In this case, Gershgorin's theorem cannot be used to access the stability criteria, and there are no other strategies to solve the fourth-order polynomial analytically. In order to 'show' that we have stability for our set of parameter values, we provide an empirical argument based on the numerical approximation of the eigenvalues.
The remaining eigenvalues follow from the 5 × 5-matrix. The eigenvalue λ = 0 has algebraic multiplicity 3. The other two eigenvalues follow from the upper left 2 × 2-block matrix and are in addition to that determined by We note that A 56 = A 65 . Hence, the remaining eigenvalues are real-valued. Solving the above equation with the abc-formula then gives Here, the discriminant is always non-negative. For stability, in this case, a necessary condition is that Squaring the left-and the right-hand-side gives Substitution of A 55 , A 66 and A 56 gives which reduces to Hence, for all μ 1 , μ 2 ≥ 0 and all l, p ∈ Z, the stability constraint is satisfied. We summarize these results in Theorem 3.
Remark 3 Note that δ c ≥ k c a I I c ρ , for k = 0 (constant states). Hence, if constant perturbations are stable, then wavelike perturbations are stable. In case δ c is not large enough, fast oscillating perturbations will vanish, while slow oscillating perturbations will not vanish and can amplify. Further, the mathematical model is actually not suitable for ε 11 = ε 22 = 0.5; however, this is still a consequence of the above analysis.
For the empirical 'proof' of the stability constraints, we only consider the eigenvalues λ of the 5 × 5-submatrix of matrix A.We keep the parameter values asinTable2 and vary ε 11 and ε 22 between − 1 and 1 with stepsize 0.01, ε 12 between − 1 and 0.5 with stepsize 0.1, and l and p (integers) between 1 and 100. We define Hence, S = 1 corresponds to stability, whereas S = 0 corresponds to instability. Figure 1 shows the results of S for some values of ε 12 .
If ε 12 increases, then the region for stable ε 11 and ε 22 values grows. Further, there is symmetry in the line y = x. Given the complexity of the fourth-order polynomial, it is unclear what defines these boundaries of stable regions. One can predict that there will not be any stable values for the effective strain in the shown domain if ε 12 becomes (much) smaller. However, one does not expect the effective strain to become that negative (the model breaks down for these values). Focusing on the origin (ε 11 , ε 22 ) = (0, 0), one can see that for small |ε 12 |, say |ε 12 | ≤ 0.2, there is always a stable region around the origin. For values of |ε 12 | that are larger, which are not interesting from a physical point of view, the stability region may no longer include the origin.
where the model is linearly stable.

Stability of the discrete problem
The stability of the continuous problem does not always automatically imply the stability of the discrete counterpart of the problem. Therefore, we assess the stability of the semi-discrete problem, which we use to analyze the stability of the whole discrete system. Lax's Equivalence Theorem states that a consistent, stable method converges. Under stability, the global truncation error tends to zero as the step size tends to zero (as h → 0) if the local truncation error (i.e., the difference between the derivatives and difference ratios) tends to zero as the step size goes to zero. We consider a unit rectangle = [0, 1] 2 that we divide into small rectangles with sides h such that (n + 1)h = 1. At the intersections of the grid lines, we have nodal points where we approximate the solutions of the variables. We denote the unknowns at node (x i , y j ) by z i, j , z ∈ {N , M, c, ρ, v 1 , v 2 , ε 11 , ε 12 , ε 22 } and apply finite differences on the eigenvalue problem. Then the finite difference method (FDM) of the spatial part of the linearized problem (13)-(15) gives: for the chemical part of the model, and writing v 1 i, j for v 1i, j and ε 11 i, j for ε 11i, j etc. gives for v 1 (again, the equation for v 2 is similar), and for the effective strains ε 11 , ε 12 (again, the equation for ε 22 is similar). We perform von Neumann eigenvalue and stability analysis. Let for z ∈ {N , M, c, ρ, v 1 , v 2 , ε 11 , ε 12 , ε 22 }, where i again represents the imaginary unit number. Substitution of (47) in Eqs. (44)- (46), multiplication by e 2πlkhi e 2π pjhi , using Euler's formula, 2 − 2 cos(2π x) = 4 sin 2 (π x) and e i x e iy − e i x e −iy − e −i x e iy for the chemical part of the model and for the displacement velocity v 1 (the equation for v 2 is similar and yields an equal result), and for the effective strains ε 11 and ε 12 (the equation for ε 22 is similar and yields an equal result for ε 22 ). As an example, the derivation of Eq. (49) is given in the Appendix. The Eqs. (48)-(50) are in the form λz = Cz with the matrix C as in (33). Hence, the eigenvalues are found the same way as in the continuous case. Note that, since the overall system has the form y + Ay = 0, the discrete system is linearly stable if and only if the real part of the eigenvalues is non-negative, hence we need: for the chemical part of the model. To guarantee linear stability, the first requirement of Eq. (51) states qδ N ≤ κ F r F N 1+q , given δ N = r F (1 − κ F N )N q . The second and fourth eigenvalues meet the stability condition independent of the chosen values for the parameters, given that the parameters are positive. Finally, the third requirement states δ c ρ ≥ k c a I I c . These statements remain the same when the horizontal and vertical step sizes are unequal ( x = y).
Here, we used sin(2π xh) = 2 sin(π xh) cos(π xh). Note that the subtractions by the cosines are bounded from above. Therefore, for all μ 1 , μ 2 and all l, p ∈ Z, Eq. (53) is satisfied. To conclude, we have demonstrated that if the equilibrium is stable in the continuous problem, it is also stable in the semi-discrete problem. A consistency exists between the continuous problem's stability criteria and the discrete problem's stability criteria. We show this by writing sin(x) and sin 2 (x) as a Taylor series. Substitution into the third equation in (51) yields: Comparison to the third equation in (34)

Numerical validation
We need to validate whether the linear stability conditions we have derived also hold in a finite element setting where we consider the fully nonlinearly coupled model. We approximate the solution to the model equations by the finite element method using bilinear basis functions (Van Kan et al. 2014). We refer to Appendix 2 in Koppenol and Vermolen (2017) for deriving the general finite-element approximation. We consider a rectangle = [0, L] 2 that we subdivide into small rectangles (quadrilaterals) with sides x = y. Then, we convert the regular mesh to a triangulation. The quadrilateral mesh faces are converted to triangles by splitting the faces into triangles according to a cross-division of the quadrilateral. Let ≈ h = n e e=1 e the global finite element subspace and a j , j ∈ {1, . . . , n}, n ∈ N the coordinates of these vertices of the elements. We choose the piecewise linear Lagrangian basis functions ϕ i ∈ X h (t) with ϕ i (a j , t) = δ i j , i, j ∈ {1, . . . , n} as basis functions for the finite-dimensional subspace h , where δ i j denotes the Kronecker delta function. Note that the following holds for the chosen subspace h ⊂ x,t : Dϕ i Dt = 0 for all ϕ i (Dziuk and Elliot 2007). We simplify the Galerkin equations using this property. We approximate the integrals over the interior of the elements by a Newton-Cotes rule based on linear basis functions.
Regarding the computation of the integrals that arise during the discretization, we apply the Newton-Cotes integration rule and backward Euler time integration. We use a monolithic approach with inner Picard iterations to account for the non-linearity of the model equations. We approximate the local displacements of the dermal layer (u) with Further, we update the mesh and keep track of the mesh quality by computing min e k J e k / max e k J e k , e k ∈ , The reference is the solution in which h = 0.05 with J the Jacobian. In case min e k J e k / max e k J e k < 0.5, we perform remeshing. In all simulations, this ratio was at least 0.9915; remeshing was unnecessary. We use mass lumping to avoid loss of monotonicity (i.e., oscillations). We experimentally evaluate the convergence of the numerical method in a domain of [0, 3.2] 2 cm 2 with a wound between [0, 1.2] cm 2 . This domain represents a quarter of the domain of the modeled skin on which we perform computations, which is possible because of the model's symmetry. The transition from healthy to injured skin is steep, and we account for this steepness of gradients through an interval of 0.8 cm. In this transition, the initial solutions vary between the equilibria and the initial wound densities. Within the wound, we assume that 2000 fibroblast cells/cm 3 , 10 −8 g/cm 3 signaling molecules, and 0.01125 g/cm 3 collagen are present. We model the slope of the variables with sine functions.
We divide the computation domain into n x × n y = (3.2/h) 2 elements with h ∈ {0.2, 0.1, 0.05}. In order to have mesh convergence of the numerical solution of the second order, we choose t = h 2 . We simulate post-burn contraction for one day and report the densities of the variables (the solutions). We compute the convergence order results using the L 2 error norm. Let lim h→0 z h (x, 1) = z(x, 1) denote the true density of variable z on day one and z 0.05 (x, 1) =: z h ref the solution in the last simulation (i.e., the reference, the one computed using the highest numerical resolution). We approximate the errors with the following error definition: where the grid-points x i,n correspond to the grid-points in the simulation with h = 0.2 (n = 289 nodes). Hence, we evaluate the solution to the equations between simulations on a fixed set of initial nodes. Table 1 shows the results for the L 2 error. All the L 2 errors decrease consistently as h becomes smaller (figures not shown), and the values in the table show an order of convergence above O(h 2 ).
We perturb the initial conditions around equilibria using sine functions and vary the parameter δ c . Hence, we fix all parameters except for the signaling molecule decay rate. Table 2 shows the fixed parameter values. For the time integration, we use an initial step of t = 10 −2 days until half a day is simulated, after which we increase the time step by t new = min{2, 1.1 × t}. We use a domain of [0, 1] 2 cm 2 that we first divide into equilateral elements (rectangles) with h = 0.05; then, we convert this mesh to an equilateral triangulation by cross division of quadrilaterals.
For the initial conditions, we vary the number of waves k using two levels (1 and 2). We perturb the initial condition for the fibroblasts and collagen by using a product of sine functions with an amplitude of 10 cells/cm 3 and 10 −2 g/cm 3 , respectively. Since the fibroblasts' and collagen's biological equilibria are non-zero, we do not need to worry about negative values in these initial conditions. For the initial condition of the myofibroblasts and the signaling molecules, we use a product of uniform splines with 2k + 1 knots. On the boundaries, the knots have zero value, and between the values are 3 and 1 cells/cm 3 for the myofibroblasts and 2 × 10 −7 and 5 × 10 −8 g/cm 3 for the signaling molecules. These setups ensure that the myofibroblast distribution and the signaling molecule density values are positive. We note that these initial amplitudes become smaller because of the product of splines. The initial amplitudes of the displacement velocity v 1 is 5 × 10 −1 . To ensure symmetry, we set the initial condition v 2 (x, 0) = −v 1 (x, 0). We do not perturb the effective strain densities and set ε 11 (x, 0) = ε 12 (x, 0) = ε 22 (x, 0) = 0.
For stability, Theorem 3 further requires that δ c ρ ≥ k c a I I c in case l = 0. We choose to vary the signaling molecule decay rate δ c using three levels: 2 × 10 −4 , 3 × 10 −4 and 5 × 10 −4 cm 6 /(cells g day), where the first two values do not meet the stability condition.
In the first simulation, we take δ c = 5 × 10 −4 cm 6 /(cells g day), k = 1 and simulate over a time interval of 200 days. We note that this parameter value meets the stability criterium. Figures 2, 3 and 4 show the results.
The displacement velocity densities move to equilibria at a rapid pace. In the first time integration, the densities drop two orders from 5 × 10 −1 to 5 × 10 −3 within 0.01 day (≈ 15 min). The plots show that the peaks in the regions (x, y) = {0 ≤ x ≤ 0.5 ≤ y ≤ 1 ∧ 0 ≤ y ≤ 0.5 ≤ x ≤ 1} merge (i.e., the top left and bottom right corners, negative values for v 1 and positive values for v 2 ). The peaks in the regions (x, y) = {0 ≤ x, y ≤ 0.5} (i.e., the bottom left corners) shift to the edges of the computational domain, the one in v 1 shifting to the horizontal axis of symmetry, and the one in v 2 shifting to the vertical axis of symmetry. These variations in the displacement velocities are still visible on day 0.02. The peaks in the regions (x, y) = {0.5 ≤ x, y ≤ 1} (i.e., the top right corners) disappear within this time. The densities drop one order further towards equilibria in the next quarter of an hour. Following the stability theory, the displacement velocity densities converge to equilibria within ten days. Figure 3 shows the evolution of the effective strains for the values of input parameters within the stability regime. The effective strain densities change from the equilibria to perturbations on the first day of the simulation because of the initial perturbations in the other variable densities. We see that diagonal tensions arise in the effective strain densities ε 11 (x, 1) and ε 22 (x, 1). At the same time, it is more circular for the effective strain density ε 12 (x, 1). Here, positive and negative values alternate. These diagonal tensions disappear gradually between day ten and day 200 as the effective strain densities move to the equilibria ε 11 (x, 200) = ε 12 (x, 200) = ε 22 (x, 200) = 0. Note that the theory (Thrm. 3 part 1) states that the constant state equilibria ε 11 , ε 12 , ε 22 ∈ R are stable.  Table 2 shows the values of the other parameters. In subfigure a the upper plots show the displacement velocity v 1 , and the lower plots show the displacement velocity v 2 . The shown domains are (0, 1) 2 cm 2 , and the color bars show the displacement velocities in cm/day. For both v 1 and v 2 , the first two plots have different color bars, and the last three plots share the same color bar shown on the right. In subfigure b the left and right plots show the time-dependent change of v 1 at points p 1 ≈ (0.5, 0.5) and p 2 ≈ (0.75, 0.75), respectively (see the black stars in the last plot for v 1 in (a). Note that these points are subject to displacement Figure 4 shows the evolution of the chemicals for the values of input parameters within the stability regime. All the plots show that the perturbed chemical densities move gradually toward the equilibria. For the fibroblast cell density, the perturbation leaves a few fibroblasts in the origin of the computational domain on day 25. It takes up to 200 days to move the cell density toward the equilibrium N = 10 4 cells/cm 3 . For the myofibroblast cell density, no cells are present on day 120 as the perturbed cell density moves quickly and gradually toward the biological equilibrium M = 0 cells/cm 3 . The perturbed signaling molecule density moves even quicker toward the biological equilibrium c = 0 g/cm 3 , having the perturbations almost vanish on day 25. The perturbed collagen density takes longer to move to equilibrium as it takes more than 25 days, and on day 200, a slight perturbation is still visible. From this figure, we can conclude that the signaling molecule density first moves toward the biological equilibrium, after which the myofibroblast cell density moves toward the biological equilibrium. It takes longer for the fibroblast cell and collagen density to move toward the biological equilibria, taking more time for collagen for the current input values. Overall, the perturbations disappear gradually with stable parameter values, and the numerical method behaves stably.
From a biological perspective, minor fibroblast cell and collagen variations already initialize long-term wound healing. Fibroblasts move toward the center of the wound,  Table 2 shows the values of the other parameters. In subfigure a the upper, middle, and lower plots show the effective strains ε 11 , ε 12 , and ε 22 , respectively. The shown domains are (0, 1) 2 cm 2 . The color bars show the effective strains (no unit) and correspond to the four plots shown on the left of the color bars. In subfigure b the left and right plots show the time-dependent change of ε 11 at point p 1 ≈ (0.5, 0.75) and of ε 12 at point p 2 ≈ (0.25, 0.75), respectively (see the black stars in the last plots for ε 11 and ε 12 in (a) and collagen regeneration takes over half a year. Variations arise in the effective strain, after which the tensions disappear.
In the next simulation, we take δ c = 3 × 10 −4 cm 6 /(cells g day) with k = 2. This parameter value does not meet the stability criterium. We only present a few figures to avoid too many figures in this manuscript.
Likely, as in Fig. 2, the perturbed displacement velocity densities move gradually toward equilibria in this simulation (figure not shown). In the first step of time integration, the densities align in a similar pattern as we have seen before, decreasing from  Table 2 shows the values of the other parameters. In subfigure a the upper, middle, and lower plots show the effective strains ε 11 , ε 12 , and ε 22 , respectively. The shown domains are (0, 1) 2 cm 2 , and the color bars show the effective strains. In subfigure b the left and right plots show the time-dependent change of ε 11 and ε 12 , respectively, at point p ≈ (1, 0.5) (see the black star in the last plot for ε 11 in (a). Here, the (initial) y-limit is adjusted from order 10 −3 to 10 −5 to show the effects order O(10 −1 ) to O(10 −3 ). However, in this simulation, in the next time step, the densities decrease to order O(10 −5 ) in contrast to the order O(10 −4 ) in Fig. 2. We see a different evolution for the effective strains and the chemicals. We distinguish between evolution in the simulation's early and later stages of the simulation. We see a similar pattern in the early stage of the simulation for the effective strains, as shown in Fig. 3. For the effective strains ε 11 and ε 22 , the tensions in the densities are diagonal and peak on day 5. The peaks diminish in magnitude in the first 51 days (figure not shown). The effective strain ε 12 density shows the same intensity of variations on day 5, which are circular and alternating between positive and negative values, like in the stable simulation. However, Fig. 5 shows this pattern changes after 51 days.
Unlike the equilibria found in the simulation with stable parameter values, the effective strain densities increase intensely in variation up to day 223 (the second  . 6 Evolution of the signaling molecule density for δ c = 3 × 10 −4 cm 6 /(cells g day) in the first 13 days. Table 2 shows the values of the other parameters. The shown domains are (0, 1) 2 cm 2 , and the color bar shows the signaling molecules in g/cm 3 column plots). For the effective strain ε 11 , this is an increase on the right edge of the computational domain; for the effective strain ε 22 on the top edge, and the effective strain ε 12 around the top right corner. In opposite directions, these densities decreased. For example, we see a decrease on the vertical axis of symmetry for the effective strain ε 11 . After these peaks of intensities in the effective strain densities on day 223, the densities gradually decrease until day 351 and increase in intensity until they reach equilibrium on day 1200. Although it is difficult to see, in this 2D figure, the effective strain densities oscillate around the (new) equilibria. Compared to the simulation with stable parameter values, we see an increase in the intensity of the same order for ε 12 , albeit with more significant numbers. We note that the order of magnitude may also result from the larger wavenumber (k = 2) in the initial perturbations. The early evolution of the chemicals for δ c = 3 × 10 −4 cm 6 /(cells g day) is comparable to the evolution of the chemicals for δ c = 5 × 10 −4 cm 6 /(cells g day) (figure not shown). In the first 51 days, the perturbed fibroblast and myofibroblast cell densities move gradually toward the biological equilibria. The perturbed collagen density moves gradually toward the biological equilibrium in the first 119 days. However, the perturbed signaling molecule density does not reach the expected biological equilibrium. Figure 6 shows the early evolution of the signaling molecule density.
Unlike the perturbed signaling molecule density evolution for stable parameter values shown in Fig. 4, the perturbations in the signaling molecule density do not disappear in the first 13 days for an unstable signaling molecule decay rate. The initial peaks of about 4 × 10 −14 cm 6 /(cells g day) decrease in the first few days. At the same time, these peaks merge and shift toward the origin as they decrease further in the first five days. The peaks continue merging, completed within 13 days; however, the signaling molecule density increases strongly in the origin of the computational domain. In the beginning, this increase does not significantly affect the other chemicals; however, after day 51, it causes a considerable difference. Figure 7 shows the evolution of the chemicals in the later stage of the simulation for δ c = 3 × 10 −4 cm 6 /(cells g day).
On day 51, it seems that the fibroblast cell density is in the equilibrium N = 10 4 cells/cm 3 , the myofibroblast cell density in the equilibrium M = 0 cells/cm 3 , the signaling molecule density in the equilibrium c = 0 g/cm 3 , and the collagen density is around the equilibrium ρ = 0.1125 g/cm 3 . However, these densities do not stay in and around equilibria. Note the orders of the signaling molecule concentration: 10 −10 g/cm 3 on day 200 compared to the order 10 −14 g/cm 3 on day 13 (see Fig. 6). After day 13, in the origin of the computational domain, the signaling molecule density increases  (2)), stimulate the production of collagen, and inhibit the decay of collagen (Eq. (4)). Further, myofibroblasts also stimulate collagen production. The myofibroblast cell density reaches a maximum on day 213, the collagen density on day 235, and the fibroblast cell density on day 245. After the signaling molecule density reaches a minimum on day 317, we see that the myofibroblast cell density reaches a minimum on day 337, the collagen density on day 375, and the fibroblast cell density on day 448. After these days, such an oscillating effect around new equilibria is visible, which converges on day 1200. The result is a permanently reduced number of fibroblasts, a permanently increased number of myofibroblasts, and a permanently elevated concentration of signaling molecules and collagen at the origin of the computational domain (i.e., the center of the burn). Taken together, with an unstable signaling molecule decay rate not too low, the numerical method initially behaves like a stable regime. This stable behavior changes at a later stage of simulation time, where the numerical method behaves stable enough to let the chemical densities reach new equilibria in an oscillatory way. From a biological perspective, an increased expression of signaling molecules (because of their reduced decay) can lead to a period in which a wound fluctuates in contraction. This contraction fluctuation is because the number of (migrating) myofibroblasts increases and decreases. In the beginning, the wound can heal well. However, because of the continued signaling, the scar will fluctuate in thickness and stiffness because of the present collagen concentration. The scar is also highly subject to contraction because of the abundance of myofibroblasts present. The abundance of myofibroblasts and the increased collagen concentration may signify hypertrophy. Figure 8 shows an example of the contraction intensity as a function of relative wound/scar area over time.
We may associate the excessive collagen deposition with keloids, and hypertrophic scars (Tuan and Nichter 1998). Abnormal TGF-β signaling in myofibroblasts is associated with the formation of hypertrophic scars (Zhang et al. 2020). Given our study, it is likely that such a situation arises because of a lower decay rate of signaling molecules. Furthermore, hypertrophic scars develop 1 to 2 months after injury, while keloids develop months to years after the initial injury. This period is consistent with our simulations showing that the abundance of myofibroblasts and collagen occurs after a few months, while the increased expression of signaling molecules occurs within a few weeks. Furthermore, experiments suggest that the hyper-proliferation of fibroblasts in hypertrophic scars can be reversed once their stimulation, such as the abundance of growth factors and cytokines, is abolished (Tuan and Nichter 1998). Our simulation partially reflected this. When we turn off the stimulation of signaling molecule expression at a later stage by setting the density to the biological equilibrium, we see the (myo)fibroblast cell densities and the collagen density change. The myofibroblasts seem to disappear, the collagen density recovers, and the fibroblast cell density recovers. However, with three cells/cm 3 myofibroblasts left in the center of the scar after 47 days after this reset, the fibroblast cell density does not increase above 9865 cells/cm 3 , and the collagen density does not go below 0.1234 g/cm 3 . The numerical method does not converge and decreases the time step. We, therefore, set the myofibroblast cell density to the biological equilibrium on day 1247 and see the collagen density move to equilibrium within 411 days and the fibroblast cell density within 436 days. Thus, according to this simulation, restoring the fibroblast cell and collagen density is possible when the signaling molecules and myofibroblasts disappear. Though, it then still takes over a year to repair the defects. Hence, given that the overexpression of signaling molecules occurs in the first weeks, we recommend monitoring this expression to intervene early when necessary.
In the last simulation, we take δ c = 2 × 10 −4 cm 6 /(cells g day) and k = 2. We note that this parameter value does not meet the stability criterium. While running the simulation, we see that in the first 40 days, the perturbed fibroblast cell density moves to equilibrium. After 40 days, the fibroblasts persistently differentiate into myofibroblasts, increasing the myofibroblast cell density over the whole computational domain from day 50. This differentiation happens because the signaling molecule density increases, and therefore, also the collagen density increases. Unlike in the last simulation, the displacement velocity does not vanish. For that reason, remeshing is necessary around day 80 of the simulation. Within a few days, remeshing is necessary again, and at some point, the Picard iterations do not yield convergence anymore because of continuous remeshing. Therefore, we ended this simulation. We see the same for k = 1.
From a biological point of view, the human body protects against lowering the signaling molecule decay rate to this extent to prevent such a non-realistic occurrence where collagen will cause the tissue to rupture because of excessive production.
To conclude this section, the two-dimensional morphoelastic model for skin contraction is stable under the condition that the signaling molecule decay rate is not reduced too far below the limit δ c ≥ k c a I I c ρ .

Discussion and conclusion
For the discussion, we focus on adapting the model to bring it closer to reality and on applying neural networks to bring the model into practice.
Since hypertrophy is one common complication after burn injuries, it would be interesting to incorporate hypertrophy into this morphoelastic model. Koppenol made some progress in incorporating hypertrophy ), though this model considers a compressible neo-Hookean solid instead of a morphoelastic solid. Therefore, these models need to be merged so that a three-dimensional morphoelastic model simulates contraction and hypertrophy. In addition, hypertrophy depends highly on angiogenesis, which can also be considered. References (Valero et al. 2012; Vermolen and Javierre 2011) provide a good start for this.
In a three-dimensional setting, computation times increase rapidly. In order to decrease computation, the computational domain could be reduced by only incorporating the wound or scar region. If so, we need to adjust the boundary conditions and compare the results to the results of the whole computation domain where the boundaries are 'sufficiently far away' from the wound boundaries. We could model the wound boundaries as elastic springs. Furthermore, rotational symmetry and isogeometric analysis offer solutions to the curse of dimensionality and the decreasing quality of the moving mesh.
Further, post-burn contraction and the persistence of myofibroblasts in hypertrophic scars depend on stretching, pulling, and tissue stiffness. If we incorporate these features into the model, we can study pulling and stretching forces because of children's growth and as a consequence of the joints' daily (repetitive) motion. We can incorporate these forces in the body force in Eq. (5) and in the boundary conditions.
As far as we know, it is uncertain if myofibroblasts proliferate. At least, we know myofibroblasts proliferate at a lower rate than fibroblasts. There is a difference between proto-myofibroblasts and fully differentiated myofibroblasts. One might discuss that proto-myofibroblasts either proliferate or contract the tissue and that fully differentiated myofibroblasts only contract the tissue. Although it may not have been demonstrated that fully differentiated myofibroblasts do not proliferate, we could nevertheless assume that these cells only contract if we, for example, add an equation for the proto-myofibroblasts that do proliferate. We can also easily consider that myofibroblasts, in response to TGF-β, move more slowly than fibroblasts do (Thampatty and Wang 2006) and that in vitro myofibroblasts can differentiate back to fibroblasts under the influence of Prostaglandin E2 (Garrison et al. 2013). Myofibroblast dedifferentiation has not been observed in vivo; therefore, this has not been accounted for in the current model formulation. The impact of the inclusion of this feature needs to be quantified.
The state change between fibroblasts and myofibroblasts is more gradual than has been modeled here. The intermediate state of proto-myofibroblasts could also be implemented into the modeling framework. If one seeks to incorporate this intermediate state, then it is important to know the behavior of proto-myofibroblasts. This model extension also needs additional input parameters, which could be prone to uncertainty. The current model allows the modification of the constitutive law between stress and strain. One could revise the elasticity part, which is currently modeled as linear (Hookean) and hyperelastic, such that large strains are allowed from a physical perspective. This extension adapts Eq. (7) (Hall 2008;Koppenol 2017), adding more nonlinearity to the model and making the stability analysis even more complicated.
Furthermore, this nonlinearity will also add some complexity to the numerical algorithms. Another extension concerns the implementation of different collagen types. This addition could make sense since myofibroblasts produce collagen type III in large quantities, whereas fibroblasts produce collagen I (embryonic). Collagen III is known to be aligned anisotropically (directed), whereas collagen I possess a more isotropic (random) alignment. This difference in alignment causes these collagen types to exhibit significantly different mechanical properties and behavior. In order to implement the distinction between the two types of collagen, more quantitative information should be available regarding the production and mechanical properties of these collagen types. In general, model extensions such that the underlying biology and physics are described in more detail lead to the need for additional input parameters, which may have yet to be measured or even be hard to measure. In addition, these extensions often lead to additional uncertainty; hence, such extensions only sometimes improve the quality of the model. However, adding the distinction between the two collagen types will lead to a better description of the development of the contraction, particularly the development of hypertrophic scars.
Given that the finite element simulations take much time and many Monte Carlobased predictions are needed, for example, to reach personalized healthcare, we prefer the application of neural networks. We can train a neural network on the data we generate with the finite element method and then make predictions about contraction. Training a neural network takes time; however, it gives predictions at ultra-high computational speed once trained. This quick speed will make running an application on a smartphone or tablet possible.
There is much variation between burns and patients, which entails uncertainty. These uncertainties are twofold. On the one hand, some parameters are not or hardly measurable. For this purpose, we must apply inverse modeling to estimate parameter values and link the results to (clinical) experiments. On the other hand, parameter values are patient-dependent, so we need to perform uncertainty quantification. Therefore, inter-parameter dependency and patient-specific factors need to be investigated.
For the sake of parameter sensitivity analysis, we need to perform Monte Carlo simulations, and for higher-dimensional frameworks, this can become a challenge because of computational complexity. Hence, we require implementing the finite element method in a high-level programming language like C++. We can also consider clever Monte Carlo techniques based on many simulations with low numerical resolution and a few with high numerical resolution. For these Monte Carlo simulations, it is essential to consider both the model's sensitivity and stability. However, we need highquality finite element simulations to train a neural network. Therefore, for Monte Carlo simulations only, we can use the method of many low numerical resolution simulations and a few simulations with high numerical resolution. Training a neural network in higher-dimensional frameworks requires only high numerical-resolution simulations.
Besides training a neural network on the generated data from simulations, it is also an option to train directly on clinical data. We can think of datasets that provide burn images as input, supplemented with the patient's data, the burn's location and cause, and the course of contraction, contracture, and hypertrophy. Creating such datasets requires a systematic follow-up or a hybrid approach that works with variable data when the data is limited.
For such a neural network, we can fit the wound using a convolutional neural network that takes images of the initial wound. For these images' contours and other features, we can use pixel-based methods, shape similarity (Andreou and Sgouros 2005), and shape matching (Veltkamp 2001).
This study presented a stability analysis for the fully continuous and semi-discrete version of the two-dimensional model for post-burn contraction. We could analytically determine the eigenvalues, which is possible because the linearised equations (13)-(15) leave out other variables after accounting for the equilibria values. As a result, some eigenvalues meet the stability constraints independent of the chosen value for the parameters, given that the parameters involved are positive and realistic. We have illustrated that for the parameter range that we are interested in, there is a region around the origin of the (ε 11 , ε 22 , ε 12 )-space where the model is stable with respect to perturbations with all wave frequencies. This condition is sufficient, meaning the current analysis does not exclude any other stable equilibria. Nonlinear effects will, at most, be able to induce constant-state instabilities for the effective strain. Further, an important stability constraint states that the model is stable for low signaling molecule decay rates, though not too low.
We have shown consistency between the semi-discrete model's eigenvalues and the continuous model's eigenvalues. If the equilibrium solution to the continuous problem is stable, then the equilibrium to the semi-discrete problem is stable under the present discretization. The obtained eigenvalues of the system establish the convergence rate towards the equilibrium. We have assessed the convergence of the numerical method experimentally, in which the order of convergence is above O(h 2 ). Since the difference between the chemical eigenvalues from the continuous and semi-discrete problem is of the order O(h 2 ), the convergence rates towards the equilibrium differ by an order O(h 2 ). This result is better than expected since the discretization method should have local truncation errors of order O(h 2 ).
Using numerical simulations, we validated the stability constraints derived from the analysis. If the input values satisfy the stability criterion, the model behaves stable given these stable parameter values.
The model can numerically be unstable if the parameters do not meet the signaling molecule stability constraint. If δ c < k c /(a I I c ρ) is not too far below the bound, then, initially, the model seems stable and the healing proper. The displacement velocity perturbations vanish quickly, and the signaling molecule perturbation shifts such that the density peak moves to the center of the wound in the first 13 days. One would expect this density to decrease from that day on, given that the other chemicals seem to reach equilibria. However, shortly after 13 days, the signaling molecule density moves away from equilibrium, affecting all the variables except displacement velocity. The distributions and densities move away from the expected equilibria and oscillate around new equilibria, where the densities remain. We have linked this situation to real-life occurrences of hypertrophic scars and keloids. By reverting the stimulation of matrix production and myofibroblast differentiation by setting the signaling molecule density to (healthy) equilibrium and removing the myofibroblasts, we have provided experimental evidence, from a mathematical point of view, that one can indeed restore the fibroblast cell and collagen density to healthy equilibria.
If the signaling molecule decay rate is too far below the stability limit, then the numerical method does not converge and loops over Picard iterations while remeshing. Taken together, the numerical model fully reproduces the stability constraints.