Sensitivity and feasibility of a one-dimensional morphoelastic model for post-burn contraction.

We consider a one-dimensional morphoelastic model describing post-burn scar contraction. Contraction can lead to a limited range of motion (contracture). Reported prevalence of burn scar contractures are 58.6% at 3-6 weeks and 20.9% at 12 months post-reconstructive surgery after burns. This model describes the displacement of the dermal layer of the skin and the development of the effective Eulerian strain in the tissue. Besides these components, the model also contains components that play a major role in the skin repair after trauma. These components are signaling molecules, fibroblasts, myofibroblasts, and collagen. We perform a sensitivity analysis for many parameters of the model and use the results for a feasibility study. In this study, we test whether the model is suitable for predicting the extent of contraction in different age groups. To this end, we conduct an extensive literature review to find parameter values. From the sensitivity analysis, we conclude that the most sensitive parameters are the equilibrium collagen concentration in the dermal layer, the apoptosis rate of fibroblasts and myofibroblasts, and the secretion rate of signaling molecules. Further, although we can use the model to simulate significant distinct contraction densities in different age groups, our results differ from what is seen in the clinic. This particularly concerns children and elderly patients. In children we see more intense contractures if the burn injury occurs near a joint, because the growth induces extra forces on the tissue. Elderly patients seem to suffer less from contractures, possibly because of excess skin.

Introduction

Burns are debilitating, life threatening, and difficult to assess and manage (Lang et al. 2019). Complications after a burn can include, among others, a shock, infection, and long-lasting distress. Further, almost all deeper burns will lead to scarring. The post-burn scars may be immature/mature, atrophic/hypertrophic/keloid, stable/unstable, depigmented (vitiligo)/hyperpigmented and may turn malignant as well (Goel and Shrivastava 2010). In addition, scars are subject to change. For example, an immature scar can mature, and an atrophic scar can become hypertrophic. Post-burn scars are dry and itchy, and need to be prevented from exposure to sunlight.

One of the common complications in post-burn scars is contraction. Contraction is an active biological process that decreases an area of skin loss in an open wound because of a concentric reduction in the wound’s size (Goel and Shrivastava 2010). Starting in the proliferative phase of wound healing, wound contraction processes until full scar maturation, after which wound contraction can become active again. Depending on the extent of contraction and the wound dimensions, this can cause limited range-of-motion of joints if the scar is on or in a joint (contracture). This can lead to immobility and is an important indication for scar revision (Egberts et al. 2020). Scar reconstruction can be necessary.

Burns are so different from other types of wounds that there is a very separate discipline for this class in medical care. There are many different classifications for burns. For example, a burn can be thermal, electrical, or chemical. Burns come with a generalized increase in capillary permeability due to heat effect and damage, and this increase in capillary permeability is not seen in any other type of wound (Tiwari 2012). Burn wound healing consists of three overlapping phases: inflammation (reactive), proliferation (reparative), and maturation (remodeling). During inflammation, the wound is cleaned and prepared for further protection from bacterial infection. Subprocesses in proliferation are reepithelialization, angiogenesis, fibroplasia, and wound contraction. The last phase, in which the scar maturates and attains a balanced structure, can take years. This results in a scar that, on average, has 50% strength of unwounded skin (within three months), and 80% on the long-term (Enoch and Leaper 2008; Young and McNaught 2011).

Within the proliferative phase, fibroplasia encompasses the sub-processes that cause the restoration of the presence of fibroblasts and the production of a new extracellular matrix (ECM) in the injured area (Koppenol 2017). Fibroblasts can differentiate to myofibroblasts which are responsible for pulling forces in the skin and stimulate, like fibroblasts, both the production of the constituents of the new collagen-rich ECM and the release of matrix metalloproteinases (MMPs). The differentiation of fibroblasts to myofibroblasts is stimulated by transforming growth factor (Desmoulière et al. 1993). The fibroblasts, the myofibroblasts, and collagen deposition play an important role in the wound contraction.

For a long time, mathematical models have been developed that simulate the processes involved in wound healing. These models, (Barocas and Tranquillo 1997; Dallon et al. 1999; Koppenol 2017; McDougall et al. 2006; Olsen et al. 1995; Tranquillo and Murray 1992) to name a few, predict the behavior of experimental and clinical wounds, and gain insight into which elements of the wound healing response might have a substantial influence on the contraction. The majority of the models can be placed into one of three categories: continuum hypothesis-based models, discrete cell-based models, and hybrid models (Koppenol 2017). One subcategory of the continuum hypothesis-based models comprises the mechano-(bio)chemical models. These models, together with hybrid models, served as a basis for the morphoelastic model that we use in this study and that is developed by Koppenol and Vermolen (2017).

This model is based on the following principle (Hall 2008): the total deformation is decomposed into a deformation as a result of growth or shrinkage and a deformation as a result of mechanical forces. In a mathematical context, one considers the following three coordinate systems: , , and , which, respectively, represent the initial coordinate system, the equilibrium at time t that results due to growth or shrinkage, and the current coordinate system that results due to growth or shrinkage and mechanical deformation. Assuming sufficient regularity, the deformation gradient tensor is written byin which the tensor Z represents the deformation gradient tensor due to growth or shrinkage, and A represents the deformation gradient due to mechanical forces (Goriely and Amar 2006; Hall 2008; Rodriguez et al. 1994).

The beauty of this model is that we can simulate a permanent deformation, which results from the contraction process in burn wound healing. The main variable in this model is the displacement of the skin (u), i.e., the variable that makes us able to determine the surface of the wound, and in later stages, the scar. Besides these results, we can also determine the degree of ‘discomfort’ that the patient experiences. We do this by integrating over the entire tissue, including the undamaged part, yielding the total strain energy density. With this, we show to what extent there is an elongation in the entire tissue compared to the situation in which there would be no burn. This elongation, which in principle is simply determined by variations in local displacement, may signal nerves, which may cause the patient to experience a nagging sensation. This leads to discomfort in the patient.

The morphoelastic model compromises many parameters. We know some parameter values, while others are unknown and which we need to estimate. Although Koppenol has provided a great overview of parameter values, parameters vary between patients and even along a piece of skin sample. Hence we are interested in both the sensitivity of the parameters and the feasibility of the model. For the feasibility study, we might choose patient-specific information as input, such as the genetic background, gender, age, the location of the wound on the body, the depth of the wound, or any other. In our search for parameter values, we have seen that a great subset of parameters values is estimated for human skin tissue of different age. For example, the average fibroblast doubling time decreases with age (Simpson et al. 2009), indicating a decrease in fibroblast apoptosis rate with age. We therefore choose to use the patient’s age as input for the feasibility study. We summarize the change in parameter values that come with aging, and we use the sensitivity analysis results to vary the parameter values along the domain of computation. In our feasibility study, we define distinct groups of patients of different age for which we simulate many burns. The results show the variations in the relative surface area (RSA) density and the total strain energy density, both for patients of different age.

The organization of this paper is as follows. Section 2 presents the mathematical model and Sect. 3 presents the numerical method that we used to approximate the solutions. Subsequently, Sect. 4 presents the parameter values, Sect. 5 presents the sensitivity analysis, and Sect. 6 presents the feasibility study. Finally, Sect. 7 presents the conclusion and discussion.

The mathematical model

For the sake of completeness, we present the model that was also used in some of our earlier studies. We model the contraction in burn wounds and scars using partial differential equations that consider the displacement of the dermal layer (u), the displacement velocity of the dermal layer (v), the effective Eulerian strain present in the dermal layer (), and the changes in distributions of cells and densities of molecules. This morphoelastic continuum hypothesis-based modeling framework was introduced by Koppenol and Vermolen (2017). The model follows the assumption from Ramtani (2004); Ramtani et al. (2002), which states that the Young’s modulus of skin depends on the collagen concentration. The model incorporates the evolution in the distributions of fibroblasts (N) and myofibroblasts (M), the concentrations of signaling molecules (c) such as cytokines, chemokines, and growth factors, and the collagen concentration ().

We discuss the conservation laws for mass and linear momentum, and the evolution equation that describes how the infinitesimal effective strain changes. Since the domain deforms because of forces exerted by cells, the points in the domain are subject to displacement. We incorporate this local displacement rate by passive convection, which is reflected in the second term in the left-hand side of the equations.

The signaling molecules play an important role in the immune and inflammation response after wound healing. While these molecules migrate to the wound, they induce directed chemotactic migration of cells. We assume that the diffusion is according to normal Fickian diffusion. Furthermore, we take into account the enhanced secretion by fibroblasts and we assume that a portion of myofibroblasts (Barrientos et al. 2008) secrete and consume signaling molecules as well. It is known that MMPs remove signaling molecules from the dermal layer (Mast and Schultz 1996; Sternlicht and Werb 2001). The MMP concentration depends on the distribution of (myo)fibroblasts and the collagen concentration (Lindner et al. 2012), and the signaling molecule concentration (Overall et al. 1991). Hence taking the proteolytic breakdown by MMPs into account, we have the following:Here, is the Fickian diffusion coefficient of the signaling molecules, is the maximum net secretion rate of the signaling molecules, is the ratio of myofibroblasts to fibroblasts in the maximum secretion rate of the signaling molecules, is the concentration of the signaling molecules that causes the half-maximum net secretion rate of the signaling molecules, is the proteolytic breakdown rate parameter of the signaling molecules, is the ratio of myofibroblasts to fibroblasts in the secretion rate of the MMPs and represents the inhibition of the secretion of the MMPs. Note that the MMP balance is assumed to be instantaneous. One can find the derivation of this equation in the work of Koppenol and Vermolen (2017), and partly in the work of Olsen et al. (1995).

We assume that the migration of the cells through the wound bed and surrounding undamaged tissue is determined by random walk and by chemical stimulation by the signaling molecules (Boon et al. 2016; Dallon et al. 2001; Postlethwaite et al. 1987). This is modeled by a minimal model for chemotaxis (Hillen and Painter 2008), and a cell density-dependent Fickian diffusion. We incorporate the proliferation of the cells into the model using logistic growth models, dependent on the signaling molecules (as an activator–inhibitor (Murray 2011)), and inhibition because of crowding (Vande Berg et al. 1989). Further, the equations represent differentiation of fibroblasts to myofibroblasts (Tomasek et al. 2002), enhanced by signaling molecules, and apoptosis of the cells:Here, represents (myo)fibroblast random diffusion and is the chemotactic parameter that depends on both the binding and unbinding rate of the signaling molecules with its receptor, and the concentration of this receptor on the cell surface of the (myo)fibroblasts, is the cell division rate, is the maximum factor of cell division rate enhancement because of the presence of the signaling molecules, is the concentration of the signaling molecules that cause half-maximum enhancement of the cell division rate, represents the reduction in the cell division rate because of crowding, q is a fixed constant, is the signaling molecule-dependent cell differentiation rate of fibroblasts into myofibroblasts, is the apoptosis rate of fibroblasts and is the apoptosis rate of myofibroblasts. Myofibroblasts only proliferate in the presence of the signaling molecules, hence the difference between equation (3) and equation (4). Although myofibroblasts can differentiate back to fibroblasts under the influence of Prostaglandin () (Garrison et al. 2013), we do not take into account the re-differentiation of myofibroblasts to fibroblasts.

For collagen, we assume that there is no active transport present, because collagen molecules are large, which reduces their diffusivity. Since collagen is extracellular, it is, next to diffusion, not subject to further active migration mechanisms. Collagen is produced by (myo)fibroblasts (Baum and Arpey 2006), and enhanced by the secretion by signaling molecules (Ivanoff et al. 2005). The proteolytic breakdown of collagen is like for the signaling molecules related to the MMP concentration:Here, is the collagen secretion rate, is the maximum factor of secretion rate enhancement because of the presence of the signaling molecules, 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 is assumed always to be at a local equilibrium concentration. Reasoning for this modeling choice has been to avoid even more complexity and additional unknown parameter values.

In the equation for the displacement velocity, the Cauchy stress tensor is related to the effective Eulerian strain and displacement velocity gradients by a visco-elastic constitutive relation. The body force is generated by an isotropic stress and a pulling force on the extracellular matrix (ECM) by myofibroblasts, proportional to the product of the cell density of the myofibroblasts and a function of the concentration of collagen:Here, represents the total mass density of the dermal tissues, is the viscosity, represents the Young’s modulus (stiffness), is the generated stress per unit cell density and the inverse of the unit collagen concentration, R is a constant.

To incorporate plastic deformation in the equation for the effective strain (7), a tensor-based approach is used that is also commonly used in the context of growth of tissues (such as tumors). We assume that the rate of active change of the effective Eulerian strain is proportional to the product of the amount of effective Eulerian strain (Hall 2008), the local MMP concentration, and the local signaling molecule concentration:Here, is the rate of morphoelastic change (i.e., the rate at which the effective strain changes actively over time).

The domain of computation, initial and boundary conditions

We define the domain of computation by with , the closed interval. Similarly, we define the wounded area by the subspace , . Furthermore, we define the steepness of the boundary of the wound by s that counts for the slope of the constituents on the boundary of the wound. The dimension x is in centimeters and t in days.

We use the following functions for the initial fibroblast density and the initial signaling molecule concentration:, , . Here, are the fibroblast density and signaling molecule concentration in healthy dermal tissue and in the wound, respectively. Examples of possible initial densities are shown in Fig. 1. We use these functions in order to avoid steep changes in the densities. We have assumed that some fibroblasts are present in the wound and used that signaling molecules are present in the wound due to the secretion by, for instance, macrophages during inflammation.

Fig. 1

Example of the initial fibroblast density and initial signaling molecule density, with values of the parameters: cm, cm, cm, cells/, cells/, g/ and g/

In line with Koppenol, we assume that initially no myofibroblasts are present, that a fixed collagen concentration is present, and that the displacement of the dermal layer, the displacement velocity, and effective Eulerian strain initially are zero. Hence for all :We impose the following boundary conditions. For all and :The boundary condition for the displacement velocity follows from the assumption that the boundary of computation is sufficiently far away from the wound boundary. We impose no boundary conditions on the collagen concentration, the displacement of the dermal layer, and the effective Eulerian strain.

The numerical method

We approximate the solution to the model equations by the finite-element method using linear basis functions. For more information about this method, we refer to Van Kan et al. (2014). We multiply the equations (2)-(7) by a test function , integrate over the domain of computation (integration by parts), apply the application of the Gauss’ theorem, and apply the Leibniz–Reynold’s transport theorem.

To construct the basis functions, we subdivide the domain of computation into sub-domains (i.e., the elements). Let the finite element subspace and the vertices of the elements. We choose , as the linear basis functions, where denotes the Kronecker delta function.

Note that the following holds for the chosen subspace : for all (Dziuk and Elliot 2007). The Galerkin equations are simplified using this property. We solve the Galerkin equations using backward Euler time integration and we use a monolithic approach with inner Picard iterations to account for the non-linearity of the equations. To avoid loss of monotonicity (i.e., oscillations), we use mass lumping.

We approximate the local displacements of the dermal layer (u) with

Parameter values

In our search for parameter values, we used various sources from the literature. One of the most important sources is Koppenol’s thesis, in which several parameter values have been estimated that we did not find in existing literature. In our study, we conduct a sensitivity analysis and a feasibility study. In the Appendix, one finds the parameter values that we used in this study. In this chapter, we describe how we chose these values.

Equilibrium values

Taking into account, the reaction term for the signaling molecules and the equilibria and , the equilibrium density of the signaling molecules should be g/.

The estimation of the equilibrium distribution of fibroblasts differs per study. One estimates the number to be about cells/ (Olsen et al. 1995), and the other estimates the number to be about cells/ (Miller et al. 2003). The estimation of the number of cells also differs for the papillary and the reticular dermis, where there exist much more fibroblasts in the papillary dermis (Harper and Grove 1979; Randolph and Simon 1998). In our simulations, we have seen that the model works best with the equilibrium distribution of cells/. We note that some other parameter values () depend on the chosen order for , since we need to take into account the density of MMPs. Furthermore, research has found that among ages 1-10 the number of fibroblasts is nearly twice as high as in any other postnatal age group (Gunin et al. 2011). Therefore, we choose the mean value cells/ and let the value decrease with age.

The number of myofibroblasts present in the skin depends on the condition of the skin. Myofibroblasts result from the differentiation of fibroblasts. We assume that myofibroblasts are not present at the beginning of proliferation, since healthy skin contains almost no myofibroblasts. Hence, cells/.

Olsen et al. estimate the equilibrium collagen concentration as follows. Roughly 75% of the 15% of other substances than water and fat in 1 ml of human dermal tissue is collagen (Olsen et al. 1995). This yields g g . Furthermore, in human skin, the collagen content decreases at about 2% per year (Farage et al. 2015). Therefore, we choose the mean value g/ and let the value decrease exponentially with age.

Initial values

Because of growth factor supply in the inflammatory phase, the initial signaling molecule concentration is unequal to zero. The value should not exceed 15-50 ng (Olsen et al. 1995) and is therefore chosen to be g/.

The thermal injury causes sudden death of cells. The dead cells lose their solid integrity, which causes the release of cytokines. These cytokines trigger the immune response where several types of immune cells clear up the debris and release signaling molecules, which trigger the fibroblasts to migrate to the damaged region. Since we simulate from the onset of proliferation, we assume that several fibroblasts are present. We let this number be 20 percent of the equilibrium number. So the mean value is cells/.

Flux values

Sillman et al. vary the migratory rates of fibroblasts depending on the experimental medium used: in serum-containing medium the average velocity was as low as 0.23 mm/min, while in serum-free keratinocyte medium the average velocity was as high as 0.36 mm/min (Sillman et al. 2003). Hence, in serum-containing medium, the rate was /day and in serum-free keratinocyte medium the rate was /day. All the reported values together yield a mean value of /day and standard deviation /day. However, other estimates are /day and /day (Olsen et al. 1995) & (Simpson et al. 2017). We therefore estimate the value of /(cells day). Furthermore, we assume that the diffusion of (myo) fibroblasts decreases with age.

For the chemotactic parameter, we adopt /(cells day) from Murphy et al. (2012). For the diffusion parameter, we adopt /day from Haugh (2006). Furthermore, we assume that the diffusion of signaling molecules decreases with age.

Chemical kinetics values

Olsen et al. relate the inhibitor of TGF- to the initial concentration of the growth factors so that g/ (Olsen et al. 1995). We adopt this value.

Myofibroblasts produce roughly twice the collagen that is synthesized by fibroblasts Rudolph and Vande Berg (1991). Hence the constant .

The half-life of TGF- is about 2 minutes (Wakefield et al. 1990), and the half-life of PDGF is about 2 minutes as well (Bowen-Pope et al. 1984). So, signaling molecules have a decay rate of /day. However, Olsen et al. decrease the value for two reasons: not all signaling molecules may bind, for example because of insufficient levels of binding protein present at the wound site, and the bound complex may be recognized by (myo)fibroblasts leasing to internalized and metabolized signaling molecules (Olsen et al. 1995). Therefore, the estimated decay rate is 0.5/day. Other estimates for TGF- are /day (Javierre et al. 2009) and 0.354/day (Murphy et al. 2012; Yang et al. 1999). Given our equilibrium parameter values, the MMP density has order of magnitude . Hence, taking care of the equilibrium dimensions of the model, we end up with a range of /(cells g day). We take the value /(cells g day).

From our previous stability analysis, it follows that (Egberts et al. 2020). Given the parameter values, we set g/(cells day).

We estimate the constant , which is a small deviation from the constant estimated in Koppenol (2017).

Overall et al. estimate /g (Overall et al. 1991). We choose the lower limit, hence /g. Furthermore, the production of MMPs increases with age (Ashcroft et al. 1997). Given the equation for the MMPs (Koppenol and Vermolen 2017), that isto let the production of MMPs increase, we must decrease the inhibiting factor . Hence, we let the inhibition factor decrease with age.

Cell doubling time can be calculated using the growth rate (amount of doubling in one unit of time) in the following way: doubling time = ln(2)/growth rate. The average doubling time for fibroblasts is approximately 18–20 h (Alberts et al. 1989; Gosh et al. 2007). This gives the range for the proliferation rate of . We choose the upper limit, hence /( day). Furthermore, the percentage of PCNA-positive fibroblasts decreases with age, and PCNA can be considered as a marker for the proliferating cells (Gunin et al. 2011). We therefore let the cell division decrease with age.

TGF- increases fibroblast proliferation by 2–3 times (Strutz 2001). We choose the upper limit, hence .

The chemical concentrations required to enhance fibroblast proliferation are higher than those for chemotactic responses (Grotendorst 1992). Experimental evidence indicates that half-maximal enhancement corresponds to concentrations about 10 ng per ml (Olsen et al. 1995). We adopt this value and take g/.

The carrying capacity of fibroblasts is approximately /cells (Vande Berg et al. 1989). This value is adopted in this study. Furthermore, skin becomes thinner with age and therefore we assume that crowding occurs faster in elderly. Hence we let the division rate reduction value increase with age.

We need to have a stable chemical reaction in case the cell distributions and molecules densities are in equilibrium. The constant q gives us the opportunity to have a stable reaction in equilibrium for equation (3). Given the equilibria, solving for q yields:In Desmoulière et al. (1993), culturing fibroblasts in the presence of TGF- increased the percentage of cells expressing -SMA from 7.5 to 45.3%, representing an activation of 37.8% of myofibroblast type cells. This experiment occurred over a one week period, with a TGF- dose of 5–10 ng per ml. Suppose the activation of myofibroblasts follows a linear equation. Then given , we have /day. A dose of 5–10 ng per ml yields and /(g day), giving the range /(g day). We choose the upper limit. Furthermore, Simpson et al. demonstrated a failure of fibroblast-myofibroblast differentiation and showed that this is associated with in vitro aging (Simpson et al. 2009). Hence, we let the differentiating parameter decrease with age.

The average fibroblast doubling time (DT) ranges between 18 and 20 h (Alberts et al. 1989; Gosh et al. 2007), and the average lifespan of fibroblasts varies between 40 and 70 population doublings (PD) (Azzarone et al. 1983; Moulin et al. 2011). Using the formula we end up with the range . We choose the value /day and let this value decrease with age, since on average, the doubling time of fibroblasts decreases with age (Simpson et al. 2009).

The apoptosis rate of myofibroblasts was estimated in a previous study for hypertrophic scars (Koppenol et al. 2016). Within this study, it was found that a value of /day corresponds to hypertrophic scars and that a value of /day corresponds to normal scars. Further averages are: 8.85% for normal scars and 1.06% for hypertrophic scars (Moulin et al. 2003). Combination of these results yield the range for normal scars and for hypertrophic scars. For our study, we use the lower value /day for normal scars.

The secretion rate of collagen gives us the opportunity to have a stable reaction in equilibrium for equation (5). Given the equilibria, solving yields .

The synergistic effects of growth factors may accelerate collagen biosynthesis up to tenfold (Olsen et al. 1995). Hence, .

Data suggest that the half-maximal enhancement of collagen synthesis occurs at TGF- concentrations of the order of 1 ng per ml (Roberts et al. 1986). We adopt this value, hence g/.

For the decay rate of collagen, we let /(cells g day) (Koppenol et al. 2016). Furthermore, the collagen turnover decreases with age (Farage et al. 2015). Hence, we let the proteolytic breakdown of collagen decrease with age.

Mechanical values

Koppenol et al. estimated the viscosity value to be of order for the two-dimensional morphoelastic model (Koppenol and Vermolen 2017). In our previous study, the stability analysis showed that must hold for the one-dimensional morphoelastic model (Egberts et al. 2020). Given other parameter values we can adopt the value (N day)/. Furthermore, since the viscosity is constant for patients up to their 40s and increases a little after turning 40 (Xu and Tianjian 2011), we let the viscosity increase with age.

We estimate that the constant E in the Young’s Modulus is 350 N/((g cm)) for the one-dimensional morphoelastic model and let this value increase with age (Pawlaczyk et al. 2013; Pond et al. 2018).

For the parameters in the body force, we adopt (N g)/(cells ) (Maskarinec et al. 2009) & (Wrobel et al. 2002), g/, and /(cells g day). We set /(cells g day) and let this value increase with age, because the skin’s ability to recover after stretching decreases over lifetime (Krueger and Luebberding 2017).

Last, but not least, g/ for human skin (Wrobel et al. 2009). We assume this density does not change with age.

Sensitivity analysis

The model contains 34 parameters of which we vary the following independent 30 to study the sensitivity of these parameters:We also vary the length of the initial wound . The analysis is organized as follows. For each chosen parameter we vary the value by decreasing or increasing it by . This means we perform 341 simulations. Namely, for each parameter, we perform 11 simulations while leaving the values of the other parameters at the mean value. The mean values are given in Tables 3a–d, where in Tables 3c, d the mean values are given in the third column ().

Table 3

Parameter values. In subtables (c) and (d), denotes the mean value in group i of patients of different age

Parameter	μ	Dimension
(a)Parameter values that are constant among age groups and constant along domain
M¯	0	cells/cm^3
c¯	0	g/cm^3
c~	10^-8	g/cm^3
δ_M	6×10^-2	/day
R	0.995	g/cm^3
rFmax	2	–
η^I	2	–
η^II	0.45	–
kρmax	10	–

the equilibria and , and the initial conditions and ;

the apoptosis rates and , and the decay rates and ;

the parameters responsible for enhancement of cell division and molecule secretion and , and the inhibition of MMP secretion ;

the ratios from myofibroblasts to fibroblasts and , and chemokine dependent differentiation rate ;

the diffusion and chemotaxis rates and ;

the proliferation and secretion rates and , and the maximum factors and ;

the crowding factor ;

the parameters and R that influence the force;

the viscosity , Young’s Modulus factor E, morphoelastic factor , and the total mass density of dermal tissues .

We are interested in the contraction during wound healing and scar formation, and in the total strain energy density. The latter is assumed to be a measure for the discomfort that the patient experiences and is defined by:Here, we used the symmetry of the domain.

The results show the minimum of the relative surface area () (i.e., maximum contraction value) in a time of one year, the day on which the minimum relative surface area is reached () (i.e., the day after which the wound/scar retracts), the relative surface area on day 365 (), the maximum value of the strain energy density (), and the day on which the maximum value of the strain energy density is reached () (i.e., the day after which the patient experiences a reduction in discomfort due to the internal stress in the skin). In our simulations we used cm for the boundary of the domain of computation and the mean value cm for the wound.

Given the values in and for a variation , we compute the z-scores for the parameter . The basic z-score for a sample is , where is the sample mean and is the sample standard deviation.

We define the measure for sensitivity by summing over the absolute values of the z-scores: where is the z-score of the data in r for parameter i in variation j. For example, represents the z-score of the relative surface area on day 365 for parameter in the simulation where the value for is increased with 15%.

Table 1 gives an overview of the sensitivity values in terms of z-scores for the 31 parameters that we varied. In the last column, we rounded the sum of the values. From this table, we can see that the parameter that represents the equilibrium collagen concentration () with score 129 is the most sensitive. It is therefore interesting to study the equilibrium collagen concentrations in human skin, since collagen concentrations decrease with age (Farage et al. 2015) and we use this value for our age study in the next section. Parameters that are the least sensitive are the diffusion rate of (myo)fibroblasts () with score 10 and the initial collagen concentration () with score 11. Concerning the diffusion , we note that the mean value is of order , which differs from the order used by Koppenol (2017), where it is of order . This may lead to stating that variations of this parameter have no major impact on the simulations, while in a different geometry it might be much more sensitive. Concerning the initial collagen concentration (), we note that the value is varied when the equilibrium collagen concentration is fixed to the mean value (). In case the equilibrium collagen concentration is varied, the parameter for the initial collagen concentration is fixed to 20% of 0.1125, which is the mean value of and not the variation.

Table 1

Sensitivity of the varied parameters in terms of z-scores

Param.	Dimension	S^RSA_min	S^RSA_day	S^RSA_365	S^SED_max	S^SED_day	S^total
N¯	cells/cm^3	9.207	8.835	9.496	10.349	7.58	45
ρ¯	g/cm^3	20.351	33.916	15.383	29.19	30.21	129
c~	g/cm^3	3.687	3.015	2.065	2.721	4.544	16
N~	cells/cm^3	3.347	2.959	1.979	3.813	3.321	15
ρ~	g/cm^3	1.418	3.107	1.765	2.991	2.048	11
δ_c	g/cm^3	12.001	9.09	11.558	5.1	16.287	54
δ_N	/day	22.549	11.747	23.891	17.851	15.05	91
δ_M	/day	20.52	11.794	23.787	17.245	11.423	85
δ_ρ	cm^6/(cells g day)	8.905	15.125	6.879	13.054	13.552	58
acI	g/cm^3	4.828	2.057	4.779	2.823	4.499	19
acII	cm^3/g	8.929	10.753	8.912	11.091	8.917	49
acIII	g/cm^3	4.496	9.211	5.731	10.169	2.257	32
acIV	g/cm^3	2.306	3.107	1.529	2.991	3.037	13
η^I	–	8.175	3.364	8.227	4.821	8.248	33
η^II	–	6.551	5.713	9.24	6.304	5.607	33
D_c	cm^2/day	2.154	3.015	1.595	2.882	2.297	12
D_F	cm^5/(cells day)	1.954	2.057	1.604	1.867	2.427	10
χ_F	cm^5/(g day)	1.817	3.107	1.598	2.991	2.512	12
k_c	g/(cells day)	16.993	15.121	14.234	10.024	22.49	79
k_F	cm^3/(g day)	1.329	7.903	2.324	8.49	2.259	22
r_F	cm^3q/(cells^q day)	4.649	1.486	3.69	2.301	4.07	16
rFmax	–	17.546	9.879	18.152	13.054	12.097	71
kρmax	–	8.268	12.692	6.455	11.591	12.036	51
κ_F	cm^3/cells	1.89	3.107	1.573	2.991	2.551	12
ξ	(N g)/(cells cm^2)	14.892	2.057	14.571	1.18	11.62	44
R	g/cm^3	20.174	10.753	21.508	13.568	13.648	80
μ	(N day)/cm^2	1.782	3.107	1.6	2.991	2.423	12
E	N/((g cm)^1/2)	13.306	4.157	6.397	4.94	8.485	37
ζ	cm^6/(cells g day)	3.845	5.098	1.728	2.991	10.851	25
ρ_t	g/cm^3	1.817	3.107	1.598	2.991	2.512	12
L^w	cm	2.038	9.173	11.521	9.795	4.011	37

The third to seventh columns show the scores on the minimum of the relative surface area (RSA) , the day on which the minimum RSA is reached , the RSA on day 365 , the maximum value of the strain energy density (SED) , and the day on which the maximum of the SED is reached . The last column shows the total of the scores .

Other parameters that seem significantly sensitive are the apoptosis of fibroblasts () with score 91, the apoptosis of myofibroblasts () with score 85, the constant R that influences the force with score 80, and the secretion of signaling molecules () with score 79. However, we note that Koppenol estimated the mean values for , R and for a two-dimensional setting in which he used other values for the parameters (Koppenol 2017). Furthermore, the value of is based on a stability constraint that was found in our previous stability analysis (Egberts et al. 2020), and since the parameter for the equilibrium collagen concentration is sensitive, it is not a surprise that this secretion parameter is also sensitive. The value for the secretion rate of signaling molecules is not that straightforward. The secretion rate of cytokines differs from the secretion rate of growth factors and yet we model these together in one variable c representing signaling molecules. To prevent the model from unnecessary complicated computations, we continue modeling with this simplification and bear in mind the sensitivity of the parameter . It is therefore also interesting to study the rates of apoptosis of fibroblasts, since the doubling time of fibroblasts decreases with age (Simpson et al. 2009), and we use this value for our age study in the next section.

To get a visual insight into the sensitivity of the parameters, we present the effect of the variations on the parameters on both the relative surface area density of the wound/scar, and the strain energy density in Fig. 2. To have a clear distinction between the sensitivity of the parameters, there is no legend, and we labeled the most important lines with different styles. Within the figure, the sub-figures show the effects on the minimum of the relative surface area (a), the effects on the day on which the relative surface area reaches minimum (b), the effects on the relative surface area on day 365 (c), the effects on the maximum of the total strain energy density (d), and the effects on the day on which the total strain energy density reaches maximum (e).

Fig. 2

Effects of the variations in parameters for the relative surface area and the total strain energy density. Shown are the effects on the minimum of the relative surface area a, the effects on the day on which the minimum of the surface area is reached b, the effects on the relative surface area on day 365 c, the effects on the maximum of the total strain energy density d, and the effects on the day on which the maximum of the total strain energy density is reached e. In the figures, , , and are the fibroblast and myofibroblast apoptosis rates, and the signaling molecule and collagen decay rates, respectively, is the collagen equilibrium density, is the signaling molecule secretion rate, and R is a constant

From Fig. 2a, b, we see that the parameters that have the most influence on the minimum of the relative surface area and the day on which the relative surface area reaches minimum, are the apoptosis rates of the (myo)fibroblasts, the equilibrium collagen density, and the constant R that together with myofibroblasts and collagen is an important parameter in the body force in equation (6). In particular, a reduction in the fibroblast cell-density yields a reduction in the production of collagen, whereby a reduction in the skin’s stiffness (). Note that a smaller relative surface area corresponds with more intense contraction. Hence, from Fig. 2a, we see that there is less contraction for a smaller fibroblast apoptosis rate, a larger myofibroblasts apoptosis rate, a larger collagen equilibrium density, and a larger value for R. The latter three are no surprise, given the body force term in equation (6). Considering the fibroblast apoptosis rate, a reduction in apoptosis rate means more cells survive, hence a relative increase in proliferation. We can see this from equation (12): a smaller fibroblast apoptosis rate directly correlates with a smaller value for q, and given equation (3), we see that if q becomes smaller, the production of fibroblasts increases. This also relates to the skin’s stiffness, as in an increase in the fibroblast cell-density yields an increase in collagen production, and hence an increase in stiffness. The tissue’s strength therefore increases, and the effective strain decreases, i.e., becomes smaller. The result is less contraction during proliferation (Fig. 2a) and maturation (Fig 2c), where the effect during proliferation is larger. Further, from Fig. 2d, e we conclude that a decrease in fibroblast apoptosis rate results in less discomfort, over a longer period in time.

In Fig. 2b, the equilibrium collagen density stands out. We see that an increase in the equilibrium collagen density results almost linearly in a significant reduction in the time where the contraction is maximal. In contrast, a decrease of 25% relative to mean values results in a more exponential-like increase in the maximal contraction time. This property also relates to the total strain energy density, which we see in Fig. 2d, e. From the first, we see that the discomfort that a patient can experience increases more intensively for a 5–20% reduction in the myofibroblast apoptosis rate and the body force-inhibiting constant R. If the equilibrium collagen density decreases with 25%, then the effect is larger. In reality, it is not likely to change the equilibrium collagen density, however, we can use collagen dressings and, for example, vitamin C supplements to reduce contraction and healing time.

Figure 2b, c also feature the secretion rate of the signaling molecules. From the figures, we see that a lower signaling molecule secretion rate reduces both the period and intensity of contraction. Especially during maturation, when inflammatory responses are not favorable.

Taken together, targeting contraction intensity in the proliferative phase of wound healing is most likely effective in case fibroblast survival and collagen density are considered. Targeting contraction, and contractures, during maturation is more likely to be effective when inhibition of signaling molecules and collagen production are considered.

We use the results from the sensitivity analysis to perform a feasibility study in the next chapter.

Feasibility study of modeling age dependent scar/skin contraction

To study the feasibility of age-dependent uncertainty quantification, we focus on the effect of aging of skin on contraction, the final contracture, the total strain energy density, and the maximum of the total strain energy density. Just like any other organ, aging also affects the skin. Aging has a delaying effect on wound healing and immune responses. Intrinsic aging is the effect of generic and internal influences, such as hormones or metabolic substances. Extrinsic aging is the effect of external influences, such as UV radiation and environmental toxins (Wiegand et al. 2017). Clear general signs of aging are wrinkles, sagging skin and pigmentary irregularities, and increased tendencies to injuries and the faster opening of healing wounds. These symptoms result from physiological changes such as decreased cell replacement rate. We review various sources from literature to find suitable values for the parameters of the model. In this way, we can perform simulations for patients of different ages. Based on the results found, we have chosen the groups that are presented in Table 2.

Table 2

Groups of patients of different age

Group	Age
1	0–15
2	16–40
3	41–70
4	71+

In this study, there are five groups of parameters: To assess the uncertainty in the input data, we use a basic Monte Carlo method in which we sample input data from predefined statistical distributions. Regarding spatially heterogeneous parameters, we use sampling from a log-normal distribution. Each sample is a one-dimensional realization, and is based on the heterogeneous sampling through a normalized truncated Karhunen–Loéve expansion of a zero-mean stochastic process, by where , hence denotes a set of iid stochastic variables that follow the standard normal distribution, is the length of the domain of computation, and . From the stochastic variable , we show the regeneration of, for example, bytherewithHence, is a realization of a lognormal distribution with mean (expected value) and standard deviation . These values can be expressed by the arithmetic (sample) mean and arithmetic standard deviation :In the same way, we can create heterogeneous, stochastic inputs for other parameters as well.

parameters that are constant among age groups and not varied along the domain of computation,

parameters that are constant among age groups and varied along the domain of computation,

parameters that are varied between age groups and not varied along the domain of computation,

parameters that are varied between age groups and varied along the domain of computation,

parameters that depend on other parameters.

To test the model’s feasibility, we vary the parameter values based on the results found in literature on aging skin. We are interested in the differences in the intensity of contraction and the total strain energy between the distinct age groups. To quantify these differences, we test the null-hypothesis versus a two-sided alternative for groups A and B of patients using the following t-statistic:where and are the mean values of the results in distinct age groups A and B, is the estimated standard error of , and are the standard deviations in the age groups A and B, and is the number of samples in the age groups. Here, we assume that the number of samples in the age groups are equal. We reject the null-hypothesis if , with .

To reduce the computation time, we performed simulations on half a domain with cm. In all simulations cm, hence the initial wounded area is given by .

We simulated burns per age group. Hence, in total we simulated 7800 burns. We used parallel computing with 12 processors, three processors responsible per group, on a 64 bit Windows 10 Pro system with 16 GB RAM and 3.59 GHz AMD Ryzen 5 3600 6-Core Processor. The total computation time was 13.5 hours; hence, the mean computation time per simulation takes less than half a minute. A major advantage of the one-dimensional implementation is its short computation time, which allows to do many of them within a reasonable time interval. For the test statistic, we used . The standard deviations for the age groups for the minimum of the relative surface area (i.e., maximum contraction) areThe standard deviations for the age groups for the relative surface area on day 365 areTo get some insight into the effects of the four groups of patients of different age shown in Table 2, we present confidence intervals for the contraction of burns in different age groups in Fig. 3. In Figs. 4 and 5, we show probability density functions of the minimum relative surface area and the relative surface area on day 365, respectively, together with the corresponding cumulative distribution function. Finally, we present confidence intervals for the total strain energy in the healing of burns in different age groups in Fig. 6, and probability density functions of the maximum of the total strain energy density and its cumulative distribution function in Fig. 7. We computed the probability density functions in Figs. 4a, 5a and 7a using the kernel density estimation method. The values for the parameters of the model, together with the standard deviation values for the variation over the domain of computation (if applicable), are shown in Table 3 in the Appendix.

Fig. 3

Confidence intervals for the contraction of burns in different age groups. The intervals show the mean values and the 95% confidence values of the mean

Fig. 4

Estimated probability density function a and cumulative distribution function b of the minimum relative surface area

Fig. 5

Estimated probability density function a and cumulative distribution function b of the relative surface area on day 365

Fig. 6

Confidence intervals for the strain energy in the healing of burns in different age groups. The intervals show the mean values and the 95% confidence values of the mean

Fig. 7

Estimated probability density function a and cumulative distribution function b of the maximum total strain energy density

Figure 3 shows four 95% confidence intervals for the mean of the size of the scar, each confidence interval corresponding to a group of patients. The range of the contraction values comes from the variability of the parameters over the domain of computation. We can see that the maximum contraction value is about the same in the first two groups of ages, and from group 2, a higher age group gives a larger reduction of the size of the scar, and therewith a larger intensity of the contraction. Further, for higher ages, it takes more time to reach the maximum intensity of contraction.

There seems to be more variability in the permanent deformation in the elderly patients. In the elderly patients, it takes longer before the wound healing cascade reaches equilibrium than in younger patients. Minima of the relative surface area were mostly reached on days 34, 61, 74, and 95, with values of 76.7, 76.0, 74.4, and 72.9% for groups 1, 2, 3, and 4, respectively. Unfortunately, these results do not correspond fully to the observations in the clinic. Normally, contraction is of less order in elderly patients and, in general, retraction takes a longer period. We note that the longer retraction period is visible in the two-dimensional model results by Koppenol and Vermolen (2017), and that this ‘stretched retraction period’ is handled with the parameter . Furthermore, in clinic one sees more contractures in younger patients in case the injury was in or near a joint. In contrast, in the clinic, the elderly seem to experience less discomfort because of contraction. One reason for this could be that the skin of the elderly is less tight than the skin of the young. Looser tissue can move more compared to skin that is already tight and is therefore less likely to cause movement restriction when it contracts. The reason our simulation results do not meet the clinical observations is because of the variety of factors that have not been modeled in our mathematical model yet.

Figure 4 shows the results on the minimum relative surface area (i.e., maximum of contraction). Although we see that there is an overlap between all the groups in Fig. 4a, the maximal contraction is significantly different () between the groups (see Table 4(a) in the Appendix). The differences in the minimum of the relative surface area between consecutive age groups are largest between ages 16–40 (group 2) and 41–70 (group 3), and smallest between ages 0 and 15 (group 1) and 16–40 (group 2). Given that in reality the evolution of the size of the scar differs from our results, we expect that the differences between group 1 and group 2 are larger in reality, possibly largest.

Table 4

The values of the t-statistic comparing the minima of the relative surface area (a), the maxima of the total strain (b) and the relative surface area on day 365 (c)

Groups	t-value
(a) The values of the t-statistic comparing the minima of the relative surface area between age groups
1&2	17
1&3	58
1&4	94
2&3	43
2&4	80
3&4	36
(b) The values of the t-statistic comparing the maxima of total strain between age groups
1&2	23
1&3	63
1&4	91
2&3	49
2&4	83
3&4	32
(c) The values of the t-statistic comparing the relative surface area on day 365 between age groups
1&2	167
1&3	273
1&4	385
2&3	126
2&4	241
3&4	108

The overlap between the age groups is also visible in the estimated cumulative distribution function plot in Fig. 4b, where the functions of the first and second group, and the third and fourth group, almost intersect. From the cumulative distribution functions, we can estimate the probabilities of reaching a certain amount of maximum contraction. For example, this figure suggests that with 70% probability a patient from groups of patients 1, 2, 3, and 4, respectively, can reach 22.6%, 23.4%, 24.9%, and 26.4% contraction. Using such functions in the future can help to predict the probability of developing a contracture, to intervene in time when possible.

Figure 5 shows the results on the relative surface area on day 365 (i.e., permanent contraction). Eventually, the scar maturates and because of the morphoelastic behavior of the skin, the size of the scar almost never reaches its initial configuration again. Here, the initial configuration represents the initial size and geometry of the burn wound. Like for minimal contraction values, the overlap between consecutive groups for the intensity of contraction after one year in Fig. 5a is significant (, see Table 4(b) in the Appendix). The differences in the relative surface area after one year between consecutive groups are similar to the differences in the minimum of the relative surface area. The probability density functions show a correlation between aging and the spread, confirming the observation found in the confidence intervals, that there is larger variability in the intensity of contraction in elderly people than in children.

From the cumulative distribution functions in Fig. 5b, we can estimate the probability of a certain contraction intensity. Here we see mean values for the final contraction about 97.6%, 95.9%, 94.5%, and 93.3% for groups 1, 2, 3, and 4, respectively. This final contraction intensity is an indicator of a possible contracture, in case the scar is in or near a joint. Given the location of the scar, one might estimate probabilities that this scar will develop a contracture of certain extent.

Figure 6 shows four confidence intervals for the mean of the total strain energy density. The total strain energy density is a measure for the discomfort that a patient experience. Each confidence interval corresponds to a group of patients of different age. The range of the total strain energy density comes from the variability of the parameters over the domain of computation. We can see that the maximum of the total strain density is about the same in the first two groups, and from group 2, a higher age group gives a larger maximum of the total strain energy density. Further, for higher ages, it takes more time to reach the maximum total strain energy density. Note that all these results relate to the relative surface area densities shown in Fig. 3.

The total strain energy densities reach maxima on days 36, 64, 78, and 95, with values of 62, 66, 74 and 79 for groups 1, 2, 3, and 4, respectively. We note that the total strain energy density reaches maxima a few days later than the maximum contraction in almost all groups. Because the contraction data in Fig. 3 does not relate to what doctors see in the clinic, we assume this is the same case for the total strain energy density. This means that in real-life children might experience more discomfort than the elderly, in contrast to what we see in Fig. 6.

Figure 7 shows the results on the maximum of the total strain energy density. Like in Figures 4 and 5, although the probability density functions in Fig. 7a show an overlap between all the groups, the maximum of the total strain energy density is significantly different () between the groups (see Table 4(b) in the Appendix). This overlap is also visible in the cumulative distribution function plot in Fig. 7b, where the functions of the first and second group intersect on the top, and the third and fourth group almost intersect. From the cumulative distribution functions, we can estimate the probabilities of reaching a certain amount of maximum total strain. For example, this figure suggests that with 80% probability a patient from groups of patients 1, 2, 3, and 4, respectively, can reach a maximum of 67.9, 70.4, 78.5, and 83.6 total strain. Since the figures show a strong correlation between the contraction and total strain, we conclude that targeting the contraction directly targets the total strain as well.

Conclusion and discussion

In this study, we worked with the morphoelastic model for contraction and contractures in burn scars. This model was developed by Koppenol, and based on the principle of morphoelasticity (Hall 2008). We have provided a one-dimensional version of this model and focused on the parameter values, the sensitivity of the parameter values, and the feasibility of the model for patients of different ages.

We comprehensively described the (ranges of the) values of the parameters. Most of the variety in the parameter values we have found in literature sources. We estimated some parameter values and adopted other parameter values from Koppenol. In case there were ranges of values found in literature, we chose upper or lower bounds, or a fixed value in between. We have used these values as mean values for the sensitivity analysis, and for the parameters we have used for the second age group in our feasibility study.

For our sensitivity analysis, we varied 30 parameters by %, and we varied the length of the initial wound. We showed results for the minimum of the relative surface area, the day on which wound healing reached this minimum, the relative surface area on day 365, the maximum value of the strain energy density, and the day on which wound healing reaches this maximum. The most sensitive parameter is the equilibrium collagen concentration present in the dermal layer. Other parameters that seem significantly sensitive are the apoptosis rate of fibroblasts, the apoptosis rate of myofibroblasts, the constant R in Eq. (6) that influences the body force, and the secretion rate of signaling molecules. The parameters we found to be the least sensitive are the fluxes parameters, the crowding factor, the viscosity, the mass density of dermal tissue, and the initial collagen concentration. We note that we let the initial collagen concentration depend on the equilibrium collagen concentration, which can influence the sensitivity value of this parameter.

Furthermore, we performed a feasibility study for the model to investigate the effect of aging on contraction, contractures, and discomfort in burn wound healing. We have chosen for four groups of patients in age groups: 0–15, 16–40, 41–70, 70+. We varied the parameters of the model according to observations from the literature, so that there was a variation between the groups of patients. We furthermore varied the parameters of the model using Karhunen–Loéve expansions, to model the heterogeneity of human skin, and in the framework of our Monte Carlo method we performed sampling from statistical distributions to assess the impact of uncertainty in the data on the behavior of contraction. The model is feasible for this approach, showing increased extent of contraction with age, a delayed maximum amount of contraction in elderly people (showing delayed healing), increased contracture in elderly people, and increased variety of contracture formation in elderly people compared to children. Next to these results, we see that the extent of discomfort is highly related to the contraction in wound healing. The figures show that there is a larger amount of discomfort in elderly patients and that the maximum discomfort is experienced significantly () earlier in younger children than in other age groups.

This study shows that contraction increases with age and shows that there is a significant difference () in maximum amount of contraction between the different age groups. The least significant difference is found between ages 0–15 and 16–40, which can also be seen in Fig. 3. Further, the differences in the amount of contraction on day 365 in consecutive groups is least significant between 41–70 and 70+ years (see Table 4c in the Appendix). The most significant difference in the amount of contraction on day 365 in consecutive groups is found to be between 0–15 and 16–40 years. Given the cumulative distribution functions, we can give probabilities of a certain amount of contraction in specific age-dependent groups of patients. We have seen that the differences in the maximum contraction and the contraction on day 365 are of a few percentages (less than 10) of order. For the maximum discomfort that a patient might feel, we have seen that there is significant difference () between all groups, of which in consecutive groups the difference between ages 16–40 and 41–70 is most significant. From the figures, we can conclude that these patients experience the same amount of discomfort, although this happens much quicker in children.

However, we note that we obtained the results with a mathematical model for which it is hard to find validation data. In the clinic, contraction seems to be of less order in elderly patients, and the retraction takes a longer period. The elderly have excess skin, meaning that they suffer less from contraction. There is less tension and less stretch in the skin of the elderly and hence we assume that there is less contraction. An explanation for our results might be that in fact such an extent of contraction is present in the elderly, however, this is not seen and noticed because of above reasons. On the other hand, the process of contraction is affected by growth and mobility. These two kinds of forces yield a different process of contraction in agile, growing children compared to grown ups. Generally, children seem to suffer less from contractures because of the common chest burn. However, when a child burns near a joint, the contractures seem to be more intense than in grown ups. This is because the growth induces extra forces on the damaged tissue, which keeps the myofibroblasts active. Therefore, the myofibroblasts remain in the scar to pull on the surrounding tissue. For the age groups in between the children and the elderly, we are not sure if our results are representative. There is a lack of patient-specific data because these types of studies are insufficiently performed. Therefore, we conclude that the model could be feasible for an age-related study for contraction in dermal wound healing and for that, the model needs further adjustments to account for children’s growth and the excess skin of elderly patients.

In this study, we defined the discomfort a patient may experience as the total strain energy. In equation (14), we defined the total strain energy density for epsilon between and . However, it could well be that the definition where we take epsilon between and is more appropriate as a measure of the discomfort the patient is experiencing. We plan to consider this option in future studies, where we take into account two-dimensional effects.

Further research is needed to understand the differences of extent and timing of contraction in patients. To account for growth and movement forces, we plan on adding a factor for growth in equation (7), and wound boundary forces for movement. To simulate for aging patients, we might put time-dependence on (some of) the parameters and simulate for a much longer time period. For example, the skin collagen content decreases at about 2% per year (Farage et al. 2015), which has the most influence on the time healing period and amount of contraction. To distinguish between superficial and severe wounds, we can adapt real-life behavior in the initial conditions of the variables. In fact, in wounds where the environment still is vital, the release of cytokines is because cells lose their solid integrity. However, in severe wounds with a damaged environment, the release happens because of two reasons. Namely, the injury boiled signaling molecules, and the environment is dead. Here, circulation is absent and therefore little happens. Cytokines are present, but from vital edges where circulation still occurs, penetrating the tissue. The question is how deep, with an answer we can consider in the initial conditions. To incorporate more geometrical matters, it is necessary to extend the model to a more-dimensional framework. Though more-dimensional frameworks allow to assess geometrical issues, these more-dimensional frameworks will require more simulation time and the use of a more advanced computer infrastructure if the objective is to carry out Monte Carlo simulations for the assessment of the likelihood that contraction of a particular intensity occurs. To reduce computational times, we want to model the boundaries of the wound as elastic springs. Further, it would be an improvement to code the finite element solution to the model in a high level programming language such as C++ and use an artificial intelligence framework such as neural networks. In finite element analysis, the necessary mesh generation procedure is widely regarded as the weakest link in the chain of the analysis. To avoid the failure of a mesh to be analysis suitable, we want to use isogeometric analysis (IGA) in the more-dimensional framework. We plan to assess these issues in future work.