A biomechanical mathematical model for the collagen bundle distribution-dependent contraction and subsequent retraction of healing dermal wounds

A continuum hypothesis-based, biomechanical model is presented for the simulation of the collagen bundle distribution-dependent contraction and subsequent retraction of healing dermal wounds that cover a large surface area. Since wound contraction mainly takes place in the dermal layer of the skin, solely a portion of this layer is included explicitly into the model. This portion of dermal layer is modeled as a heterogeneous, orthotropic continuous solid with bulk mechanical properties that are locally dependent on both the local concentration and the local geometrical arrangement of the collagen bundles. With respect to the dynamic regulation of the geometrical arrangement of the collagen bundles, it is assumed that a portion of the collagen molecules are deposited and reoriented in the direction of movement of (myo)fibroblasts. The remainder of the newly secreted collagen molecules are deposited by ratio in the direction of the present collagen bundles. Simulation results show that the distribution of the collagen bundles influences the evolution over time of both the shape of the wounded area and the degree of overall contraction of the wounded area. Interestingly, these effects are solely a consequence of alterations in the initial overall distribution of the collagen bundles, and not a consequence of alterations in the evolution over time of the different cell densities and concentrations of the modeled constituents. In accordance with experimental observations, simulation results show furthermore that ultimately the majority of the collagen molecules ends up permanently oriented toward the center of the wound and in the plane that runs parallel to the surface of the skin.


Introduction
The overall healing of deep dermal wounds that cover a large surface area is an extremely complex process (Enoch and Leaper 2007;O'Toole and Mellerio 2010). Looking more closely at the healing process, the process can be split up into an interdependent series of relatively simpler subprocesses. One of these subprocesses is wound contraction. Wound contraction is a biomechanical process that causes the circumferential inward movement of surrounding uninjured tissue toward the wounded area (Baum and Arpey 2005;Monaco and Lawrence 2003). Due to this inward movement of the uninjured tissue, the exposed surface area of the wound can be decreased substantially and relatively fast without the production of new tissues. For instance, due to wound contraction, typical full-thickness wounds in humans may undergo a reduction in wound surface area of up to 40 % over a period of several weeks, while rapid wound contraction in rats may even be responsible for up to 90 % of the closure of the wounded area (Li et al. 2007;McGrath and Simon 1983).
Given the possible extent of wound contraction and the severity of related complications such as the development of permanent shortenings of scar tissue, a lot of resources have been allocated to the research on wound contraction over the last few decades. This research has resulted in the production of much knowledge about the biomechanical mechanisms underlying wound contraction. However, there is still much that remains understood incompletely about both the mechanisms underlying wound contraction and the etiology of the complications that may develop as a result of wound contraction. Given that fully adequate treatment plans for the prevention of, for instance, the development of permanent shortenings of scar tissue do not exist at the moment and a better understanding of the mechanisms underlying wound contraction most probably would aid in the development of better treatment plans, this is an unsatisfactory situation.
In order to help with gaining more insight into the mechanisms underlying wound contraction, Tranquillo and Murray (1992) formulated the first mathematical framework for the modeling of the mechanical component of the wound healing response. This framework is based on the continuum hypothesis and the conservation of both mass and linear momentum. In the resulting model, the dermal tissues are modeled as homogeneous, isotropic, linear viscoelastic solids and the fibroblasts that are present in these tissues produce an isotropic stress that works on these tissues. Interestingly, Tranquillo and Murray were able to replicate with this model the experimental data on wound contraction in rats that were collected by McGrath and Simon (1983). Subsequently, the general framework has been extended and has been adapted in several more recent modeling studies in order to investigate the impact of the addition and the adaptation of various different components of the wound healing response (Friedman et al. 2012;Javierre et al. 2009;Koppenol et al. 2016;Murphy et al. 2012;Olsen et al. 1998aOlsen et al. , 1995Tracqui et al. 1995;Ramtani 2004;Ramtani et al. 2002;Valero et al. 2013Valero et al. , 2014Vermolen and Javierre 2012).
Recently, continuum hypothesis-based models have been formulated that are based on a morphoelastic framework (Bowden et al. 2016;Murphy et al. 2011). With these models, it is possible to simulate simultaneously both the contraction and the growth of involved dermal tissues during the execution of the overall healing process. This combination makes it possible to simulate the permanent deformation of these tissues and the development of residual stresses within these tissues.
Although the mathematical models in the above cited studies contain functional descriptions of several different components of the wound healing response, they all lack a description of the dynamic regulation of the geometrical arrangement of the constituents of the extracellular matrix (ECM). Given that the behavior of dermal tissues in response to mechanical forces such as wound contraction is influenced strongly by the geometrical arrangement of the constituents of the ECM (Wilkes et al. 1973), we think it is a limitation of these models that they lack such a description.
The ECM is the non-cellular component of dermal tissues and is composed of two classes of macromolecules: proteoglycans and fibrous proteins (Jarvelainen et al. 2009;Schaefer and Schaefer 2010). Taken together these macromolecules form a relaxed network of protein fibers embedded in a hydrated gel. The most abundant type of fibrous protein is collagen. These proteins provide most of the tensile strength to the tissues they are embedded in (Rozario and DeSimone 2010). The majority of the collagen molecules is produced by fibroblasts and is organized into interconnected sheets and bundles by these same fibroblasts (Monaco and Lawrence 2003). Due to this organization of the collagen molecules into interconnected sheets and bundles, the geometrical arrangement of the collagen bundles, in particular, has a huge impact on the response of dermal tissues to mechanical forces (Jor et al. 2011). Barocas and Tranquillo (1997) formulated the first mathematical framework that includes a description of the dynamic regulation of the geometrical arrangement of collagen bundles in tissues. This framework is also based on the continuum hypothesis and the conservation of both mass and linear momentum. The ECM is modeled as a heterogeneous, anisotropic, biphasic medium consisting of a fibrillar network and an interstitial solution. The two phases display intraphase viscoelasticity and interphase frictional drag due to the relative motion of the two phases. The biphasic medium both orients and is oriented by fibroblasts. Toward a validation of the resulting model, Barocas and Tranquillo demonstrated qualitative agreement between model predictions on the one hand and various outcomes of experiments performed on tissue-equivalent systems on the other hand.
Subsequently, Olsen et al. (1998bOlsen et al. ( , 1999 also formulated a continuum hypothesis-based general framework that includes a description of the dynamic regulation of the geometrical arrangement of either one type or two types of bundle. In this framework, the ECM is modeled as a heterogeneous, anisotropic, viscoelastic medium in which the bundle type(s) both orient(s) and is (are) oriented by fibroblasts. Dallon et al. (1999Dallon et al. ( , 2000Dallon et al. ( , 2001 and McDougall et al. (2006) also developed a theoretical framework in order to study the dynamic regulation of the arrangement of collagen bundles. However, they took a different approach: cells are modeled as discrete entities, while the ECM is modeled as a continuum. Here too, the direction of movement of the cells is influenced by the local orientation of the collagen bundles, while the reorientation of the bundles is in turn dependent of the polarity of the cells in the vicinity. A couple of years ago, this hybrid framework has been extended by Cumming et al. (2010). More recently, we have extended the model devel-oped by Cumming et al. on our part by adding a mechanical component to this model (Boon et al. 2016).
Although the frameworks and the associated models mentioned in the last two paragraphs are very elegant, they either do not contain a mechanical component at all or they lack an incorporation of the effect of the geometrical arrangement of the collagen bundles on the bulk mechanical behavior of the involved dermal tissues. Hence, these models are not able to simulate the direct influence of the composition and the topology of the constituents of the ECM on the behavior of dermal tissues in response to mechanical forces such as wound contraction.
Recently, two models have been developed that do incorporate the effect of the geometrical arrangement of the collagen bundles on the bulk mechanical behavior of the tissues these bundles are embedded in (Valero et al. 2015;Yang et al. 2013). However, in our opinion, each of these two models has a serious limitation. Due to the fact that the model developed by Yang et al. uses a hybrid framework similar to that developed by Dallon et al. (1999), the domain of computation, and hence the created wounds, has to be small so that the computation times and the computer memory requirements remain acceptable. Therefore, they used a domain of computation of solely 4 mm 2 and created circular wounds with a radius of 400 µm. However, if dermal wounds cover such a small surface area, then wound contraction actually does not play a very substantial role during the healing process. Because of the fact that the model developed by Valero et al. is a continuum hypothesis-based model, it is possible to simulate the healing of deep dermal wounds that cover a large surface area with this model. However, this latter model lacks a dynamic regulation of the geometrical arrangement of the collagen bundles. Given that, for instance, the proportion of the collagen bundles that runs parallel to the surface of the skin in general increases considerably due to the execution of the wound healing processes (Welch et al. 1990), we consider it a serious limitation that the dynamic regulation of the geometrical arrangement of the collagen bundles is not included in the model developed by Valero et al. Because of the fact that we wish to simulate the contraction of deep dermal wounds that cover a large surface area where the contraction process is influenced by the geometrical arrangement of the collagen bundles, and given the aforementioned limitations associated with the models developed by Yang et al. and Valero et al., we developed a new, fully continuum hypothesis-based model for the simulation of the contraction and the subsequent retraction of healing dermal wounds. In this model, the bulk mechanical behavior of the involved dermal tissues is dependent on the geometrical arrangement of the collagen bundles. For this end, we use a tensorial approach similar to the one proposed by Barocas and Tranquillo (1997) and Cumming et al. (2010) to repre-sent the collagen bundles, and we let the bulk mechanical properties of the tissues such as the Young's moduli and Poisson ratios depend locally on both the local concentration and the local geometrical arrangement of these bundles. Furthermore, we incorporate into this new model the dynamic change of the geometrical arrangement of the collagen bundles similar to how this process is incorporated in the model developed by Olsen et al. (1998b). A detailed description of the full model is presented next in Sect. 2. In Sect. 3, we give a short overview of the applied numerical algorithm for obtaining simulation results. The simulation results are presented in Sect. 4. Finally, the model and the simulation results are discussed in Sect. 5.

Development of the mathematical model
In order to simulate the contraction and subsequent retraction of healing dermal wounds, we include into the model some of the subprocesses that take place during the proliferative and remodeling phase of the wound healing cascade (Enoch and Leaper 2007). With respect to the subprocesses that are executed during the proliferative phase, we select the following subprocesses: wound contraction and fibroplasia. (Fibroplasia encompasses the subprocesses that cause the restoration of the presence of fibroblasts and the restoration of a collagenrich ECM in the injured area.) Since wound contraction mainly takes place in the dermal layer of the skin, solely a portion of this layer is included explicitly into the model.
The dermis is modeled as a heterogeneous, anisotropic continuous solid with bulk mechanical properties that are locally dependent on both the local concentration and the local geometrical arrangement of the collagen bundles. Due to the fact that the collagen bundles are represented by means of a symmetric tensor, the dermal layer has material properties that differ locally along three mutually orthogonal twofold axes of rotational symmetry. (The symmetry axes coincide with the lines that pass through the individual material points of the dermal layer and run parallel to the individual eigenvectors of the tensor that represents the collagen bundles. See Sect. 2.3 for further details.) Therefore, the dermal layer is actually modeled as an orthotropic material (Lai et al. 1999). With respect to the mechanical component of the model, the displacement of the dermal layer (u) is chosen as the primary model variable. Furthermore, we assume that it is appropriate to apply the infinitesimal strain theory in this study. Hence, we take where e is the Eulerian strain tensor and is the infinitesimal strain tensor. Finally, we select the following four constituents of the dermal layer as primary model variables: fibroblasts (N ), myofibroblasts (M), collagen bundles ( ρ ), and a generic signaling molecule (c). The general continuum hypothesis-based mathematical framework of Tranquillo and Murray (1992) is used as basis for the model. This framework consists of the following general set of conservation equations in local form: Equation (2a) is the mass conservation equation for the cell density/concentration of constituent i of the dermal layer and Eq. (2b) is the reduced conservation equation for the linear momentum of the dermal layer. It is assumed that the inertial forces that work on the dermal layer are negligible, and therefore, the conservation equation for the linear momentum of the dermal layer reduces to the above force balance equation. Within the above equations, z i represents the cell density/concentration of constituent i, v represents the displacement velocity of the dermal layer, J i represents the flux associated with constituent i per unit area, R i represents the (bio)chemical kinetics associated with constituent i, σ represents the Cauchy stress tensor associated with the dermal layer, and f represents the total body force working on the dermal layer. Given the chosen primary model variables, we have i ∈ {N , M, c, ρ jk } with j, k ∈ {1, 2, 3}. In the remainder of this text, we replace z N by N , z M by M, z c by c, and z ρ jk by ρ jk .

The cell populations
The functional forms for the movement of the (myo)fibroblasts and the functional forms for the biochemical kinetics associated with these cells are identical to functional forms used previously . For completeness we present these functional forms here as well. More details about the functional forms can be found in the cited article. The functional forms for the cell fluxes are where Parameter D F is the cell density-dependent random motility coefficient of the (myo)fibroblast, and χ F is the chemotactic coefficient. The functional forms for the biochemical kinetics associated with the (myo)fibroblasts are where r F is the cell division rate, r max F is the maximum factor with which the cell division rate can be enhanced due to the presence of the signaling molecule, a I c is the concentration of the signaling molecule that causes the half-maximum enhancement of the cell division rate, κ F F represents the reduction in the cell division rate due to crowding, q is a fixed constant, k F is the signaling molecule-dependent cell differentiation rate of fibroblasts into myofibroblasts, δ N is the apoptosis rate of fibroblasts, and δ M is the apoptosis rate of myofibroblasts.

The generic signaling molecule
The functional form for the dispersion of the generic signaling molecule and the functional forms for the release, the consumption, and the removal of the generic signaling molecule are also identical to functional forms used previously : The parameter D c represents the random diffusion coefficient of the generic signaling molecule, k c 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, a I I c is the concentration of the signaling molecule that causes the half-maximum net secretion rate of the signaling molecule, δ c is the proteolytic breakdown rate of the signaling molecules, 1/(1+a I I I c c) represents the inhibition of the removal of signaling molecules and collagen molecules (See Sect. 2.3) due to the presence of the signaling molecule, and tr ( ρ ) represents the concentration of the collagen molecules (See Sect. 2.3).

The collagen bundles
As has been mentioned in the introduction, we use a tensorial approach similar to the one proposed by Barocas and Tranquillo (1997) and Cumming et al. (2010) to represent the collagen bundles. The orientation of the collagen bundles and the concentration of the collagen molecules at location x and time t within the dermal layer are represented together by the symmetric tensor with (p(θ, ϕ)) T = [sin(ϕ) cos(θ ), sin(ϕ) sin(θ ), cos(ϕ)] the unit vector in the direction of the azimuthal angle θ and polar angle ϕ, and ρ(x, θ, ϕ, t) the concentration of the bundles at location x and time t with angle θ and angle ϕ.
Due to the symmetry of the above tensor, the tensor is orthogonal diagonalizable. Hence, the tensor can be decomposed as a sum of weighed outer products of orthonormal eigenvectors: where (λ i , q i ) are the eigenpairs of the tensor. These eigenpairs can be used to measure the degree of anisotropy of the dermal layer. That is, the larger the difference between the different (positive, real) eigenvalues, the more anisotropic the dermal layer. If all eigenvalues are equal, then the dermal layer is perfectly isotropic (See also Sect. 2.4). The eigenvector corresponding to the largest eigenvalue provides the dominant direction of the collagen bundles. The concentration of the collagen molecules at location x and time t can be recovered from the above tensor by either adding its eigenvalues or determining its trace. In the remainder of the text, we order the eigenvalues such that λ 3 ≥ λ 2 ≥ λ 1 . We assume that the secreted collagen molecules are attached to the ECM instantaneously. Hence, the flux associated with the ( j, k)th entry of the tensor ρ is Furthermore, we incorporate into the model the production of collagen molecules by both fibroblasts and myofibroblasts. Similar to the mechanism proposed by Olsen et al. (1998b), we assume that a portion of the collagen molecules are deposited and reoriented in the direction of movement of the (myo)fibroblasts. The remainder of the newly secreted collagen molecules are deposited by ratio in the direction of the present collagen bundles. The ratio of the amount of molecules that are deposited in the direction of movement of the cells to the amount of molecules that are deposited in the direction of the present collagen bundles is determined by the walking speed of the cells (i.e., the magnitude of the cell fluxes). We assume that the secretion rate of the molecules is enhanced in the presence of the signaling molecule. Finally, we use the same general functional form for the removal of the collagen molecules as was used previously . This removal takes place by ratio. Taken together, we obtain for the ( j, k)th entry of the tensor ρ where k ρ is the collagen molecule secretion rate, k max ρ is the maximum factor with which this rate can be enhanced due to the presence of the signaling molecule, a I V c N is the concentration of the signaling molecule that causes the half-maximum enhancement of the secretion rate, β ρ represents the sensitivity of (myo)fibroblasts to deposit and reorient (newly secreted) collagen molecules in the direction of cell movement, η is the ratio of myofibroblasts to fibroblasts in the maximum net secretion rate of collagen molecules, and δ ρ is the degradation rate of the collagen molecules. The constant γ is a small positive value that is added to the model to prevent the division by zero in the case of no flux of either cell type. In this study, we take γ = 10 −8 cells/(cm 2 day).

The force balance
Given that we model the dermal layer as a orthotropic material, we use the following general constitutive stress-strain relationship (i.e., σ = C ): with = ν 13 ν 21 ν 32 + ν 23 ν 12 ν 31 + ν 13 ν 31 and where (σ ) T = [σ 11 , σ 22 , σ 33 , σ 23 , σ 13 , σ 12 ], ν ·· are the Poisson ratios, E · are the Young's moduli, G ·· are the shear moduli, ( ) T = [ 11 , 22 , 33 , 2 23 , 2 13 , 2 12 ] (Lempriere 1968). Here, the axes of the material coordinate system coincide locally with the principal axes of the sample (i.e., the axes of the material coordinate system run parallel to the eigenvectors of the tensor ρ ). In order to make the bulk material properties dependent on the local geometrical arrangement of the collagen bundles, we propose the following definitions for these properties: where E and ν are constants. Using these definitions have two nice consequences: The symmetry of the elasticity tensor C is guaranteed and the positivity of the stored strain energy density in the dermal layer is guaranteed. The elasticity tensor is symmetric when the equality holds for j, k ∈ {1, 2, 3} and j = k (Lempriere 1968). The stored strain energy density is positive when the elasticity tensor is positive definite. This is the case when the inequalities hold (Lempriere 1968). Together, the dynamics and the initial conditions of the modeled constituents of the dermal layer imply λ i (x, t) > 0 for all x ∈ x , for all time t and i ∈ {1, 2, 3} (with x the domain of computation in Eulerian coordinates). Combined with the proposed definitions for the bulk material properties of the dermal layer and the values chosen for the constants E and ν, these are sufficient conditions to guarantee that the above equalities and inequalities hold; checking this is straightforward. Hence, the elasticity tensor is indeed symmetric, positive definite, and consequently, the stored strain energy density in the system is always positive. Notice furthermore that if the distribution of the collagen bundles is uniform (i.e., λ 1 = λ 2 = λ 3 ), then all Poisson ratios are equal. Likewise, all Young's moduli and shear moduli are equal. This implies that the elasticity tensor becomes equal to the elasticity tensor of an isotropic material (Lai et al. 1999). This is a nice property because this is exactly what you would expect given the uniformity of the distribution of the collagen bundles in the dermal layer.
The tensors from Eq. (14) need to be transformed so that they coincide with the global coordinate system (i.e., the Eulerian coordinate system). Due to the made assumptions with respect to the derivatives of the cell densities (See Sect. 2.5), the chosen initial conditions for the distribution of the collagen bundles [i.e., Eq. (33)] and the included dynamics for the production of collagen molecules [i.e., Eq. (13)], the first axis of the local material coordinate system always runs parallel to the first axis of the global coordinate system. This implies that the following transformations can be used to transform the tensors from Eq. (14): with and where σ = C , (σ ) T = [σ x x , σ yy , σ zz , σ yz , σ xz , σ xy ], ( ) T = [ x x , yy , zz , 2 yz , 2 xz , 2 xy ], m = cos(ϕ r ), n = sin(ϕ r ), and ϕ r is the angle of the clockwise rotation that, respectively, aligns the second and the third axis of the local material coordinate system (i.e., the second and third eigenvector of the tensor ρ ) with the second and the third axis of the global coordinate system. (Please note that the following holds: (T 1 (ϕ r )) −1 = T 1 (−ϕ r ). This equality simplifies the calculation of the elasticity tensor C from Eq. (23).) Finally, we incorporate into the model the generation of an isotropic stress by the myofibroblasts due to their pulling on the ECM. This 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). No other forces are incorporated into the model. Taken together, we obtain where ψ is 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.

The domain of computation
We assume that u = ∂v/∂x = ∂w/∂x = 0 holds within the modeled portion of the dermal layer, with the yz-plane running parallel to the surface of the skin and u = (u, v, w) T . Furthermore, 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 oriented parallel to the surface of the skin. Using Lagrangian coordinates (X = (X, Y, Z ) T ), the domain of computation X is described by

The initial conditions and the boundary conditions
The initial conditions give a description of the cell densities and the concentrations at the onset of the proliferative phase of the wound healing cascade. For the generation of the simulation results, we use two wound shapes: a circular shape and a square shape. For the construction of these shapes, we use the following general indicator function The values for the parameters s 1 and s 2 determine, respectively, the location and the steepness of the boundary of the wounded area. In order to construct the circular wounds, we use the following function w c (X) = I X , 16/π , 0.10 , In order to construct square wounds with the same surface area, we use the following function ×I (Y r , 2, 0.10) I (Z r , 2, 0.10) , with X r = R(θ r )X = (X r , Y r , Z r ) T and θ r = π/4 rad, where R(θ ) is the counterclockwise rotation matrix that rotates vectors by an angle θ about the X -axis. For i ∈ {c, s}, w i = 0 corresponds to completely wounded dermis and w i = 1 corresponds to unwounded dermis. Based on these general functions for the shapes of the wound, we take the following initial conditions for the cell densities and the concentration of the signaling molecules: where N and M are, respectively, the equilibrium cell density of the fibroblast population and the equilibrium cell density of the myofibroblast population of the unwounded dermis. The constant c w represents the maximum initial concentration of the signaling molecule in the wounded area. We assume that there are some fibroblasts and collagen bundles present in the wounded area; the parameter I w determines the minimum amount of fibroblasts and collagen bundles present initially in the wounded area.
For the initial distribution of the collagen bundles, we take the following: for j ∈ {X, Y, Z }, where e j is the unit vector that runs parallel to the jth coordinate axis. The value for the parameter r a determines which proportion of the collagen bundles are oriented in the direction perpendicular to the surface of the skin and which proportion of the collagen bundles are running parallel to the surface of the skin. In uninjured skin, the majority of the collagen bundles of the dermal layer run parallel to the surface of the skin, while only a small portion of the fibers are oriented out-of-plane (Holzapfel 2001; Annaidh et al. 2012). Therefore, we set r a to a relatively large value. The values for the parameter r b and the angle θ b together determine the overall distribution of the collagen bundles that run parallel to the surface of the skin. Over simulations, we vary the values for these latter two parameters in order to investigate the effect of such a variation on the contraction of wounds. The used (ranges of) values are depicted in Table 1 in the appendix. The parameter ρ is the equilibrium concentration of the collagen molecules in the unwounded dermis. With respect to the initial conditions for the displacement of the dermal layer, we observe the following. The initial cell density of the myofibroblasts is equal to zero everywhere in the domain of computation. Looking at Eq. (27), this implies f(x, 0) = 0. Therefore, See Fig. 1 for graphical representations of the initial conditions that have been used in this study. With respect to the boundary conditions for the modeled constituents of the dermal layer, we take the following Dirichlet boundary conditions where c is the equilibrium concentration of the signaling molecules in the unwounded dermis. Finally, we take the following Dirichlet boundary condition for the mechanical component of the model Table 1 in the appendix provides an overview of the dimensional (ranges of the) values for the parameters of the model. Most 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 other parameters ).

The applied numerical algorithm
In order to express the kernel of the algorithm and generate simulation results, we used MATLAB together with MATLAB's Parallel Computing Toolbox (The MathWorks Inc 2014). Furthermore, we interfaced the kernel consecutively with a slightly adapted version of the mesh generator developed by Persson and Strang (2004) for the generation of a base triangulation of the domain of computation, the element resolution refinement/recoarsement tool of the computational fluid dynamics (CFD) software package FEATFLOW2 for the adjustment of the resolution of the elements of the base triangulation (Turek 1998), and the scaling and permutation routine HSL_MC64 for the scaling and permutation of the linear systems (HSL 2013). The equations of the model were non-dimensionalized furthermore by applying the following non-dimensionalization: where L = 1 cm is the length scale of the model. The variables with the asterisks are the non-dimensionalized variables.
In order to find a solution for the system of equations from Eq. (2), we used the method of lines together with the standard fixed-point defect correction method (Van Kan et al. 2014). The two equations of the system were solved in a segregated way. That is, each time-step approximations of the solutions for the modeled constituents of the dermal layer were determined first, and subsequently, an approximation of the solution for the displacement of the dermal layer was determined. This scheme was iterated until certain standard convergence criteria were met. For the discretization of the system of equations, the first-order backward Euler timeintegration method was used together with a moving-grid finite-element method (Madzvamuse et al. 2003). On the discretized system of equations that describe the dynamics of the modeled constituents of the dermal layer, we applied a semi-implicit flux-corrected transport (FCT) limiter and a source term splitting procedure in order to enforce positivity of the approximations of the solutions for the constituents of the dermal layer Patankar 1980). Finally, we used an automatically adaptive time-stepping method to select the sizes of the individual time steps (Kavetski et al. 2002).
More details about certain aspects of the applied numerical algorithm can be found in one of our previous studies . More details about the procedure surrounding the adjustment of the resolution of the elements of the base triangulation can be found in the studies by Möller (2008); Möller and Kuzmin (2006).

Simulation results
In order to obtain some insight into the dynamics of the model, we present an overview of a simulation in Fig. 2. Figure 3 shows what the impact is of changing the initial distribution of the collagen bundles that run parallel to the surface of the skin, on the contraction of a circular wound. Due to the symmetry properties of a circular wound, it is not interesting to vary the value for the angle θ b . (In the case of a perfectly circular wound, changing the value for the angle θ b will merely rotate the solution counterclockwisely through an angle θ b .) Hence, solely the value for the ratio r b is changed over simulations. Note that changing the value for the ratio r b not only influences the degree of overall contraction of the wounded area (as measured by the evolution over time of the relative surface area of the wound compared to the surface area of the wound at day 0), but also the evolution over time  Table 1. The first 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 cells/cm 3 . The last two rows show the evolution over time of the concentration of, respectively, the signaling molecules and the collagen molecules. The color scales represent the concentrations, measured in g/cm 3 . Within the subfigures, the scale along both axes is in centimeters   Figure 4 shows what the impact is of changing the initial distribution of the collagen bundles that run parallel to the surface of the skin, on the contraction of a square wound. (Note that if r b = 1, then the solution is not dependent on the value for the angle θ b . Looking at Eq. (33), this becomes clear straight away.) Contrary to the case where the wound is circular, it is interesting to change both the value for the angle θ b and the value for the ratio r b in the case of a square wound. Due to the symmetry properties of a square wound, it is sufficient to vary the angle θ b between 0 rad and π/4 rad. (A solution for an angle larger than π/4 rad can be obtained from a solution for an angle between 0 rad and π/4 rad by applying a proper reflection and/or rotation on this latter solution.) Similar to the case of a circular wound, changing the value for the ratio r b influences both the degree of overall contraction of the wounded area (as measured by the evolution over time of the relative surface area of the wound compared to the surface area at day 0) and the evolution over time of the shape of the wounded area (as measured by the evolution over time of the ratio of the lengths of the line segments OC and O E, the ratio of the lengths of the line segments O D and O F, and the surface area of the wound relative to the surface area of the quadrilateral C EG H).
Furthermore, Fig. 4 shows that the evolution over time of the shape of the wounded area is also influenced by the orientation of the collagen bundles relative to the position of the wound in the case of a square wound. Interestingly, this relative orientation of the collagen bundles has hardly any influence on the evolution over time of the surface area of the wound relative to the surface area at day 0. Although not depicted in Fig. 4, we want to mention here that similar to the case of a circular wound, the evolution over time of the cell densities and the concentrations of the different modeled constituents is hardly influenced by changing the values for the ratio r b and the angle θ b .
Finally, Figs. 5 and 6 show the evolution over time of the geometrical arrangement of the collagen bundles during the healing of a circular wound. In order to demonstrate the effect of either including the deposition/reorientation of collagen molecules in the direction of cell movement or not, we set  Table  1 in the appendix. See Fig. 1b Table 1 for the generation of the simulation results in Fig. 6. (See the text next to Eq. (13) for a description of the parameter β ρ .) If β ρ = 0 (cm 2 day)/cells, then this will result ultimately in newly generated tissue with a collagen bundle distribution that is exactly equal to the collagen bundle distribution of the surrounding uninjured tissue. If β ρ = 10 −1 (cm 2 day)/cells, then ultimately the majority of the collagen molecules ends up permanently oriented toward the center of the wound and in the plane that runs parallel to the surface of the skin.

Discussion
We have presented a continuum hypothesis-based model for the simulation of the collagen bundle distribution-dependent contraction and retraction of healing dermal wounds that cover a large surface area. In order to simulate the con-traction and subsequent retraction of these wounds, some of the subprocesses that take place during the proliferative and remodeling phase of the wound healing cascade were incorporated into the model. Since wound contraction mainly takes place in the dermal layer of the skin, solely a portion of this layer is included explicitly into the model. This portion of dermal layer is modeled as a heterogeneous, orthotropic continuous solid with bulk mechanical properties that are locally dependent on both the local concentration and the local geometrical arrangement of the collagen bundles.
The following four constituents of the dermal layer were selected furthermore as primary model variables: fibroblasts, myofibroblasts, generic signaling molecules, and collagen bundles. The functional forms 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 the generic signaling molecule are identical to the functional forms used previously ).

Fig. 5
An overview of the evolution over time of the geometrical arrangement of the collagen bundles during the healing of a circular wound. For the generation of these simulation results, we took r b = 1, θ b = 0 rad and β ρ = 0 (cm 2 day)/cells. The values for the remaining parameters are equal to those depicted in Table 1. λ 1 , λ 2 , and λ 3 are the eigenvalues of the tensor ρ and θ q3 is the angle between the eigenvec-tor related to the largest eigenvalue λ 3 (where the third element of this eigenvector is larger than/equal to zero) and the positive horizontal axis. Within the subfigures on the bottom row, the transparency (i.e., α(x, t)) is set to either zero (opaque) or one (fully transparent) based on the fol- With respect to the way that the cell differentiation of fibroblasts into myofibroblasts has been incorporated into the presented model, we want to place some remarks here. Similar to Olsen et al. (1995), we assume that the rate of cell differentiation is dependent on the concentration of the signaling molecule with no differentiation taking place in the absence of the signaling molecule. Like others such as Murphy et al. (2012) and Valero et al. (2014), we are aware of the fact that this cell differentiation can only take place under conditions of sufficient mechanical stiffness. However, as is also clearly pointed out by Van de Water et al. (2013), it is unclear at present what the actual stiffness is that is perceived by fibroblasts. Furthermore, recent experimental studies have shown that the differentiation of fibroblasts into myofibroblasts is also critically dependent on the presence of particular isoforms of fibronectin (Van de Water et al. 2013). Taken together, these issues imply that the incorporation of the cell differentiation mechanism into a mathematical model in a realistic way is basically impossible at present. Hence, while being aware of the fact that the cell differentiation process is a very complex process in reality, we decided to keep things relatively simple in this modeling study.
For the representation of the collagen bundles, we used a tensorial approach similar to the one proposed by Barocas and Tranquillo (1997) and Cumming et al. (2010). Com-pared to the approaches developed by Dallon et al. (1999Dallon et al. ( , 2001Dallon et al. ( , 2000; McDougall et al. (2006) and Olsen et al. (1999Olsen et al. ( , 1998b, this approach has some advantages. Dallon et al. use a vector-based representation for the representation of collagen bundles in their models. Using this vector-based representation has two major disadvantages (Cumming et al. 2010). Firstly, it does not provide any information about the degree of isotropy of the collagen bundles at individual material points within the dermal layer. Secondly, due to the use of a vector representation, collagen bundles are treated basically as unidirectional entities. Given that fiber alignment is bidirectional, this may lead to duality when one wishes to determine the degree of alignment of collagen bundles. These two disadvantages are overcome by using a tensorial approach. Furthermore, Olsen et al. model the alignment of collagen bundles in their models by considering two orthogonal configurations and assuming that the transition between these two configurations is a dynamic and reversible process. The major disadvantage of this latter approach is that this representation cannot provide an accurate representation of the orientation of the collagen bundles when this orientation is continuously distributed. The tensorial approach does not suffer from this drawback.
With respect to the dynamic change of the geometrical arrangement of the collagen bundles it is assumed that a por- Fig. 6 An overview of the evolution over time of the geometrical arrangement of the collagen bundles during the healing of a circular wound. For the generation of these simulation results, we took r b = 1, θ b = 0 rad and β ρ = 10 −1 (cm 2 day)/cells. The values for the remaining parameters are equal to those depicted in Table 1. λ 1 , λ 2 , and λ 3 are the eigenvalues of the tensor ρ and θ q3 is the angle between the eigenvector related to the largest eigenvalue λ 3 (where the third element of this eigenvector is larger than/equal to zero) and the positive horizontal axis. Within the subfigures on the bottom row the transparency (i.e., α(x, t)) is set to either zero (opaque) or one (fully transparent) based on the following rule: tion of the collagen molecules are deposited and reoriented in the direction of movement of the (myo)fibroblasts (Olsen et al. 1998b). The remainder of the newly secreted collagen molecules are deposited by ratio in the direction of the present collagen bundles. The ratio of the amount of molecules that are deposited in the direction of movement of the cells to the amount of molecules that are deposited in the direction of the present collagen bundles is determined by the walking speed of the cells.
With the developed model, it is possible to simulate some general qualitative features of the dermal wound healing response (Baum and Arpey 2005;Enoch and Leaper 2007;Li et al. 2007;Monaco and Lawrence 2003). The restoration of the presence of fibroblasts in the wounded area can be simulated. Furthermore, the initial expansion and subsequent reduction of the myofibroblast population in the wounded area during the execution of the wound healing response can be simulated, and it is possible to simulate the restoration of a collagen-rich ECM in the recovering wounded area. In accordance with experimental observations (Zuijlen et al. 2003;Welch et al. 1990), it is also possible with this model to simulate the permanent increase in the proportion of the collagen bundles that runs parallel to the surface of the skin as a consequence of the execution of the wound healing process.
Finally, the contraction and the subsequent retraction of the wounded area can also be simulated.
With respect to the simulation of the collagen bundle distribution-dependent contraction and subsequent retraction of healing dermal wounds, the following can be observed. Figures 3 and 4 show clearly the impact of changing the initial distribution of the collagen bundles that run parallel to the surface of the skin. The distribution of the collagen bundles influences the evolution over time of both the shape of the wounded area and the degree of overall contraction of the wounded area. Interestingly, Fig. 3 shows that these effects are solely a consequence of alterations in the initial overall distribution of the collagen bundles, and not a consequence of alterations in the evolution over time of the different cell densities and concentrations of the modeled constituents. Furthermore, it is very interesting to observe in Fig. 4 that the evolution over time of the shape of the wound is also influenced by the orientation of the collagen bundles relative to the position of the wound, while this relative orientation does not influence the evolution over time of the relative surface area of the wound compared to the surface area of the wound at day 0. Figures 5 and 6 show clearly the effect of either including into the model the dependence of the geometrical arrange-ment of collagen bundles on the movement of cells or not. If it is included, then an increased portion of the collagen bundles ultimately ends up permanently in the plane running parallel to the surface of the skin, and the majority of these bundles are oriented toward the center of the wound. If the deposition/reorientation of collagen molecules in the direction of cell movement is not included, then all the newly secreted collagen molecules are deposited by ratio in the direction of the present collagen bundles. As demonstrated in Fig. 5, this will result ultimately in newly generated tissue with a collagen bundle distribution that is exactly equal to the collagen bundle distribution of the surrounding uninjured tissue.
As has been mentioned before, the fact that an increased portion of the collagen bundles ultimately ends up permanently in the plane running parallel to the surface of the skin and the fact that the majority of these bundles are oriented toward the center of the wound when the dependence of the direction of deposition/reorientation of collagen molecules on the movement of cells is included into the model, are in accordance with experimental observations (Zuijlen et al. 2003;Welch et al. 1990). This is an interesting observation. To the best of our knowledge, it is unknown at present which wound healing mechanisms cause these experimental observations. In the presented model, the dynamics related to the geometrical arrangement of the collagen bundles are dependent on the speed and the direction of movement of the (myo)fibroblasts. This direction and speed of movement are dependent subsequently on the gradient of the concentration of the signaling molecule and the evolution over time of the distribution of the cell densities. Taken together, the results obtained with the presented model suggest that the geometrical arrangement of collagen bundles in scar tissue might be altered by changing the gradient of the concentration of the signaling molecule and/or the evolution over time of the distribution of the cell densities and therewith the direction and speed of movement of the (myo)fibroblasts, during wound healing. Given that the geometrical arrangement of the collagen bundles has a huge impact on the response of dermal tissues to mechanical forces (Jor et al. 2011), this is an interesting suggestion because we think that its offers new ideas for a better treatment of deep dermal wounds that results in scar tissue that is more akin to the original tissue. We can imagine that it must be possible to investigate this suggestion further by means of an experimental study with a tissue equivalent.
Our ultimate goal is to construct mathematical models with which the structural properties of healing wounds can be predicted accurately. Since collagen bundle distributiondependent contraction and subsequent retraction are important components of the wound healing response in the case of deep dermal wounds that cover a large surface area, we consider the development of the presented model as a step toward the construction of such predictive models.
Obviously, the presented model can be expanded in different ways. One relatively simple addition to the model would be, for instance, to make the direction of cell movement also dependent on the orientation of the collagen bundles. Previous experimental studies have demonstrated that fibroblasts align with collagen bundles, which subsequently influences the direction of movement of these cells (Guido and Tranquillo 1993). In order to incorporate this effect into the presented model, we could, for example, replace Eqs. (3) and (4) with where c = ρ tr ( ρ ) .
In addition, Verhaegen et al. (2012) have demonstrated that the stretching of both healthy skin and scar tissue induces permanent adaptation of the orientation of the collagen bundles in these tissues. Therefore, it might be interesting to incorporate into the presented model the reorientation of collagen bundles due to forces. This could be accomplished perhaps by absorbing into the presented model some of the ideas that were used for the development of the model by Olsen et al. (1999). Furthermore, it might also be interesting to add to the presented model a morphoelastic component. Adding this component to the model will make it possible to simulate the often present permanent deformation of recovering dermal tissues. The incorporation of such an effect into the current model could be accomplished, for instance, by using the morphoelastic framework developed by Hall (2009). Finally, it is probably also very interesting to investigate in a three-dimensional portion of dermal layer what would happen to the shapes of wounds and the geometrical arrangement of collagen bundles during healing when the assumptions made in Sect. 2.5 are removed.