Comparison between a phenomenological approach and a morphoelasticity approach regarding the displacement of extracellular matrix

Plastic (permanent) deformations were earlier, modeled by a phenomenological model in Peng and Vermolen (Biomech Model Mechanobiol 19(6):2525–2551, 2020). In this manusctipt, we consider a more physics-based formulation that is based on morphoelasticity. We firstly introduce the morphoelasticity approach and investigate the impact of various input variables on the output parameters by sensitivity analysis. A comparison of both model formulations shows that both models give similar computational results. Furthermore, we carry out Monte Carlo simulations of the skin contraction model containing the morphoelasticity approach. Most statistical correlations from the two models are similar, however, the impact of the collagen density on the severeness of contraction is larger for the morphoelasticity model than for the phenomenological model. Supplementary Information The online version contains supplementary material available at 10.1007/s10237-022-01568-3.


Introduction
Generally speaking, wound healing contains four overlapping phases for the skin to cure itself after getting injured. It is a complicated process with many biological elements involved, for instance, cells, cytokines and tissue molecules. They collaborate in the form of biological events, which contribute to the resurfacing, reconstitution and restoration of the tensile strength of the injured skin. Superficial skin injuries, where only the epidermis is damaged, heal without any problem. However, deeper wounds, where the damage extends into the dermal layer, cause more complications. The severeness depends on the size and depth of the wound.
In many cases, deep tissue injury causes contraction of the damaged skin tissue. Severe contractions that cause loss of mobility of joints are commonly referred to as contractures. Hence, contractures usually occur with disabilities of the joints of the patient. Contractures take place due to the large traction forces exerted by the (myo)fibroblasts on the extracellular matrix. Skin contraction occurs in the proliferation phase, which continues for two to four weeks post wounding. The contraction process is characterized by active proliferation and differentiation of fibroblasts. The main healing mechanisms in this phase are re-epithelialization, fibroplasia, angiogenesis and the development of granulation tissue. During this stage, fibroblasts are attracted to the damaged region by a number of chemokines like platelet-derived growth factor (PDGF) and transforming growth factor-beta (TGF-beta). Induced by the high concentration of TGF-beta, some fibroblasts differentiate into myofibroblasts, which is a cell phenotype having the properties of both fibroblasts and smooth muscle cells (Darby et al. 2014;Grinnell 1994;Phan 2008;Chaudhari 2015). Both fibroblasts and myofibroblasts are responsible cell phenotypes for the regeneration of collagen bundles and the reestablishment of the extracellular matrix (ECM). Compared to fibroblasts, the myofibroblasts exert even larger forces on the ECM and secret more collagen bundles, and as a result of myofibroblasts, the damaged skin contracts more. For acute wounds, myofibroblasts are 1 3 often found to persist, which are otherwise supposed to be extinct as a result of apoptosis (programmed death) (Frangogiannis 2017; Xue and Jackson 2015). The amount of contraction is related to the size, shape, depth and anatomical location of the wound. For instance, tissues with stronger laxity contract more than loose tissues, furthermore squareshaped wounds contract more than circular ones (Enoch and Leaper 2008;Li and Wang 2011). Usually, a reduction by 5 − 10% of the original wound volume is observed in clinical studies. For a more detailed biological description on wound healing, we refer to Enoch and Leaper (2008); Li and Wang (2011).
Many different models have been developed for all the phases of wound healing, see for instance Boon et al. (2016); Koppenol (2017); Cumming et al. (2009a); Murray (2003); Ziraldo et al. (2013); Cohen and Mast (1990). These models can be categorized as agent-based models and continuum-based models. Since we are interested in the cause of contractions and contractures, an agent-based model is preferred, since all the cellular events and cell positions can be tracked explicitly. To model the traction forces exerted by the (myo)fibroblasts, we divide the cell boundary into line segments, then use a point force at the midpoint of every boundary segment. Regarding the deformation of the ECM, there are two different approaches, namely a phenomenological approach and a morphoelasticity approach.
In Peng and Vermolen (2020), we modelled the traction forces exerted by the (myo)fibroblasts with a phenomenological approach, in which the virtual forces are applied at the line segments of the mesh elements. An alternative approach to model plastic deformation is by the use of morphoelasticity, which is a more advanced model. Morphoelasticity is widely used in biological modelling to describe elastic growth, for instance, the growth of tumors (Goriely and Moulton 2011), the growth of seashells (Rudraraju et al. 2019) and contractions of scars after an injury (Koppenol 2017;Ben Amar et al. 2015) etc. Morphoelasticity provides a description of the evolution of plastic (permanent) deformation.
This manuscript is structured as follows: we firstly introduce both approaches in Sect. 2. Section 4 describes the investigation of the impact of various parameters on the morphoelasticity approach. Subsequently, we show some numerical results in two dimensions to compare these two approaches in Sect. 5. In Sect. 6, we perform a Monte Carlo based parameter sensitivity analysis on the basis of the morphoelastic approach for the displacement of the extracellular matrix (ECM). Finally, Sect. 7 presents the conclusions.

Mathematical models
In this section, mainly the force balance part of the skin contraction model will be discussed. For the sake of completeness of this manuscript, we present a brief summary of the skin contraction model developed in Peng and Vermolen (2020), which contains all the other mechanics like chemotaxis and interaction between cells etc.
To model the force exerted by a (myo)fibroblast, we divide the cell boundary into line segments, then a point force is applied at the midpoint of each line segment, which points inward to the center of the cell; see Figure 1 where we use a square to approximate the cell and the point forces are applied on each line segments. The point force is decribed by the Dirac delta distribution, which is defined in any dimensionality by and constrained to satisfy the identity that Here, Ω is an open subset of ℝ n . Therefore, we use the expression in Peng and Vermolen (2020) to describe the total forces exerted actively by all the cells in the computational domain: where T N (t) is the number of cells at time t, N i S is the number of line segments of cell i, P(x, t) is the magnitude of the pulling force exerted at point x and time t per length, n(x) is the unit inward pointing normal vector (towards the cell centre) at position x , x i j (t) is the midpoint on line segment j of cell i at time t and ΔΓ i,j N is the length of line segment j. Note that the above equation represents the cellular forces in case of a (1) Fig. 1 A schematic of the cellular traction forces exerted by a cell, which is approximated by a square general polygonal (in ℝ 2 ) or a general polyhedral (in ℝ 3 ) division of the cell boundary. This method is also used in fluid dynamics and is known as the immersed boundary method. Theoretically, as N i S → ∞ , i.e. ΔΓ i,j N → 0 , Eq (1) becomes (Boon et al. 2016) Here, x i (t) is a point on the cell boundary of cell i at time t. In the current manuscript, we do not consider any shape changes of cells, and therefore, we divide the boundary of each cell in four segments on which a force towards the center is exerted. This has been visualized in Figure 1. This approach has been used in both the earlier modelling framework in Peng and Vermolen (2020) and the current modelling framework.

Phenomenological approach
In Peng and Vermolen (2020), a phenomenological model was developed to simulate permanent deformations of skin. The idea was based on the inclusion of artificial forces on the faces (that is, boundary segments) of the mesh elements. This approach is motivated by the fact that it is known that myofibroblasts shorten the polymeric chains that constitute the extracellular matrix. The phenomenological forces on the boundary of mesh elements are also referred to as 'virtual forces' or 'plastic forces'. The latter qualification is motivated by the permanent nature of these forces. In Peng and Vermolen (2020); Koppenol (2017), the magnitude of the virtual forces in a mesh element is determined by the duration that myofibroblasts occupy the region of the mesh element. The principle behind the determination is based on solving a ordinary differential equation that leads to a chemical equilibrium that reflects the shortening of the polymeric chains. In this manuscript, our focus lies on the qualitative comparison between two approaches, hence, we simplify the phenomenological approach here by assuming the magnitude of the virtual force to be constant and to be equal to Q. For the plastic forces in the phenomenological model, we set where N E is the total number of mesh triangular elements, N i e is the total number of edges of mesh elements (hence for triangular elements, N i e = 3 ), Q is the magnitude of the plastic force density of myofibroblast, n(x) is the unit inward (hence directed towards the center of the mesh element) pointing normal vector, x i e (t) is the midpoint of the boundary segment in consideration and ΔΓ i,e E is the length of the boundary segment, respectively. The virtual force ( f p (x;t) ) has the similar form as the expression of the total force exerted by the cell (see Eq (1)), however, the virtual force is applied on the mesh line segments to describe the permanent force, which simulates the residual forces even after the cell dies and has disappeared from the computational domain. Therefore, this virtual force only occurs in the phenomenological model but not in the morphoelastic model. We note that this virtual force framework is sufficiently general to extend to different element shapes. Furthermore, we assume that the displacement of the tissue vanishes far away and hence a homogeneous Dirichlet boundary condition is used. Hence, the boundary value problem is expressed by Here, contains both elastic and viscoelastic part, that is, Furthermore, is the infinitesimal strain tensor as defined by and tr ( ) is the trace of the tensor, and f t and f p are defined by Eq (1) and Eq (3), respectively. Furthermore, ̇ and u represent the time derivative of and u , respectively.
To model a wound, a subdomain Ω w ⊂ Ω (strict subset) is chosen. Further, zero initial displacement is assumed, that is u(x, t) = 0.

Morphoelasticity approach
Whereas the phenomenological approach is based on directly solving the equations for the local displacement of the tissue, the morphoelastic approach represents partial differential equations for the displacement velocity and the effective strain . Note that the effective strain does not necessarily satisfy Eq (6). To this extent, the equations for the balance of momentum and the evolution of the strain tensor are given by: where is the density of the extracellular matrix, f t is the temporary force described in Eq (1), L = ∇v and is a nonnegative constant. Here, Dy Dt = y t + v∇y is the material derivative for any tensor field y and v is the displacement velocity. Furthermore, skw (L) = 1 2 (L − L T ) is the skewsymmetric and sym (L) = 1 2 (L + L T ) is the symmetric part of the tensor, respectively, and tr ( ) represents the trace of . The definition of the stress tensor in the morphoelasticity approach has the same components as in the phenomenological approach in Eq (5). Since one solves for the displacement velocity and the effective strain in the morphoelastic approach, the viscous part in Eq (5) is replaced with The initial boundary value problem for morphoelasticity is nonlinear and requires the solution of the displacement velocity v and the strain tensor . This is in contrast to the boundary value problem in Eq (4). Since we are interested in the contraction of the scar region Ω w , we approximate the displacement by integrating the velocity over time:

A summary of the skin contraction model
For the sake of the completeness of this paper, we briefly present the model developed in Peng and Vermolen (2020). The model contains various biological segments, for instance, signalling molecules, cells and tissue bundles. For a more detailed description of the skin contraction model, we refer to our previous work (Peng and Vermolen 2020).

Signalling molecules
There are three categories of signalling molecules in the model. The concentration of each signalling molecules is expressed by the convection-diffusion equation: with Robin's boundary condition: where c is the concentration of the signalling molecule, D is the diffusion rate, v(t, x(t)) is the displacement velocity of the ECM (solved from one of the force balance models in Section 2), F is the reaction term depending on the category of the signalling molecules, and c is a positive constant.
To distinguish between the cytokines, we add the subscript "PDGF", "TGF" or "tPA" if it is necessary (i.e. c PDGF , c TGF and c tPA ). Speaking for the initial condition, PDGF has higher concentration in the wound region and lower in the undamaged region; there is no TGF-beta over the domain initially; the concentration of tPA is higher over the edge between the wound and undamaged skin and lower in the rest of the computational domain. Dallon et al. (1999) and Cumming et al. (2009b) developed a tensor-based representations of tissue bundles, such that the density and the orientation of the tissue bundles can be described. The tensor of the orientation of the tissue bundle is given by where k indicates the type of the tissue: k = f for fibrin and k = c for collagen, p( ) T = [cos , sin ] is the unit vector in the direction , and (x, t, ) is the density of collagen at position x and time t. Since the tensor is symmetric positive definite, it can be decomposed as the sum of weighed outer products of orthogonal eigenvectors, which in the twodimensional case gives:

Tissue bundles
where 1 (x, t) and 2 (x, t) are eigenvalues, and v 1 (x, t) and v 2 (x, t) are corresponding eigenvectors. We assume that the fibrin bundles are degradated by the concentration of tPA, hence, for every entry in f , it is determined by for any i, j ∈ {1, 2} and represents the degradation rate of fibrin bundles. The collagen bundles are deposited by the (myo)fibroblasts in the direction of active migration. Since the cells are much smaller than the computational domain, we use Dirac Delta functions to model the secretion of the collagen bundles. Furthermore, the secretion rate depends on the amount of total density including fibrin and collagen: vdt‖ is the unit vector of the direction of active displacement of (myo)fibroblast k at time t for any i, j ∈ {1, 2} , and x k N (t) is the centre position of (myo)fibroblast with index k. Initially, it is assumed that the collagen and fibrin are isotropic. Therefore, we have and where Ω w is the indicator function that was defined as where Ω w is the wound region as a subdomain in the computational domain. and is a positive constant.

Cellular activities and displacement of cells
In this model, we assume that only regular fibroblasts can proliferate and differentiate to myofibroblasts, while macrophages and myofibroblasts will only be subject to apoptosis (programmed death). These cellular events are modelled by stochastic process, which follow the exponential statistical distribution.
Cells migrate due to various mechanisms, for instance, chemotaxis, random walk, passive convection and interaction between cells. The collagen bundles also have an immediate impact on the migration of (myo)fibroblasts. For the immune cells, chemotaxis is associated with PDGF, then the displacement of the center position of i-th macrophage is given by where c is the constant representing the weight of chemotaxis, which is expressed by here, v is the speed of biased movement of cells, c PDGF is the concentration of PDGF which is initially high in injured region and low in uninjured region, is a small positive constant to prevent the denominator from being zero and v is the displacement velocity of the substrate, which follows from solving the momentum balance displacement.
For the (myo)fibroblasts, similar to the displacement of microphages, the new positionof the i-th (myo)fibroblast is updated via where dr i (t) is calculated from where c is the same expression as before in Eq (11), and c TGF is the concentration of TGF-beta secreted by macrophages.
For a detailed description of the model and the explanation of all the parameters, we refer to our previous work Peng and Vermolen (2020).

Comparison between two approaches
Both the phenomenological and the morphoelastic approaches can be used to simulate the permanent contraction of post-wounded skin. The two models have the agentbased character for the cell dynamics in common. The cellular dynamics regarding migration, division, differentiation and the proliferation of collagen, including its orientation are identical. The momentum balance differs between the two models in the sense that the phenomenological model explicitly computes the displacement at all positions in the computational domain. Furthermore, the permanent displacements are modelled by the use of 'apparent', or 'virtual' forces that arise depending on the exposure of an mesh element to myofibroblasts. The physics behind this assumption is motivated by the assumption that the myofibroblasts induce shortening of the polymeric chains, which leaves residual stresses around these fibers. The morphoelastic model is based on exposure to strain. A large exposure to strain, induces significant built-up of permanent deformation. The formalism that we used here is inspired from the model in Koppenol (2017), where morphoelasticity is considered in the context of post-trauma skin contraction. In that model, the strain and the number of Matrix metalloproteinases (MMPs) and myofibroblasts determine the extent of permanent contraction. The current modelling framework does not account for MMPs and hence uses a somewhat simplified relation between permanent displacement and local strain.
The similarities and differences between the morphoelastic model and the phenomenological model are summarized in Table 1. Firstly, the unknown parameters to be solved in these two models are different: one can obtain the displacement of the ECM directly by solving the linear system from the phenomenological model, however, one needs to approximate the displacement in the morphoelastic model, which is a nonlinear system solving both the velocity and the strain tensor of the domain. Secondly, the phenomenological model does not consider inertia, which the morphoelastic model does. Thirdly, the permanent displacements in the phenomenological model are solely determined by the history path of local presence of myofibroblasts (in a mesh element), and this is modelled in the balance of momentum directly, whereas the permanent displacements in the morphoelastic model follow from the history path of the local strain as a result of the forces exerted by both the myofibroblasts and the fibroblasts. In other words, the exposure to mechanical triggers invokes the permanent deformations in the morphoelastic formalism, whereas the actual presence of myofibroblasts triggers the permanent deformations in the phenomenological model. Speaking of the similarities, both models in this manuscript are agent-based models, therefore, they are mainly suitable for the microscale due to the high computational cost. Furthermore, both of them contain various stochastic processes regarding random walk, cell activities (i.e. proliferation, apoptosis and differentiation).

Sensitivity test of morphoelasticity approach
As the model mainly deals with the deformation of the wound, we define a subdomain (as scar region) in the centre of the computational domain (see Fig. 2) and we investigate the relative scar area at different times. Here, we define the relative scar (or wound) area by where Ω w is the scar (or wound) region, A(Ω w ;t) is the area of Ω w at time t. Since we are interested in the impact of the input parameters on the contraction of the skin tissue, we carry out sensitivity tests where we only change the value of one parameter. The parameter values are displayed in Table 2.
In this section, we assume that the shapes and sizes of the cells are identical, and that they are at fixed locations that are strictly within the wound region Ω w ⊂ Ω . The traction forces are released at all times from t = 5 day. For reasons of comparison and presentation, we define t min as the time point at which the scar exhibits the maximal contraction, defined by A Ω w (0) , Table 1 Similarities and differences between the phenomenological model and the morphoelastic model , that is, the area of the scar is minimal here: Note that t min > 0 . We are further interested in how the scar flows back towards its (new) equilibrium. Since this adjusted equilibrium is an asymptotic value, we define t end by In this manuscript, we use a finite-element method to solve all the boundary value problems. Regarding the timeintegration, we use a backward Euler method. From the theory, it is known that smooth solutions would be subject t end ∶= arg min t>t min to errors of the order O(h 2 ) in the L 2 norm of the numerical approximation and O(Δt).
The numerical solutions of most boundary value problems are naturally smooth and convergent, except for the force balance model. The numerical stability and convergence of morphoelasticity model has been discussed in our other work Peng, Vermolen and Weihs (2021): with the refinement of the mesh, the convergence rate of the L 2 − norm of the solution to Eq (7) (i.e. the velocity) is computed, which is 1.899828112 that is close to the theoretical value 2.
The Young's modulus E s denotes the stiffness of the ECM. According to clinical observations (O'Leary et al. 2008), softer skins will develop a more serious skin contraction compared to stiffer skin. In other words, since the amount of contraction and displacement decreases with increasing stiffness ( E s ), the scar area stays larger if the stiffness is larger. We show the results in Fig. 3  values of E s . In Fig. 3a, the relative scar area is plotted as a function of E s and various colors of the curves stand for different time points. Figure 3b shows the relative scar area as a function of time and various colors of the curves stand for different values of E s . It is concluded that the model confirms that a stiffer skin develops less contraction. Furthermore, Fig. 3a shows a hyperbolic correlation between the wound area and the stiffness. Skin is taken as a viscoelastic material (Kuwahara et al. 2008), that is, skin has both viscous and elastic properties. Viscosity is related to delayed recovery (or damping) from deformation, while elasticity is related to rebounding and quick recovery from deformation. Hence, we are interested in the impact of the damping property of the skin on the wound contraction. As it is shown in Fig. 4, with a larger viscosity , the wound reaches its minimal area later and develops less contraction. In particular, there is an a weakly nonlinear correlation between the minimal wound area and the value of . Similarly, for the recovering phase, a smaller viscosity results into faster recovery, that is, the model reaches the equilibrium earlier, i.e. at a smaller t end . It can be seen from Fig. 5 that scar contraction takes place over a longer period with a linear tendency as increases. However, less significant differences appear in the final stages of scar contraction, and this tendency may be attributed to the fact that in the final stages relaxation takes place, where the equilibrium does not depend on the viscous part of the model.

using various
The Poisson ratio indicates how much a material will deform in the direction perpendicular to the force direction. We vary here between (0.118, 0.495). Note that for incompressible materials, we have = 0.5 . As the ECM becomes more incompressible, the wound area stays constant and this yields small contractions regardless of the time point post wounding. Furthermore, the curves in Fig. 6a show a hyperbolic behaviour. In general, a more compressible ECM indicates a higher degree of skin contraction and a larger healing time; see Fig. 6b.
Different from the phenomenological approach, the degree of the permanent deformation in the morphoelasticity approach is determined by the value of in Eq (7). Figure 7 displays simulation results with different values of . Generally speaking, a smaller value of is beneficial to the scar since this leads to smaller final contractions. Note that when = 0 , theoretically the wound will recover to its original size. However, due to numerical errors, in Fig. 7b, the final area ratio is slightly above 1. In fact the case = 0 reflects the viscoelastic case with the addition of inertia.   (14) when the scar reaches its equilibrium

Comparison between the phenomenological approach and the morphoelasticity approach
As both approaches are capable of reproducing the permanent plastic deformation during skin contraction, we are interested in whether they are numerically consistent, when both boundary value problems are solved by the finite-element methods with Lagrangian piecewise linear basis functions. We altered the skin contraction model in Peng and Vermolen (2020) by replacing the phenomenological model for plastic deformation (see Eq (4)) with the morphoelastic formulation (see Eq (7)). Most of the parameter values in the current study have been adopted from Peng and Vermolen (2020). These parameter values involve the properties of the cells and mechanical parameters of the ECM like the elasticity, viscosity and Poisson ratio. However, the current (new) model differs in the use of morphoelasticity. The parameters that were needed for this inclusion have not been taken from our earlier work. Hence, we only display the parameter values used in the force balance of the skin contraction model in Table 3. We have attempted to use realistic values for the input variables. However, some values of input variables are subject to large variations. The reasons for this variability are the patient-specific nature and poor availability of measurements that support these values. A typical example is the stiffness (Young's modulus) of skin, which varies very much with the location in the body and varies a lot from patient to patient. Variations may even exceed a factor of 1000 (Liang et al. 2010; Trotta and Annaidh 2019; Wang and Lakes 2002). As a follow-up of our previous work in Peng and Vermolen (2020), we inherit the most settings and the initial domain is shown in Fig. 8. The results of the relative wound area against time are shown in Fig. 9. We run several simulations with each approach since the model contains random processes. Figure 9 presents a case without and with permanent deformations (left and right subfigures, respectively). Figure 9a shows the case when plastic forces are not active in the phenomenological approach and the case that = 0 in the morphoelastic approach. In other words, after the forces being removed from the computational domain, the wound will recover to its original volume, that is, no contraction occurs if f p = 0 in the phenomenological approach and = 0 in the morphoelastic approach. Figure 9b shows the results using the corresponding parameter values from Tables 2 and 3. In Fig. 9a, the curves of both colors are mostly overlapping, which indicates that both approaches can be used to describe scar contraction and as long as the main parameter values are the same, the models are more or less consistent. The inclusion of permanent deformations in both models by setting and Q to positive values shows the same qualitative trends. However, it has turned out to be very difficult to make the curves from the two different models overlap, that is, make them agree quantitatively. We have carried out multiple attempts to reduce the quantitative difference between the results of both modelling frameworks. Either the minimal scar area from the two models agreed with other, or the final scar area did. We think that from a qualitative point of view the models predict the same kind of behavior. However, this difference shows that the two modelling frameworks are not the same. In case of the morphoelastic formulation, which was also used in the work by Koppenol (2017), the permanent deformation is solely modelled as a consequence of the strain (history path) that the tissue endures. The history path of the strain is caused by the forces that were exerted by both the fibroblasts and myofibroblasts. In case of the phenomenological approach, the permanent deformation is modelled by the actual presence of the myofibroblasts only, and not by the strain history path. This gives a different model formulation, though we observe similar evolution of skin. A next step would be to compare the outcomes from both modelling frameworks to small scale (in vitro) observations. Another difference between both modelling frameworks is in the balance of momentum. The balance  of momentum in the morphoelastic approach involves inertia, which is more interesting from a numerical point of view, and also closer to physics, whereas inertia has not been accounted for in the phenomenological approach. The difference in the treatment of the balance of momentum does not give rise to significant differences between the modelling outcomes, see Fig. 9b. We finally remark that the differences between both modelling frameworks could possibly be reduced if the morphoelastic would only depend on the local presence of myofibroblasts. This was beyond the scope of our study. Regarding computational efficiency, these two approaches take more or less the same CPU time for each iteration (around 4 − 7 s an Intel(R) Core(TM) i7-6500U CPU @ 2.50GHz computer). Hence in terms of computation time, there is no preference for either two models.

Monte Carlo simulations with morphoelasticity
In the current simulations, we use the cell-based model as described in Peng and Vermolen (2020). Cells, treated as separate entities, are subject to migration (chemotaxis and random walk), proliferation, differentiation and apoptosis. Random walk is modeled by means of a vector Wiener process, and chemotaxis is modeled by an additional term migratory term that is proportional to the gradient of the TGF-beta. Furthermore, the random cellular processes such as apoptosis, cell differentiation and cell proliferation are modeled by the use of sampling from exponential distributions. More information can be found in Peng and Vermolen (2020). Therefore, the complete model contains a large amount of randomness due to the random processes like cell division and migration, and due to the uncertainty in the values of the large number of input parameters. The input parameters are descriptive for the cells, signalling molecules and tissue bundles. The outputs of interest include the scar area and collagen density at different times. We are, for instance, interested in the time at which the scar area reaches its (modified) equilibrium (i.e. when the scar area does not change significantly anymore; see Eq (14)), the time at which the contraction is maximal and the time at which the total amount of collagen reaches a certain threshold value. Note that, in this manuscript, the size of the cells and the scar have been set to academic values due to the limitation of computational resources. The values, as well as the statistical distributions of the input values are listed in Table 4. Other input values have been taken from Peng and Vermolen (2020). In total, 1005 samples are collected and the basic statistics, being the obtained output sample mean and output sample standard deviation, are shown in Table 4. With respect to the input variables, the stiffness has the most significant impact on the outputs, which is reasonable since stiffness directly relates the magnitude of displacement to stresses and forces. Regarding the outputs, they can be categorized as wound area, specific time point and collagen density ratio, which is defined as follows: where c is the density of the collagen and 0 c is the density of collagen in the undamaged skin. This indicator illustrates the difference between the injured, which is gradually recovering, and uninjured skin.
In Table 4, an estimate of the Monte Carlo error, defined by the standard deviation of the mean, is shown. As the number of samples increases, one can approximate the relative Monte Carlo error bŷ Fig. 9 The plots show the wound area ratio with both phenomenological and morphoelasticity approach with and without permanent deformation, respectively. Red curves represent the morphoelasticity approach and blue curves represent the phenomenological approach.
The parameter values are given in Tables 2 and 3  Table 4 Distributions of the input parameters in Monte Carlo simulations when morphoelasticity approach is collaborated where ̂ and ̂ are the sample standard deviation and sample mean, respectively, and n is the sample size. By the Central Limit Theorem (Klenke 2008), one can obtain immediately that the quantity ̂− ∕ √ n follows the standard normal distribution, where and 2 are the population mean and variance, respectively. Given , > 0 , we suppose then We note that (n − 1)̂2 2 follows a 2 distribution with n − 1 degrees of freedom, and therewith √ n̂− follows a t-distribution with n − 1 degrees of freedom. However, as n → ∞ , ̂ converges to (consistent and unbiased estimator). Hence, for n very large, it follows that where Φ(z) is the cumulative probability function of standard normal distribution. We write the above inequality in terms of the relative Monte Carlo error (denoted by r.e.), then one obtains The above inequality in terms of an interval of significance says that the sample mean resides within this interval with estimated probability of (1 − ) . From Table 4, all the relative Monte Carlo errors are less than 0.006, which indicates that regarding the sample mean, the sample size is sufficiently large. Compared with the results of Monte Carlo simulations in Peng and Vermolen (2020) (see Table 5 in Appendix) where the phenomenological model is used, for the equilibrium and minimal scar area, the standard deviations from this manuscriptis are both one order of magnitude larger. This shows that the morphoelastic model is more sensitive with respect to the changes of input parameter values.
Most of the correlations obtained in this manuscript for morphoelasticity are consistent with the correlations that were obtained in Peng and Vermolen (2020). The current correlations seem to be a bit more significant, which indicates that the model based on morphoelasticity is more sensitive to parameter variation than the phenomenological model from Peng and Vermolen (2020), whereas the standard deviations of most responses are comparable. In other words, the use of morphoelasticity does not significantly reduce the uncertainty in the model. At least significant ( p < 0.001 ) correlations are observed between many of the output and input parameters. Therefore, the scar area after a short amount of time after wounding can be used as indicative for the behavior of the scar area at later times, such as the maximum contraction and the final contraction. Figure 10 shows some scatter plots between various parameters. From Fig. 10a, it can be seen that the correlation coefficient is (almost) equal to one, which reflects a very strong correlation between the maximal contraction and the final contraction. In other words, the more the scar contracts, the more the long-term area differs from the initial area. In other words, there is less absolute recovery in the long term if the earliest maximum contraction is more severe. Furthermore, there is a slightly stronger correlation between the final scar area and the scar area on the fourth day post wounding (see Fig. 10b), compared to the results in Peng and Vermolen (2020), which were based on the phenomenological model for morphoelasticity. This slightly stronger correlation may be a consequence from the fact that the morphoelasticity approach directly connects the scar deformation to the strain tensor. These correlations are useful for clinicians to predict the final contraction of the injured skin.
The stiffness of skin significantly determines the extent of contraction. It has been found that a wound in the buttock contracts more than a wound in the scalp region (O'Leary et al. 2008). In Fig. 10c, it can be seen that the model reproduces this biological observation, for example, it is observed in O' Leary et al. (2008) that the scar contracts less in the scalp (the stiffness is in the order of 20 − 40MPa (Trotta and Annaidh 2019)) than in the buttock (the stiffness is in the order of 1MPa (Wang and Lakes 2002)). Subsequently, more wound contraction goes with higher collagen density ratio. Therefore, there is a negative correlation between the collagen density ratio and the skin stiffness; see Fig. 10d.
Another difference of the results from the phenomenological model in Peng and Vermolen (2020) is the correlation between the final wound area and the final ratio of average collagen density. In Peng and Vermolen (2020), even though the correlation is significant (i.e. p value is less than 0.001), the scatter plot is not monotonic, therefore, the negative correlation is not convincing to some extent. However, in the current morphoelasticity study, we observe, among other correlations, a significant correlation between the scar area and the ratio of density of collagen. In Fig. 10e, the final wound area and the final ratio of average collagen density are strongly negatively correlated. In other words, higher collagen density results in a higher degree of contractions (i.e. a smaller final scar area). This also holds for the data from the 4th day post wounding; see Fig. 10f. This dependence is caused by the fact that the orientation of the collagen bundles guides the migration direction of the (myo) fibroblasts, and since the (myo)fibroblasts are responsible for the regeneration of collagen bundles, as well as for the contractile forces within the scar, the combination of proliferation of (myo)fibroblasts and remodeling of collagen lead to more contraction. This is also observed in Thomas (2001).
To determine the degree of contraction of post-burned skin, the scar area is used as a quantifier. Figure 11 shows the estimated probability density distribution of the final wound area from the Monte Carlo simulations. It is obvious from the shape of this probability density distribution that the final scar area is not normally distributed. A normality test confirms this observation quantitatively. From  Fig. 11b, it can be concluded, that considering a wound with size 20 m × 15 m , the probability that the wound will develop 10% contraction is around 40% and 15% contraction is around 10%.

Conclusions
In this manuscript, we consider a formalism to describe the permanent deformation of a scar region by combining morphoelasticity and a cell-based wound healing model. We confirm that both the phenomenological approach and the morphoelasticity approach qualitatively give similar numerical results. Hence, both model formulations can be used to describe to local plasticity of the wound. The parameter sensitivity analysis shows that the viscosity part in the viscoelastic model for the displacement has a significant damping effect on the evolution of the scar area.
Compared to the phenomenological approach, the model with morphoelasticity is more sensitive to the random processes and parameter variations. Furthermore, with the same magnitude of the cellular forces, the model with morphoelasticity predicts larger contractions. Regarding computation time, the solution of the partial differential equations is more expensive in the case of morphoelasticity than for the phenomenological model. However, for large numbers of cells, this difference decreases, which implies that the computation time becomes predominantly determined by the cellular processes. These cellular processes entail the cellular force interaction, development of cellular forces, cell proliferation and cell migration. Based on these considerations, the incorporation of morphoelasticity into the modeling framework is also promising in real-world situations where a parallel computing environment is used.
Next to this deterministic parameter variation, a Bayesian (Monte Carlo) approach has been done for the morphoelastic model. Most of the results confirm the correlations that were obtained from the Monte Carlo simulations on the basis of the phenomenological model. However, the correlation between the ratio of average density of collagen and the scar area is more negatively pronounced in the morphoelastic model. This observation is consistent with the clinical observations in O' Rourke et al. (2015) and in Thomas (2001). For this reason and because of the physics-based nature of morphoelasticity, the model based on morphoelasticity is preferred for use in future research.
Regarding further research, the focus should be directed towards application of the model to an enlargement of the scale and extension to three spatial dimensions, so that the modeling framework is applicable to clinical situations. This enlargement of the modeling scale requires a parallel computing environment. Furthermore, upscaling of the agent-based model to a fully continuum-scale computational framework is also a sensible step.

Appendix: Results of Monte Carlo simulations with the phenomenological model
To compare with the results of Monte Carlo simulations with the morphoelastic model shown in Sect. 6 (see Table 4), we exhibit the results from Peng and Vermolen (2020) here: