Biomedical implications from a morphoelastic continuum model for the simulation of contracture formation in skin grafts that cover excised burns.

A continuum hypothesis-based model is developed for the simulation of the (long term) contraction of skin grafts that cover excised burns in order to obtain suggestions regarding the ideal length of splinting therapy and when to start with this therapy such that the therapy is effective optimally. Tissue is modeled as an isotropic, heterogeneous, morphoelastic solid. With respect to the constituents of the tissue, we selected the following constituents as primary model components: fibroblasts, myofibroblasts, collagen molecules, and a generic signaling molecule. Good agreement is demonstrated with respect to the evolution over time of the surface area of unmeshed skin grafts that cover excised burns between outcomes of computer simulations obtained in this study and scar assessment data gathered previously in a clinical study. Based on the simulation results, we suggest that the optimal point in time to start with splinting therapy is directly after placement of the skin graft on its recipient bed. Furthermore, we suggest that it is desirable to continue with splinting therapy until the concentration of the signaling molecules in the grafted area has become negligible such that the formation of contractures can be prevented. We conclude this study with a presentation of some alternative ideas on how to diminish the degree of contracture formation that are not based on a mechanical intervention, and a discussion about how the presented model can be adjusted.

Introduction

In the United Kingdom, approximately 250,000 citizens get injured due to burning each year (Hettiaratchy and Dziewulski 2004). In the United States, about half a million citizens require medical treatment as a result of thermal injury each year (Gibran et al. 2013). The majority of these burns are minor and do not require specialized care. However, a small portion of the injuries is extensive and as a consequence roughly 13,000 individuals in the United Kingdom and approximately 40,000 individuals in the United States, are admitted to a hospital or burn center for treatment each year (Gibran et al. 2013; Hettiaratchy and Dziewulski 2004).

The core treatment of burns in these medical centers consists usually of two subparts; first most of the burnt skin is excised surgically and thereafter the newly created wound is covered by a skin graft. The use of a skin graft to cover a newly created wound has two widely recognized benefits compared to the situation where these wounds are left to heal by secondary intention alone; in general it reduces both the overall contraction of the grafted area and the development of hypertrophic scar tissue in these areas (Walden et al. 2000). Unfortunately, however, many times skin grafts still contract considerably after placement on their recipient bed and this may result then in substantial shrinkage of the grafts and hence the development of contractures in these tissues (Kraemer et al. 1988).

The development of contractures is a serious complication that has a significant impact on an affected person’s quality of life, and often requires substantial further corrective surgery (Leblebici et al. 2006). Therefore, therapies have been developed which aim at the prevention of the formation of contractures. The main therapy in current usage focuses on counteracting the mechanical forces generated within the contracting graft by means of static splinting of the covered wound after placement of the graft (Richard and Ward 2005). How effective splinting therapy is in preventing contracture formation is actually unclear at the moment; it is a fact that contracture formation is still a common complication despite the frequent application of splinting therapy (Schouten et al. 2012). This could be a consequence of the fact that it is unclear at present what the optimal point in time is after placement of the skin graft to start with the therapy. Furthermore, it could also be a consequence of the fact that it is also unclear how long the static splints have to be worn for the therapy to be effective.

This unsatisfactory situation is probably partly caused by the fact that actually little is known about the etiology of the formation of contractures (Harrison and MacNeil 2008). We think that a better understanding of the mechanism underlying contracture formation probably aids in the development of a better treatment plan that reduces the development of this sequela, and argue that computational modeling studies can contribute to the expansion of this understanding. Therefore we develop here a new mathematical model for the simulation of the contraction of skin grafts that cover excised burns in order to gain new insights into the mechanism underlying the formation of contractures. Based on the obtained insights, we give suggestions regarding the ideal length of splinting therapy and when to start with the therapy such that the therapy is effective optimally. In addition, we put forward some alternative ideas on how to diminish the degree of contracture formation that are not based on a mechanical intervention.

The development of the model is presented in Sect. 2. Subsequently, the simulation results are presented in Sect. 3. Here, we also show good agreement with respect to the evolution over time of the surface area of unmeshed skin grafts that cover excised burns between outcomes of computer simulations obtained in this study and scar assessment data gathered previously in a clinical study. The model and the simulation results are discussed in Sect. 4.

Development of the mathematical model

Given that contraction mainly takes place in the dermal layer of skin tissues, we incorporate solely a portion of this layer into the model. The layer is modeled as a heterogeneous, isotropic, morphoelastic continuous solid with a modulus of elasticity that is dependent on the local concentration of the collagen molecules. With respect to the mechanical component of the model, the displacement of the dermal layer (), the displacement velocity of the dermal layer (), and the infinitesimal effective strain present in the dermal layer () are chosen as the primary model variables (The latter variable represents a local measure for the difference between the current configuration of the dermal layer and a hypothetical configuration of the dermal layer where the tissue is mechanically relaxed (See also Eq. (3))). Furthermore, we select the following four constituents of the dermal layer as primary model variables: fibroblasts (N), myofibroblasts (M), a generic signaling molecule (c), and collagen molecules ().

In order to incorporate the formation of contractures (i.e., the formation of long term deformations) into the model, we use the theory of morphoelasticity developed by Hall (2009). Central to this theory is the assumption that the classical deformation gradient tensor (i.e., ) can be decomposed into a product of two tensors (i.e., ) (Hall 2009; Goriely and Ben Amar 2007; Rodriguez et al. 1994). The tensor can be thought of as the locally defined deformation from the fixed reference configuration to a hypothetical configuration (i.e., a zero stress state (Fung 1993)) wherein the internal stresses around all individual points in the dermal layer are relieved, and the tensor can be thought of as the locally defined deformation from this hypothetical configuration to the current configuration of the dermal layer. Based on this decomposition, Hall derived several related evolution equations that describe mathematically the change of the effective strain over time. Hence, these equations basically give a mathematical description of the remodeling of the dermal layer over time. In this study we assume that the effective strains are small. Therefore, we use here the evolution equation that describes the dynamic change of the infinitesimal effective strain over time (i.e., Eq. (1c) in this study and Eq. (5.64) in the PhD thesis of Hall (2009)). (The derivation of this evolution equation is rather lengthy and contains numerous subtleties. Therefore, we present here solely the finally derived equation. The full derivation of the evolution equation can be found in the PhD thesis of Hall.)

Combined with the general conservation equations for mass and linear momentum in local form, the following continuum hypothesis-based framework is used as basis for the model: where andEquation (1a) is the conservation equation for the cell density / concentration of constituent i of the dermal layer, Eq. (1b) is the conservation equation for the linear momentum of the dermal layer, and Eq. (1c) is the evolution equation that describes how the infinitesimal effective strain changes over time. Within the above equations represents the cell density / concentration of constituent i, represents the flux associated with constituent i per unit area due to random dispersal, chemotaxis and other possible fluxes, represents the chemical kinetics associated with constituent i, represents the total mass density of the dermal tissues, represents the Cauchy stress tensor associated with the dermal layer, represents the total body force working on the dermal layer, is the displacement velocity gradient tensor (i.e., ), and is a tensor that describes the rate of active change of the effective strain. The operator is the Jaumann time derivative, and the operator is the material time derivative. (If the material time derivative is applied to the effective strain tensor, then it is applied to each of the scalar elements of this tensor separately.) Given the chosen primary model variables, we have . In order to simplify the notation somewhat we replace by i in the remainder of this study. Hence, becomes N, becomes M, and so on.

The cell populations

The functional forms for the biochemical kinetics associated with the (myo)fibroblasts and the functional forms for the movement of these cells are identical to functional forms used previously (Koppenol et al. 2017b). For the biochemical kinetics, the functional forms are whereThe parameter is the cell division rate, is the maximum factor with which the cell division rate can be enhanced due to the presence of the signaling molecule, is the concentration of the signaling molecule that causes the half-maximum enhancement of the cell division rate, represents the reduction in the cell division rate due to 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. The functional forms for the cell fluxes are where is the cell density-dependent random motility coefficient of the (myo)fibroblasts, and is the chemotactic coefficient.

The generic signaling molecule

The functional form for the net production of the generic signaling molecule (i.e., the first term on the right hand side of Eq. (10)) and the functional form for the dispersion of the signaling molecule are identical to functional forms used previously (Koppenol et al. 2017b): The parameter represents the maximum net secretion rate of the signaling molecule, is the ratio of myofibroblasts to fibroblasts in the maximum net secretion rate of the signaling molecule, is the concentration of the signaling molecule that causes the secretion rate of the signaling molecule to be half of its maximum, and is the proteolytic breakdown rate of the signaling molecules. The parameter represents the random diffusion coefficient of the generic signaling molecule. An example of a signaling molecule that can stimulate processes such as the up-regulation of the secretion of collagen molecules by (myo)fibroblasts and the cell differentiation of fibroblasts into myofibroblasts, is transforming growth factor- (TGF-) (Barrientos et al. 2008).

The second term on the right hand side of Eq.(10) requires a more detailed introduction. In this study, we incorporate into the model the proteolytic cleavage of the generic signaling molecule by metalloproteinases (MMPs) (Mast and Schultz 1996; Lint and Libert 2007). MMPs are secreted by (myo)fibroblasts and are involved in the breakdown of collagen-rich fibrils during the maintenance and the remodeling of the extracellular matrix (ECM) (Chakraborti et al. 2003; Lindner et al. 2012; Nagase et al. 2006). The secretion of the MMPs is inhibited by the presence of signaling molecules such as TGF- (Overall et al. 1991). Therefore, we assume that the concentration of the MMPs is a function of the cell density of the (myo)fibroblasts, the concentration of the collagen molecules, and the concentration of the signaling molecules:The parameter is the ratio of myofibroblasts to fibroblasts in the secretion rate of the MMPs, and represents the inhibition of the secretion of the MMPs due to the presence of the signaling molecule.

The collagen molecules

The functional forms for the biochemical kinetics associated with the collagen molecules and the functional form for the transportation of these molecules are basically identical to functional forms used previously (Koppenol et al. 2017b): The parameter is the collagen molecule secretion rate, is the maximum factor with which the secretion rate can be enhanced due to the presence of the signaling molecule, is the concentration of the signaling molecule that causes the half-maximum enhancement of the secretion rate, and is the proteolytic breakdown rate of the collagen molecules.

The mechanical component

In this study, we use the following visco-elastic constitutive relation for the mathematical description of the relationship between the Cauchy stress tensor on the one hand, and the effective strains and displacement velocity gradients on the other hand: Here is the shear viscosity, is the bulk viscosity, is Poisson’s ratio, is the Young’s modulus, and is the second-order identity tensor. Like Ramtani et al. (2004; 2002), we assume that the Young’s modulus is dependent on the concentration of the collagen molecules. The parameter is a fixed constant.

Furthermore, we incorporate into the model the generation of an isotropic stress by the myofibroblasts due to their pulling on the ECM. This pulling stress is proportional to the product of the cell density of the myofibroblasts and a simple function of the concentration of the collagen molecules (Olsen et al. 1995; Koppenol et al. 2017a, b): The parameter represents the total generated stress by the myofibroblast population, is the generated stress per unit cell density and the inverse of the unit collagen molecule concentration, and R is a fixed constant.

Finally, we assume that the rate of active change of the effective strain is proportional to the product of the amount of effective strain (as suggested by Hall (2009)), the local concentration of the MMPs, the local concentration of the signaling molecule, and the inverse of the local concentration of the collagen molecules. The directions in which the effective strain changes, are determined by both the signs of the eigenvalues related to the effective strain tensor, and the directions of the associated eigenvectors. Taken together, we obtain the following symmetric tensor:where is the rate of morphoelastic change (i.e., the rate at which the effective strain changes actively over time).

The domain of computation

We assume , and hold within the modeled portion of dermal layer for all time t, with the yz-plane running parallel to the surface of the skin andFurthermore, we assume that the derivatives of the cell densities and the concentrations of the modeled constituents of the dermal layer are equal to zero in the direction perpendicular to the surface of the skin. Taken together, these assumptions imply that the calculations can be performed on an arbitrary, infinitely thin slice of dermal layer oriented parallel to the surface of the skin, and that the results from these calculations are valid for every infinitely thin slice of dermal layer oriented parallel to the surface of the skin. Therefore, we use the following domain of computation:where are Lagrangian coordinates.

The initial conditions and the boundary conditions

The initial conditions give a description of the cell densities and the concentrations immediately after placement of the skin graft on its recipient bed. For the generation of the simulation results, the following function has been used to describe the shape of the skin graft:whereHere corresponds to grafted dermis and corresponds to unwounded dermis. The values for the parameters and determine, respectively, the location of the boundary between the skin graft and the undamaged dermis, and the minimum distance between completely grafted dermis and unwounded dermis. Furthermore, with the counterclockwise rotation matrix that rotates vectors by an angle about the X-axis, and .

Based on the function for the shape of the skin graft, we take the following initial conditions for the modeled constituents of the dermal layer:Here , and are, respectively, the equilibrium cell density of the fibroblasts, the equilibrium cell density of the myofibroblasts, and the equilibrium concentration of the collagen molecules, of the unwounded dermis. Due to the secretion of signaling molecules by for instance leukocytes, signaling molecules are present in the wounded area. The constant represents the maximum initial concentration of the signaling molecule in the grafted area. Furthermore, we assume that there are some fibroblasts present in the grafted area. The value for the parameter determines how much fibroblasts are present minimally initially in the grafted area.

With respect to the initial conditions for the mechanical component of the model, we take the following initial conditions for all where is the initial domain of computation in Eulerian coordinates:See Fig. 1 for a graphical representation of the initial conditions that have been used in this study.

Fig. 1

A graphical overview of the initial conditions. Depicted are the initial shape of the skin graft and, in color scale, the initial cell density of the fibroblasts (). The scale along both axes is in centimeters. The X-axis points toward the reader. The black dots mark the material points that were used to trace the evolution of the surface area of the skin graft over time. That is, at each time point, the area of the polygon with vertices located at the displaced black material points has been determined

With respect to the boundary conditions for the constituents of the dermal layer, we take the following Dirichlet boundary conditions for all time t and for all where is the boundary of the domain of computation in Eulerian coordinates:The parameter is the equilibrium concentration of the signaling molecule in the unwounded dermis.

Finally, with respect to the boundary condition for the mechanical component of the model, we take the following Dirichlet boundary condition for all time t and for all :

Table 1

An overview of the dimensional (ranges of the) values for the parameters of the model

Parameter	Value	Dimensions	Reference
D_F	10^-7	cm^5/(cells day)	Sillman et al. (2003)
χ_F	2×10^-3	cm^5/(g day)	Murphy et al. (2012)
q	-4.2×10^-1	−	NC
r_F	9.24×10^-1	cm^3q/(cells^qday)	Ghosh et al. (2007)
rFmax	2	−	Strutz et al. (2001)
acI	10^-8	g/cm^3	Grotendorst (1992)
κ_F	10^-6	cm^3/cells	Vande Berg et al. (1989)
k_F	1.08×10^7	cm^3/(g day)	Desmoulière et al. (1993)
δ_N	2×10^-2	/day	Olsen et al. (1995)
δ_M	6×10^-2	/day	Koppenol et al. (2017b)
D_c	2.9×10^-3	cm^2/day	Murphy et al. (2012)
k_c	4×10^-13	g/(cells day)	Olsen et al. (1995)
η^I	2	−	Rudolph and Vande Berg (1991),
			Moulin et al. (1998)
acII	10^-8	g/cm^3	Olsen et al. (1995)
δ_c	5×10^-4	cm^6/(cells g day)	Olsen et al. (1995)
η^II	5×10^-1	−	TW
acIII	(2-2.5)×10^8	cm^3/g	Overall et al. (1991)
k_ρ	6×10^-8	g/(cells day)	NC
kρmax	10	−	Olsen et al. (1995)
acIV	10^-9	g/cm^3	Roberts et al. (1986)
δ_ρ	6×10^-6	cm^6/(cells g day)	Koppenol et al. (2017b)
ρ_t	1.02	g/cm^3	Liang and Boppart (2010)
μ_1	10^2	(N day)/cm^2	TW
μ_2	10^2	(N day)/cm^2	TW
E^I	3.2×10	N/((g cm)^1/2)	Liang and Boppart (2010)
ν	4.9×10^-1	−	Liang and Boppart (2010)
ξ	5×10^-2	(N g)/(cellscm^2)	Maskarinec et al. (2009)
			&Wrobel et al. (2002)
R	9.95×10^-1	g/cm^3	TW
ζ	(0-9)×10^2	cm^6/(cells g day)	TW
N¯	10^4	cells/cm^3	Olsen et al. (1995)
M¯	0	cells/cm^3	Olsen et al. (1995)
c¯	0	g/cm^3	NC
ρ¯	10^-1	g/cm^3	Olsen et al. (1995)
I^w	2×10^-1	−	Koppenol et al. (2017b)
c^w	10^-8	g/cm^3	Olsen et al. (1995)

The last column contains the references to the articles that were used for obtaining (estimates of) the values for the parameters. If (the range of) the value for a parameter was estimated in this study, then this is indicated by the abbreviation TW. If the value for a parameter is a necessary consequence of the values chosen for the other parameters, then this is indicated by the abbreviation NC

The parameter value estimates

Table 1 in Appendix 3 provides an overview of the dimensional (ranges of the) values for the parameters of the model. The majority of these values were either obtained directly from previously conducted studies or estimated from results of previously conducted studies. In addition, we were able to determine the values for three more parameters due to the fact that these values are a necessary consequence of the values chosen for the other parameters (Koppenol et al. 2017b).

Simulation results

In order to obtain some insight into the dynamics of the model, we present an overview of simulation results for the modeled constituents of the dermal layer in Fig. 2. Furthermore, we present an overview of simulation results for the displacement field and the displacement velocity field in Fig. 3, and an overview of simulation results for the effective strain in Fig. 4. For the generation of these overviews, the same set of values for the parameters of the model was used.

Fig. 2

An overview of simulation results for the modeled constituents of the dermal layer when the inhibition of the secretion of MMPs due to the presence of signaling molecules is relatively low () and the rate of morphoelastic change is relatively high (). The values for all other parameters are equal to those depicted in Table 1 in Appendix 3. The top two rows show the evolution over time of the cell density of, respectively, the fibroblast population and the myofibroblast population. The color scales represent the cell densities, measured in . The bottom two rows show the evolution over time of the concentrations of, respectively, the signaling molecules and the collagen molecules. The color scales represent the concentrations, measured in . Within the subfigures, the scale along both axes is in centimeters

Fig. 3

An overview of simulation results for the displacement field and the displacement velocity field when the inhibition of the secretion of MMPs due to the presence of signaling molecules is relatively low () and the rate of morphoelastic change is relatively high (). The values for all other parameters are equal to those depicted in Table 1 in Appendix 3. The top two rows show the evolution over time of the displacement in, respectively, the horizontal direction and the vertical direction. The color scales represent the displacements, measured in centimeters. The bottom two rows show the evolution over time of the displacement velocity in, respectively, the horizontal direction and the vertical direction. The color scales represent the displacement velocities, measured in . Within the subfigures, the scale along both axes is in centimeters. The black squares within the subfigures represent the (displaced) boundaries between the skin graft and the unwounded dermis

Fig. 4

An overview of simulation results for the effective strain when the inhibition of the secretion of MMPs due to the presence of signaling molecules is relatively low () and the rate of morphoelastic change is relatively high (). The values for all other parameters are equal to those depicted in Table 1 in Appendix 3. The separate rows show the evolution over time of the different components of the effective strain that are unequal to zero. The color scales represent the amount of effective strain. Within the subfigures, the scale along both axes is in centimeters. The black squares within the subfigures represent the (displaced) boundaries between the skin graft and the unwounded dermis

Figure 2 shows that the cell density of the myofibroblasts, and the concentrations of both the signaling molecules and the collagen molecules increase first within the skin graft. Subsequently, the concentrations of these molecules, just like the cell density of the myofibroblasts, start to decline until they reach the equilibrium concentrations and the equilibrium cell density of uninjured dermal tissue. Meanwhile, the cell density of the fibroblasts starts to increase within the skin graft until it reaches the equilibrium cell density of uninjured dermal tissue.

Figure 3 shows that the boundaries between the skin graft and the uninjured tissue are pulled inward increasingly toward the center of the skin graft while the concentration of the collagen molecules and the cell density of the myofibroblasts increase. Looking at the displacement velocity field, we observe that the boundaries are pulled inward relatively fast initially. Subsequently, the speed with which the boundaries are pulled inward diminishes fast. Looking carefully at the displacement velocity field, we observe that the inward movement actually reverses from a certain time point onward. It is nice to observe that this phenomenon coincides with the gradual increase in the surface area of the skin graft, and the gradual decrease in both the cell density of the myofibroblasts and the concentration of the collagen molecules within the skin graft, as can be observed in, respectively, Fig. 6 and Fig 2. Furthermore, we observe that the boundaries between the skin graft and the uninjured tissue hardly move anymore eventually (i.e., the individual components of the displacement velocity field become approximately equal to zero over the domain of computation), and that the surface area of the skin graft has diminished considerably after a year. This latter phenomenon is also clearly visible in Fig. 6.

Fig. 6

The evolution over time of the relative surface area of wounds (i.e., skin grafts) for some combinations of values for the rate of morphoelastic change (i.e., the parameter ), and the parameter related to the inhibition of the secretion of MMPs due to the presence of signaling molecules (i.e., the parameter ). The values for all other parameters are equal to those depicted in Table 1 in Appendix 3

Figure 4 also shows something very interesting. If we look at the effective strain at day 365, we observe that the individual components of the effective strain tensor are not equal to zero over the domain of computation. This implies that there are residual stresses present in the grafted area. Comparing the properties of the effective strain at day 180 with the properties of the effective strain at day 365, we observe that these are more or less the same. Hence, the residual stresses remain present in the modeled portion of dermal layer for a prolonged period of time.

Figure 5 shows the evolution over time of the relative surface area of skin grafts for particular combinations of values for two parameters that are directly related to the tensor (See Eq. (19)). In addition, the figure shows averages of clinical measurements over time of the relative surface areas of placed unmeshed skin grafts in human subjects after both early excision of burnt skin and late excision of burnt skin (El Hadidy et al. 1994).

Fig. 5

The evolution over time of the relative surface area of wounds (i.e., skin grafts) for particular combinations of values for the rate of morphoelastic change (i.e., the parameter ), and the parameter related to the inhibition of the secretion of MMPs due to the presence of signaling molecules (i.e., the parameter ). The values for all other parameters are equal to those depicted in Table 1 in Appendix 3. The black circles and the black squares show the evolution over time of the average of clinical measurements of the relative surface areas of placed unmeshed skin grafts after, respectively, early excision of burnt skin and late excision of burnt skin (El Hadidy et al. 1994)

Furthermore, Fig. 6 shows the evolution over time of the relative surface area of skin grafts for some more combinations of values for the aforementioned parameters related to the tensor . The figure shows that both an increase in the rate of morphoelastic change (i.e., the parameter ), and an increase in the inhibition of the secretion of MMPs due to the presence of signaling molecules (i.e., an increase in the value for the parameter ) results in a reduction of the final surface area of a skin graft. Within the chosen ranges for the values of the parameters, we observe that a change in the value for the rate of morphoelastic change has a large impact on the final surface area of a skin graft. Changing the value for the parameter related to the inhibition of the secretion of MMPs due to the presence of signaling molecules has a smaller impact on the final surface area of a skin graft. Note also that the value for the latter parameter has a relatively large impact on the total number of days that the boundaries between the skin graft and the uninjured tissue are pulled inward after placement of the skin graft before the retraction process starts.

Finally, it is nice to observe in Fig. 6 that, as expected, the surface area of a skin graft returns to its initial value when the rate of morphoelastic change is equal to zero. If this rate is equal to zero, then the tensor is equal to the zero tensor. In this case, one would expect an initial period during which the surface area of a skin graft diminishes due to the pulling action of the myofibroblasts, followed by a period during which this surface area slowly returns to its initial value due to the apoptosis of the myofibroblasts. This is exactly what can be observed in the figure.

Discussion

We have presented a continuum hypothesis-based model for the simulation of the (long term) contraction of skin grafts that cover excised burns. Since skin contraction and contracture formation mainly take place in the dermal layer of the skin, we incorporated solely a portion of this layer into the model. The dermal layer is modeled as a heterogeneous, isotropic, morphoelastic solid with a Young’s modulus that is locally dependent on the concentration of the collagen molecules. For this end, we used the theory of morphoelasticity developed by Hall (2009). In particular, we used in this study the derived evolution equation that describes the dynamic change of the infinitesimal effective strain over time. Furthermore, we used the general conservation equations for linear momentum and mass to describe mathematically the dynamic change over time of, respectively, the linear momentum, and the cell densities and concentrations of the modeled constituents of the dermal layer. For the description of the relationship between the Cauchy stress tensor on the one hand, and the effective strain tensor and displacement velocity gradients on the other hand, we used the visco-elastic constitutive relation given in Eq. (15).

Related to the mechanical component of the model, we want to remark the following. Traditionally, the dermis is modeled as a linear visco(elastic) solid in mechano-chemical continuum models for dermal wound healing (Javierre et al. 2009; Murphy et al. 2012; Olsen et al. 1995; Ramtani 2004; Ramtani et al. 2002; Valero et al. 2014a, b; Vermolen and Javierre 2012). More recently, continuum models have appeared where the dermis is modeled as a hyperelastic solid (Koppenol et al. 2017a; Valero et al. 2013, 2015). Unfortunately, it is difficult with any of these models to simulate the long term deformation of dermal tissues and the development of residual stresses within these tissues while these phenomena are often observed in the medical clinic (Schouten et al. 2012). Therefore, we adopted like Murphy et al. (2011) and Bowden et al. (2016), a morphoelastic framework in this study. With the application of such a framework, it becomes relatively simple to simulate both the long term deformation of a skin graft and the development of residual stresses within the modeled portion of dermal layer.

With respect to the constituents of a recovering injured area, we selected the following four constituents as primary model variables: fibroblasts, myofibroblasts, a generic signaling molecule, and collagen molecules. The mathematical descriptions for the movement of the cells, the biochemical kinetics associated with these cells, the dispersion of the generic signaling molecule, and the release, consumption, and removal of both the collagen molecules and the generic signaling molecule are nearly identical to the functional forms used previously (Koppenol et al. 2017b).

Furthermore, we present an overview of the applied numerical algorithm that has been developed for the generation of computer simulations in Appendix 1. The development of this algorithm was necessary to “catch” the local dynamics of the model and obtain sufficiently accurate simulations within an acceptable amount of CPU time. For this end, we combined a moving-grid finite-element method (Madzvamuse et al. 2003) with an element resolution refinement / recoarsement method (Möller et al. 2008) and an automatically adaptive time-stepping method (Kavetski et al. 2002). We present the derivation of the general finite-element approximation in Appendix 2. Furthermore, we applied both a source term splitting procedure (Patankar 1980) and a semi-implicit flux-corrected transport (FCT) limiter (Möller et al. 2008) on the discretized system of equations that describes the dynamics of the modeled constituents of the dermal layer in order to guarantee the positivity of the approximations of the solutions for these primary model variables.

With the developed model, it is possible to simulate some general qualitative features of the healing response that is initiated after the placement of a skin graft on its recipient bed (Harrison and MacNeil 2008). The restoration of the presence of fibroblasts within the skin graft and the temporary presence of myofibroblasts during the execution of the healing response can be simulated. Due to the initial presence of signaling molecules and the gradual increase in the cell density of the myofibroblasts in the grafted area, the secretion rate of collagen molecules is considerably larger than the proteolytic breakdown rate of these molecules in the grafted area for a prolonged period of time (See also Eq. (14)). Consequently, the concentration of the collagen molecules in the grafted area becomes substantially higher than the concentration of the collagen molecules in the surrounding uninjured dermal tissue before it gradually decreases toward the concentration of the collagen molecules in the surrounding uninjured dermal tissue. Furthermore, it is possible to simulate both the long term contraction and subsequent retraction of a skin graft, and the development of residual stresses within the dermal layer. These phenomena can be observed, respectively, in Figs. 3 and 4; both the displayed components of the displacement field and the displayed components of the effective strain tensor are not equal zero over the domain of computation at day 365, and the values of the individual components over the domain of computation at day 365 are roughly equal to the values of the individual components over the domain of computation at day 180. Looking at the individual components of the displacement velocity field in Fig. 3, it can be observed that these have become approximately equal to zero over the domain of computation at day 365.

Focusing on the simulation of the contraction of skin grafts and the formation of contractures we observe the following. Figure 5 shows a good match with respect to the evolution over time of the relative surface area of skin grafts between measurements obtained in a clinical study by El Hadidy et al. (1994) and outcomes of computer simulations obtained in this study. This agreement provides us some confidence about the validity of the model. Obviously, the number of models with which it is possible to produce the depicted contraction curves is infinite in theory. Therefore, we would have liked to validate the presented model against scar assessment data of a different kind such as cell density profiles and collagen molecule concentration profiles, in order to increase our confidence about the validity of the model. However, we have not been able to find more appropriate experimental measurement data in the available literature. We are not the only ones who have to deal with this issue. Unfortunately, it is a fundamental problem in the field of mathematical modeling of dermal wound healing processes to find suitable experimental measurement data for the proper validation of models (Bowden et al. 2016). In our opinion, this does not imply that we should refrain from deducing biomedical implications from the results obtained in this study. However, we do think that it is very important to be careful when doing so, and to keep in mind that these deductions are based on outcomes of a mathematical modeling study.

Having said that, we focus now on the implications of the results depicted in Fig. 5. In this study, we assumed that the rate at which the effective strain is changing actively over time is proportional to the product of the amount of effective strain, the local concentration of the MMPs, the local concentration of the signaling molecule, and the inverse of the local concentration of the collagen molecules. The directions in which the effective strain changes, are determined by both the signs of the eigenvalues related to the effective strain tensor, and the directions of the associated eigenvectors. The good match between the gathered scar assessment data and the outcomes of the computer simulations suggests that this combination of relationships might describe appropriately in mathematical terms the mechanism underlying the formation of contractures.

If the mathematical description for the mechanism underlying the formation of contractures is indeed appropriate, then this suggests the following. Looking at Eq. (19), it is clear that the effective strain can change solely when the local concentration of the signaling molecules is unequal to zero. Given the presence of signaling molecules within the grafted area immediately after placement of the skin graft on its recipient bed, this implies that the optimal point in time to start with splinting therapy is directly after surgery. It is interesting to note that this implication matches nicely with the finding that early mechanical restraint of tissue-engineered skin leads to a reduction in the extent of contraction (Harrison and MacNeil 2008). Furthermore, it is also evident that it is desirable to continue with splinting therapy until the concentration of the signaling molecules in the grafted area has become negligible such that the formation of contractures can be prevented. Given that it is unclear at present what the optimal point in time is after surgery to start with splinting therapy, and how long the static splints have to be worn for the therapy to be effective (Schouten et al. 2012), these are interesting observations.

Furthermore, Fig. 5 shows that the difference in the evolution over time of the average of the relative surface areas of placed unmeshed skin grafts between grafts that are placed after early excision of the burnt skin and grafts that are placed after late excision of the burnt skin might be caused by both a change in the rate of morphoelastic change and a change in the degree of inhibition of the secretion of MMPs due to the presence of signaling molecules. In itself, this is an interesting observation. In addition, it provides us with some alternative ideas on how to diminish the degree of contracture formation that are not based on a mechanical intervention. As demonstrated in Fig. 6, the final surface area of a skin graft can be increased by both a reduction in the rate of morphoelastic change, and a reduction in the inhibition of the secretion of MMPs due to the presence of signaling molecules. Within the investigated ranges of values, the former reduction has a huge impact on the final surface area whereas the latter reduction has a smaller impact on the final surface area. Perhaps that the reduction in the rate of morphoelastic change can be accomplished through the local inhibition of certain cross-linking enzymes. For instance, perhaps it is possible to use the chemical -aminopropionitrile (BAPN) for the inhibition of the cross-linking enzyme lysyl oxidase, which is crucial for the stabilization of collagen fibrils (Kagan and Li 2003; Wilmarth and Froines 1992). The reduction in the inhibition of the secretion of MMPs due to the presence of signaling molecules can be accomplished perhaps by influencing the regulation of the transcription of MMPs by signaling molecules such as TGF- (Overall et al. 1991). Given that it is actually unclear at present how effective splinting therapy is in preventing contracture formation (Schouten et al. 2012), we think that these suggestions, which are not based on a mechanical intervention, could be interesting alternatives for the reduction of the degree of contracture formation.

We want to conclude this section with a short discussion about how the presented model can be adjusted. We want to mention here a couple of ways in which the model can be adjusted. In this study, we made some assumptions which made it possible to perform the calculations on an infinitely thin slice of dermal layer oriented parallel to the surface of the skin. An obvious benefit of this is that it results in a serious reduction of the computation times to obtain computer simulations. However, it is probably very interesting to investigate in a three-dimensional portion of dermal layer what would happen to the (long term) contraction and subsequent retraction of skin grafts. For this end, the main adaptations to the model presented in this study would be the removal of the assumptions made in Sect. 2.5, and the introduction of additional boundary conditions for the interfaces between the epidermis and the dermis and the subcutaneous tissue and the dermis. In a three-dimensional setting, we could use, for example, the spring-like boundary conditions introduced in the study by Koppenol et al. (2017b) to describe the attachment of the dermal tissues to the subcutaneous tissue.

Unfortunately, there is no information presented about the shapes and the absolute surface areas of the excised burns in the study by El Hadidy et al. (1994). Therefore, we decided to use the same shape for the grafted area in all computer simulations. Obviously, it is easy to obtain computer simulations with the presented model for grafted areas with a different shape and / or surface area. The only element of the model that needs to be adjusted is the function given in Eq. (22). In particular, when more detailed clinical measurement data about the evolution over time of the shapes of skin grafts become available in the scientific literature, it will become very interesting to investigate to what extent the presented model in this study can replicate such clinical measurement data. Furthermore, the presented model could also be used for the simulation of the healing of dermal wounds that heal by secondary intention alone (i.e., without the placement of a skin graft). The only element of the model that would really need to be adjusted is the initial condition for the collagen molecules presented in Eq. (24).

In this study we assumed that the effective strains are small. As a consequence, we used in this study the evolution equation that describes the dynamic change of the infinitesimal effective strain over time (i.e., Eq. (1c)). However, it might actually be more appropriate to assume that the effective strains can become arbitrary large. If we make this assumption, then we can replace Eq. (1c) with the following evolution equation that gives a description of the dynamic change of the Eulerian finite effective strain () over time (Hall 2009):whereandSubsequently, we can replace the constitutive relation presented in Eq. (15) with, for instance, the following constitutive relation:Here the elastic component of the constitutive relation is equal to constitutive relation for a heterogeneous, isotropic and compressible neo-Hookean solid (Koppenol et al. 2017b; Treloar 1948).

The second way in which the mechanical component of the model can be adjusted is the following. It is known that the way that collagen molecules are organized into interconnected sheets and bundles influences the response of dermal tissues to mechanical forces (Jor et al. 2011). Therefore, we have developed a continuum hypothesis-based model in which the bulk mechanical behavior of the involved dermal tissues is dependent on the geometrical arrangement of the collagen bundles (Koppenol et al. 2017a). In this model, a tensorial approach is used to represent the collagen bundles, and the bulk mechanical properties of the tissues such as the Young’s moduli and Poisson ratios are dependent on both the local concentration and the local geometrical arrangement of these collagen bundles. Hence, the model we have presented here can be adjusted by using the tensorial approach for the representation of collagen bundles, and by replacing the elastic component of the constitutive relation given in Eq. (15) with the constitutive relation derived in the study by Koppenol et al. (2017a). Due to these adaptations to the model presented in this study, it becomes possible to study the impact of tissue anisotropy on the deformation of skin grafts during healing.