Computer Methods in Biomechanics and Biomedical Engineering

ISSN: 1025-5842 (Print) 1476-8259 (Online) Journal homepage: http://www.tandfonline.com/loi/gcmb20

Computational modeling of chemo-bio-mechanical coupling: a systems-biology approach toward wound healing A. Buganza Tepole & E. Kuhl To cite this article: A. Buganza Tepole & E. Kuhl (2014): Computational modeling of chemo-biomechanical coupling: a systems-biology approach toward wound healing, Computer Methods in Biomechanics and Biomedical Engineering, DOI: 10.1080/10255842.2014.980821 To link to this article: http://dx.doi.org/10.1080/10255842.2014.980821

Published online: 24 Nov 2014.

Submit your article to this journal

Article views: 97

View related articles

View Crossmark data

Full Terms & Conditions of access and use can be found at http://www.tandfonline.com/action/journalInformation?journalCode=gcmb20 Download by: [University of Otago]

Date: 01 November 2015, At: 20:21

Computer Methods in Biomechanics and Biomedical Engineering, 2014 http://dx.doi.org/10.1080/10255842.2014.980821

Computational modeling of chemo-bio-mechanical coupling: a systems-biology approach toward wound healing A. Buganza Tepolea* and E. Kuhla,b a

Department of Mechanical Engineering, Stanford University, Stanford, CA 94305, USA; bDepartment of Bioengineering, Stanford University, Stanford, CA 94305, USA

Downloaded by [University of Otago] at 20:21 01 November 2015

(Received 3 July 2014; accepted 22 October 2014) Wound healing is a synchronized cascade of chemical, biological, and mechanical phenomena, which act in concert to restore the damaged tissue. An imbalance between these events can induce painful scarring. Despite intense efforts to decipher the mechanisms of wound healing, the role of mechanics remains poorly understood. Here, we establish a computational systems biology model to identify the chemical, biological, and mechanical mechanisms of scar formation. First, we introduce the generic problem of coupled chemo-bio-mechanics. Then, we introduce the model problem of wound healing in terms of a particular chemical signal, inflammation, a particular biological cell type, fibroblasts, and a particular mechanical model, isotropic hyperelasticity. We explore the cross-talk between chemical, biological, and mechanical signals and show that all three fields have a significant impact on scar formation. Our model is the first step toward rigorous multiscale, multifield modeling in wound healing. Our formulation has the potential to improve effective wound management and optimize treatment on an individualized patient-specific basis. Keywords: wound healing; systems biology; chemo-bio-mechanics; finite elements; multiscale

1.

Motivation

Effective wound management is a quotidian concern in clinical practice. Abnormal wound healing can initiate hypertrophic scars associated with serious complications from deteriorated skin characteristics to psychological trauma (Bayat et al. 2003). The health-care cost related to wound treatment is jolting; wounds are common to many clinical procedures and span all patient demographics (Brown et al. 2008). Fostering a healthy tissue response is a non-trivial task. The process of wound healing is a complex sequence of interrelated events that involve mechanical cues, coordinated cell behavior, and the interaction of numerous chemical signals (Diegelmann and Evans 2004). In such a scenario, planning effective healing on a patientspecific basis becomes almost impossible. Computational systems biology has found a niche to enrich our understanding of this complex problem (Kitano 2002). However, despite intense efforts to characterize the healing process with mathematical models, simulation of wound healing in arbitrary three-dimensional geometries remains an open problem. Disrupting the integrity of skin triggers a cascade of events that are common to all inflammationbased systems in the human body (Vodovotz et al. 2008). In addition, during dermal wound healing, specialized processes take place to restitute the particular functional requirements of dermal tissue (Gurtner et al. 2008). Perhaps the most distinct feature of this system is the interaction of different key players across scales, both in space and time.

*Corresponding author. Email: [email protected] q 2014 Taylor & Francis

During the past decades, scientists have successfully identified and characterized the individual aspects of this network, but a holistic understanding of the healing process as a whole remains obscure (Vodovotz 2010).

1.1. Wound healing across the spatial scales The spatial scales of interest for the healing system range from the order of micrometers, to millimeters, centimeters, and decimeters (Buganza Tepole and Kuhl 2013). Figure 1 illustrates the multi-scale nature of the healing process with four interacting spatial scales (Qutub et al. 2009): the system level, the organ level, the tissue level, and the cell level (Hunter and Borg 2003). On the cell level, the smallest spatial scale of the order of micrometers, single cells are the individual actors, which directly affect the healing process (Southern et al. 2008). In the damaged dermal tissue and its surroundings, the following cell types are present: two types of leukocytes, neutrophils and macrophages, dispose pathogens and debris and establish gradients of growth factors; endothelial cells generate a new vasculature; keratinocytes divide and migrate across the epidermis to produce a new protective outermost layer; and fibroblasts deposit collagen and generate active stresses to contract the wound (Olutoye et al. 2005; Velnar et al. 2009). On the tissue level, the next larger scale of the order of millimeters, the actions of the individual cells are smeared

2

A. Buganza Tepole and E. Kuhl System level [10−1m]

Organ level [10−2m]

Tissue level [10−3m]

Cell level [10−5m]

Downloaded by [University of Otago] at 20:21 01 November 2015

Figure 1. Wound healing across the spatial scales. The chemo-bio-mechanical problem of wound healing spans from the cellular level via the tissue level and organ level to the system level bridging four orders of magnitude in space.

out by meso-scale patterns, which emanate from the collective cell behavior. This collective response exhibits distinctive characteristics, which cannot be extrapolated directly from single cell actions. Populations of keratinocyte generate a well-organized traveling wave inwards, and populations of endothelial cells keep strong cell –cell interactions to create fractal-like vascular networks (Pettet et al. 1996; Maini et al. 2004). From this scale upward, both the mechanical response and the reaction – diffusion response of the chemical species can be characterized through field variables using a continuum approach (Lanir 1983; Flynn et al. 2011). On the organ level, the scale of the order of centimeters, we can explore the interaction of the mechanical properties, different cell populations, and reaction – diffusion systems of chemical concentrations. The organ level provides a holistic approach to study the role of the individual key players of wound healing, and allows us to explore tissue function in health and disease. On the system level, the scale of the order of decimeters, we can study the entire system created by the interplay of different organs. Ideally, system level models

Hemostasis [min]

Inflammation [hours]

are generated on a patient-specific basis from individual clinical imaging data (Zo¨llner et al. 2012).

1.2.

Wound healing across the temporal scales

The temporal scales of interest for the healing system range from the order of minutes, to hours, days, and weeks (Buganza Tepole and Kuhl 2013). Figure 2 illustrates the multi-scale nature of the healing process with four overlapping temporal scales: hemostasis, inflammation, proliferation, and remodeling (Martin 1997). Immediately after the injury occurs, healing is critical to restitute the barrier function of skin. Unfortunately, the initially generated temporary scaffold has only poor mechanical characteristics. Accordingly, subsequent stages of the healing process gradually reconstruct the tissue to ultimately restore the constitution of the uninjured skin (Mutsaers et al. 1997). The entire healing process can last for weeks or even months. During hemostasis, within the order of minutes, the injured region fills with blood, which quickly coagulates. This results in the formation of an emergency scaffold of fibrin fibers. The only cells present in the clot are platelets, Proliferation [days]

Remodeling [weeks]

Figure 2. Wound healing across the temporal scales. The chemo-bio-mechanical problem of wound healing spans from the homeostatic phase via the inflammatory phase and proliferative phase to the remodeling phase bridging four orders of magnitude in time.

Downloaded by [University of Otago] at 20:21 01 November 2015

Computer Methods in Biomechanics and Biomedical Engineering responsible for coagulation and the release of growth factors. At the end of this stage, degranulation of the platelets floods the injured site with chemicals to attract leukocytes. During inflammation, within the order of hours, the first population of leukocytes, neutrophils, arrives at the wound site. Neutrophils remove pathogens and dispose of tissue debris from the injury. Shortly after, a second population of leukocytes, macrophages, migrates into the wound and continues the cleaning process. In addition, they establish gradients of various chemical signals to attract other cell populations (Tranquillo 1987). After one or two days, the inflammatory phase smoothly blends into the proliferative phase. During proliferation, within the order of days, the chemical signaling established by the macrophages attracts specialized cell populations that reconstruct skin. Endothelial cells generate a new vasculature that provides nutrients to the other cell populations (Herbert and Stainier 2011). Keratinocytes reconstruct the outermost protective layer, the epidermis, in a process called re-epithelialization (Simpson et al. 2011). Fibroblasts replace the temporary fibrin scaffold with a collagenous matrix that restitutes the desired mechanical properties of the healed tissue (Chiquet et al. 2009). Although the proliferation phase creates a somewhat functional tissue, the mechanical properties of the newly reconstructed skin are not nearly identical to healthy, uninjured skin: the newly generated material is stiff scar tissue, which is partly provisional and will be replaced during the final remodeling phase (Verhaegen et al. 2009). During remodeling, within the order of weeks, fibroblasts slowly tear down and deposit collagen until the matrix approaches the structure of healthy tissue. The remodeling phase can continue for months or even years. 1.3. Modeling wound healing Wound healing has been studied assiduously both experimentally and theoretically. Recent developments in computational systems biology suggest that we cannot gain a complete understanding of wound healing from studying isolated spatial or temporal scales alone (Sun et al. 2009). Rather, trends in modeling seem to converge toward assembling the individual building blocks for a holistic model that, once calibrated, can provide new insight into the baseline system (Dallon et al. 2001). Systematic perturbations of this system allow us to probe different healing scenarios to ultimately link computational tools with personalized models (Xue et al. 2009). The first mathematical model of wound healing was introduced in the early 1990s. Its initial goal was to simulate the traveling wave front of growing cell populations at the edge of a wound (Sherratt and Murray 1990). Since then, mathematical models have gained in

3

complexity and have gradually incorporated the different components that interact in synchrony to heal the damaged tissue (Xue et al. 2009). Recent models can be categorized according to two criteria: the aspect of healing they seek to analyze in detail and the simulation framework employed for the analysis. Three aspects of healing are particularly relevant: re-epithelialization and cell migration (Ben Amar and Wu 2014), angiogenesis (Vermolen and Javierre 2012), and mechanical aspects of wound healing such as collagen deposition and wound contraction (GarzonAlvarado et al. 2012). Four modeling strategies are prevalent: one-dimensional and axisymmetric continuum models; two-dimensional continuum models; two- and three-dimensional discrete models; and two-dimensional hybrid discrete/continuum models. A recent review article highlights the state of the art in systems biology approaches toward wound healing (Buganza Tepole and Kuhl 2013). Among the different variables that influence the outcome of healing, the importance of mechanical cues has recently been identified with more clarity (Agha et al. 2011). Fibroblasts have the capability to sense mechanical signals, to translate them into specific action such as active contraction, and to release chemical substances (Wong et al. 2011). We now know that increased stress in the wound site alters fibroblast phenotype by reducing their apoptotic rate and inducing the release of pro-inflammatory signals (Aarabi et al. 2007; Paterno et al. 2011). In turn, when the inflammation phase is prolonged, fibroblasts divide and migrate into the wound at higher rates, which results in an increased collagen deposition. The ultimate consequence is a poorly structured dermal tissue with thick collagen bundles instead of the smooth, inter-woven collagenous network found in healthy tissue. Visually, the result of such pathological reaction is very clear to the human eye, which can readily recognize hypertrophic scars (Wong et al. 2012). However, despite this obvious evidence, only few models have incorporated a detailed mechanical description of the wound environment. While some models have addressed collagen deposition and active wound contraction, their application is limited by the underlying simplified constitutive models for the tissue structure (Tranquillo and Murray 2007; Javierre et al. 2009). While such simplifications are adequate for baseline studies and first prototype simulations (Kuhl and Steinmann 2004), in the current state, these models are unable to bridge the gap toward arbitrary geometries, large deformations, and complex stress distributions that arise in more realistic settings. Motivated by the need for a computational framework that incorporates the state-of-the-art development in wound healing, here we present a novel finite element formulation for the chemo-bio-mechanical problem of wound healing in arbitrary geometries. The manuscript is structured as follows. In Section 2, we introduce the

4

A. Buganza Tepole and E. Kuhl

generic continuum framework of wound healing. In Section 3, we specify a particular type of the model parameterized in terms of a single chemical signal, a single biological cell density, and the mechanical deformation. In Section 4, we derive the discrete formulation of the particular model problem. In Section 5, we present sensitivity studies and selected examples to showcase the features of the model. Finally, in Section 6, we provide a discussion and a brief outlook.

chemical concentrations c and all cell populations r. They contain the information about how chemical substances are produced and degraded through chemical reactions with other chemical substances and by the different biological cells. In homeostasis, in the absence of chemical gradients 7c ¼ 0, the chemical production and degradation balance each other, f cp ¼ f cd c. 2.2.

Downloaded by [University of Otago] at 20:21 01 November 2015

2.

Chemo-bio-mechanical problem

We begin by introducing the generic equations that govern the dynamics of inflammation-based systems. In general, the underlying chemo-bio-mechanical problem can be characterized through three spatially and temporally interacting building blocks: chemical fields including substances such as growth factors and inflammation signals, here summarized in the vector cðX; tÞ ¼ ½c1 ðX; tÞ; c2 ðX; tÞ; :::; cnc ðX; tÞ
t ; biological fields including cell populations, here summarized in the vector rðX; tÞ ¼ ½r1 ðX; tÞ; r2 ðX; tÞ; :::; rnr ðX; tÞ
t ; and mechanical fields including the deformation w ðX; tÞ, which can be locally supplemented by microstructural internal variables such as microstructural directions nðX; tÞ ¼ ½n 1 ðX; tÞ; n 2 ðX; tÞ; :::; n nn ðX; tÞ
t or microstructural concentrations wðX; tÞ ¼ ½w1 ðX; tÞ; w2 ðX; tÞ; :::; wnw ðX; tÞ
t . In the following, we characterize the evolution equations of these sets of variables in a continuum setting.

2.1. Chemical problem: chemical concentrations Chemically, the evolution of the set of chemical concentrations c is balanced by the chemical flux q c and the chemical source f c , c_ ¼ 2divq c ð7cÞ þ f c ðc; rÞ;

ð1Þ

where f_+} ¼ df+}=dt denotes the material time derivative and 7f+} and divf+} denote the spatial gradient and divergence. The chemical flux q c is typically modeled as a linear function of the gradient of the chemical concentration 7c to indicate that the chemical signal can diffuse freely in the domain of interest, q c ¼ 2D cc 7c;

ð2Þ

where D cc denotes the chemical diffusion tensor. The chemical source f c typically consists of a production term f cp and a degradation term f cd , whereby the degradation typically scales linearly with the concentration c (Sherratt and Murray 1990), f c ðc; rÞ ¼ f cp ðc; rÞ 2 f cd ðc; rÞc:

ð3Þ

In general, f cp and degradation f cd can be functions of all

Biological problem: biological cell densities

Biologically, the evolution of the set of cell densities r is balanced by the biological flux q r and the biological source f r ,

r_ ¼ 2divq r ðc; 7c; r; 7r; 7w Þ þ f r ðc; r; 7w Þ:

ð4Þ

The biological flux q r typically consists of three contributions, q r ¼ 2D rr 7r 2 D rc ðc; r; 7wÞ7c 2 D rw : 7w:

ð5Þ

The first contribution, 2D rr 7r, describes the free diffusion of cells along cell density gradients 7r. It mimics the continuum representation of random walk and contact inhibition at the cellular level represented through the biological diffusion tensor D rr . The second contribution, 2D rc 7c, characterizes the phenomenon of chemotaxis. It is associated with the directed diffusion along chemical concentration gradients 7c represented through the chemotactic diffusion tensor D rc , which can either be constant or depend on chemical concentrations c, cell densities r, and deformation w . The third contribution, 2D rw : 7w , represents the phenomenon of mechanotaxis. It reflects the directed diffusion along mechanical cues 7w represented through the mechanotactic diffusion tensor D rw . The biological source consists of a mitotic contribution f rm and an apoptotic contribution f ra , which typically scales linearly with the cell density r, f r ðc; r; 7w Þ ¼ f rm ðc; r; 7w Þ 2 f ra ðc; r; 7w Þr:

ð6Þ

The mitotic and apoptotic terms f rm and f ra can be functions of all chemical concentrations c, of all cell populations r, and of mechanical cues 7w . The latter dependency mimics the effects of mechanotransduction, the impact of mechanical cues on biological phenomena. In homeostasis, in the absence of biological gradients 7r ¼ 0, the mitotic and apoptotic rates balance each other, f rm ¼ f ra r. 2.3. Mechanical problem: mechanical deformation Mechanically, we assume that the mechanical problem is quasi-static and balances the mechanical flux s with the

Computer Methods in Biomechanics and Biomedical Engineering mechanical source f w , 0 ¼ divsð7w ; n; wÞ þ f w :

ð7Þ

The mechanical flux s, the Cauchy stress, can be additively decomposed into passive and active contributions,

s ¼ s pas ð7w ; n; wÞ þ s act ð7w ; n; wÞ:

ð8Þ

The active stress accounts for tissue contraction by cells such as fibroblasts (Javierre et al. 2009). The passive stress typically consists of two contributions,

Downloaded by [University of Otago] at 20:21 01 November 2015

s pas ¼ s mat ð7w Þ þ ws fib ð7w ; nÞ:

ð9Þ

The first contribution s mat describes the isotropic waterbased matrix as a function of the deformation gradient 7w . The second contribution s fib describes the anisotropic response of fibrous constituents such as elastin, collagen, or smooth muscle as a function of the deformation gradient 7w and preferred microstructural directions n, scaled by the fiber content w. The mechanical source f w , the external mechanical force such as gravity, is typically negligible in the context of inflammation-based systems, f w ¼ 0:

ð10Þ

Biological cells continuously interact with and remodel the tissue in their immediate environment to establish a welldefined microstructural arrangement in healthy tissue. After an injury, this microstructure of the healthy skin disappears. Local remodeling by cells becomes the crucial connecting point between the chemical, biological, and mechanical fields (McDougall et al. 2006). We typically model this coupling through the internal variables n and w, which evolve in response to cell population dynamics r. In the most general setting, we allow the microstructural directions n to gradually reorient according to a set of local evolution equations (Garikipati et al. 2006), e.g., driven by chemical gradients 7c (Bowes et al. 1999), by biological gradients 7r, by mechanical gradients 7w , and by the current microstructural directions n, n_ ¼ f n ð7c; 7r; 7w ; nÞ:

ð11Þ

Similarly, the local fiber content w can evolve in time, e.g., driven by chemical concentrations c, by biological cell densities r, by mechanical gradients 7w , and by the current fiber content w, _ ¼ f w ðc; r; 7w ; wÞ: w

ð12Þ

Even though the microstructural direction n and the microstructural fiber content w are parameterized in terms of inhomogeneous fields, their evolution equations are strictly local as they do not contain any gradient or

5

divergence terms. This suggests to treat the microstructural information n and w as a set of internal variables (Menzel 2007; Himpel et al. 2008). In summary, we represent the chemo-bio-mechanical problem through a system of three sets of partial differential equations for the chemical concentrations c, the biological cell densities r, and the mechanical deformation w , locally supplemented by two sets of ordinary differential equations for the microstructural directions n and the microstructural fiber content w. In the following section, we specify these generic equations to explore a particular model problem of wound healing.

3.

Model problem of wound healing

In this section, we illustrate the features of the proposed generic framework in terms of a simple model problem of wound healing restricting attention to a few key players. We represent the chemical problem through the concentration of the inflammatory signal cðX; tÞ, the biological problem through the fibroblast density rðX; tÞ, and the mechanical problem through the deformation w ðX; tÞ supplemented by the collagen content wðX; tÞ as local internal variable. We assume that the collagen fiber orientation nðX; tÞ remains constant throughout the healing process.

3.1.

Chemical problem: inflammatory signal

Chemically, we characterize the response through the inflammatory signal c, which represents the initial recruitment of macrophages and their contribution to generate growth factor attractors for fibroblasts. In reality, the cascade of chemical signaling of inflammation is much more complex and includes several cytokines and other cell types such as endothelial cells and neutrophils (Sherratt and Murray 1990; Werner and Grose 2003). Nonetheless, previous mathematical models have shown good qualitative behavior considering only macrophages and fibroblasts (Cumming et al. 2009). We follow that approach here and synthesize the inflammatory signal into a single field variable. According to the generic chemical balance law (1), we balance its rate of change with the chemical flux q c and the chemical source f c , c_ ¼ 2div q c þ f c :

ð13Þ

For the chemical flux q c , we assume a linear isotropic function of the gradient of the chemical concentration 7c to indicate that the chemical signal can diffuse freely and isotropically in the domain of interest (Chary and Jain 1989), q c ¼ 2D cc 7c:

ð14Þ

where D cc is the isotropic chemical diffusion coefficient. For the chemical source, we assume that the inflammatory

6

A. Buganza Tepole and E. Kuhl

signal has no production component and displays a linear degradation, f c ¼ 2kc;

ð15Þ

where k is the chemical degradation rate.

3.2.

c ¼ c mat ðJ; I 1 Þ þ wc fib ðI 1 ; I 4 Þ;

Biological problem: fibroblasts

Biologially, we characterize the response through the fibroblast density r. According to the generic biological balance law (4), we balance its rate of change with the biological flux q r and the biological source f r ,

r_ ¼ 2div q r þ f r : Downloaded by [University of Otago] at 20:21 01 November 2015

n 0 (Buganza Tepole et al. 2012). For this simple model problem of wound healing, we do not consider the active stress exerted by the fibroblast cell population. We characterize only its passive constitutive response through a compressible, transversely isotropic, hyperelastic free energy function,

ð16Þ

For the biological flux, we assume that fibroblasts are motile cells, which diffuse freely along their own gradients 7r perturbed by a biased diffusion toward the gradient of the inflammatory signal 7c, q r ¼ 2D rr 7r 2 ar7c;

ð17Þ

where D rr and D rc ¼ ar denote the isotropic biological and chemotactic diffusion coefficients. For the biological source, we make the following ansatz in terms of the fibroblast density r and the intensity of the inflammatory signal c, f r ¼ k1 ½r0 2 r
þ k2 cr;

ð18Þ

where r0 is the homeostatic fibroblast concentration, k1 is the physiological mitotic and apoptotic rate, and k2 is the mitotic rate induced by the inflammatory signal c. Under healthy conditions, fibroblast mitosis and apoptosis balance one another to ensure a stable fibroblast population r0 . However, in the presence of inflammatory signals, the mitotic rate increases and creates an imbalance with respect to the steady-state r0 to increase the fibroblast density.

which consists of an isotropic part c mat for the noncollagenous isotropic matrix and an anisotropic part c fib for the collagen network weighted by the collagen content w. Here, we have introduced three kinematic invariants, the Jacobian J ¼ det F for the volumetric response, the first invariant I 1 ¼ b : I for the isotropic response, and the fourth invariant I 4 ¼ ½n^n
: I for the anisotropic response, where n ¼ Fn 0 is the preferred collagen fiber orientation in the deformed configuration. We model the matrix material as standard isotropic, compressible NeoHooke-type parameterized in terms of the Lame´ constants l and m and the collagen fibers as Holzapfel-type (Holzapfel et al. 2000), parameterized in terms of the collagen stiffness c1 , the nonlinearity c2 , and the fiber dispersion k, 1 1 c mat ¼ m½I 1 2 3
2 m lnJ þ l ln2 J 2 2 ð21Þ 1 c1 fib c ¼ ½expðc2 ½kI 1 þ ½1 2 3k
I 4 2 1
2 Þ 2 1
: 2 c2 The additive decomposition of the strain energy function translates into the additive decomposition of the Cauchy stress according to the generic ansatz (9),

s ¼ F

Mechanical problem: deformation

Mechanically, we characterize the response through the deformation w , from which we derive the deformation gradient F ¼ 7w and the left Cauchy – Green deformation tensor b ¼ FF t as key kinematic variables. According to the mechanical balance law (7), we balance the mechanical flux s characterizing the Cauchy stress and the mechanical source f w characterizing the external mechanical forces, 0 ¼ div s þ f : w

ð19Þ

Skin has a well-organized microstructure with an isotropic water-based matrix that serves as a scaffold for the anisotropic collagen network with a preferred orientation

2 ›c F t ¼ s mat þ ws fib ; J ›C

ð22Þ

with the following matrix and fiber contributions, 2 J 2 ¼ F J

s mat ¼ F 3.3.

ð20Þ

s

fib

›c mat F t ¼ m ½b 2 I
þ l ln JI ›C ›c fib F t ¼ 2c1 b þ 2c4 n^n: ›C

ð23Þ

where c1 and c4 denote the first derivatives of the fiber energy with respect to the first and fourth invariants, ci ¼ ›c fib =›I i ,

c1 ¼ c0 þ c1 k½kI 1 þ ½1 2 3k
I 4 2 1
expðc2 ½kI 1 þ ½1 2 3k
I 4 2 1
2 Þ

c4 ¼ c1 ½1 2 3k
½kI 1 þ ½1 2 3k
I 4 2 1


ð24Þ

expðc2 ½kI 1 þ ½1 2 3k
I 4 2 1
2 Þ: For the mechanical source, we assume that external forces such as gravity do not play a major role during wound

Computer Methods in Biomechanics and Biomedical Engineering healing and can therefore be neglected,

Downloaded by [University of Otago] at 20:21 01 November 2015

f w ¼ 0:

ð25Þ

It remains to characterize the evolution of the collagen content w and its reference orientation n, which we treat locally as internal variables. For the sake of simplicity, we assume that the fiber orientation remains constant throughout the healing process. This implies that new collagen is deposited with the same orientation as the surrounding, healthy tissue. We assume that healthy skin possess a homeostatic collagen content w0. After an injury, the collagen content decreases dramatically in the affected region. It is the task of the fibroblasts to deposit new collagen to restore skin’s mechanical properties. We therefore model the evolution of the collagen content w _ to depend on the chemical signal c, the fibroblast density r, and the current collagen content w itself, specifically,
 rg a cg w_ ¼ f w with f w ¼ 1 2 w þ : ð26Þ 1 þ w2 1þc where g denotes the physiological collagen deposition rate and a denotes the increase in collagen synthesis in response to inflammation c.

4. Computational modeling of wound healing To characterize the spatio-temporal evolution of the inflammatory signal c, the fibroblast density r, the deformation w , and the collagen content w, we solve the global balance Equations (13), (16), and (19) along with the local evolution Equation (26).

4.1. Continuous residuals To derive the algorithmic solution, we begin by restating the global balance Equations (13), (16), and (19) in their residual forms and introduce the chemical, biological, and mechanical residuals Rc , Rf , and Rw throughout the domain of interest B, R ¼ c_ þ div q 2 f ¼ _ 0 r r r R ¼ r_ þ div q 2 f ¼ _ 0 c

c

forms of the chemical, biological, and mechanical problems, we multiply the residual statements (27) and the corresponding Neumann boundary conditions with the scalar- and vector-valued test functions dc, dr, and dw and integrate them over the domain B.

4.2.

Discretization in times

To discretize the weak forms of the residual statements (27) in time, we partition the time interval of interest T into nstp nstep 21 subintervals ½tn ; t
as T ¼ 1.8 1.6 – 1.8 1.4 – 1.6 1.2 – 1.4 1.0 – 1.2 0.8 – 1.0 0.6 – 0.8 0.4 – 0.6 0.2 – 0.4 < 0.2

0.6 0.4 0.2 0.0

0

2

4

6

8

10

Normalized time

Figure 5. Wound healing across the temporal scales. Sensitivity of collagen content w with respect to inflammation-induced collagen synthesis rate a. Increasing the collagen synthesis rate induces an overproduction of collagen associated with hypertrophic scarring.

2009). This response is markedly different, however, from the reported overshoot in the cell population days after the peak of inflammation (Mi et al. 2007). Our model based on Equation (4) lies between these two curves and captures both trends. Figure 7(left) illustrates the collagen content. Most existing models follow a very similar profile (Mi et al. 2007; Cumming et al. 2009). The collagen density increases monotonically during the first ten days of healing and then plateaus. Our simulation based on Equation (7) nicely captures these basic features. This set of homogeneous examples provides confidence that our constitutive equations indeed capture a 1.2 Normalized concentration

Downloaded by [University of Otago] at 20:21 01 November 2015

11

1.0 k2 > 1.8 1.6 – 1.8 1.4 – 1.6 1.2 – 1.4 1.0 – 1.2 0.8 – 1.0 0.6 – 0.8 0.4 – 0.6 0.2 – 0.4 < 0.2

0.8 0.6 0.4 0.2 0.0 0

2

4

6

8

10

Normalized time

Figure 6. Wound healing across the temporal scales. Sensitivity of collagen content w with respect to inflammation-induced mitotic rate k2 . Increasing the mitotic rate induces an overproduction of collagen associated with hypertrophic scarring.

5.2. Wound healing across the spatio-temporal scales Now, we explore the spatio-temporal evolution of the chemical, biological, and mechanical variables in a heterogeneous three-dimensional setting. In contrast to the first problem, this now allows us to probe the constitutive equations for the chemical, biological, and mechanical flux terms q c , q r , and s defined in Equations (14), (17), and (22) and perform sensitivity analyses with respect to the associated material parameters. We idealize the tissue sample as a rectangular prism and model the wound as an elliptical enclosure at its center (Wyczalkowski et al. 2013). The tissue has dimensions of 4 £ 4 £ 1 cm3 . Since the problem has two planes of symmetry, we discretize a quarter of the system using trilinear brick elements. The boundary conditions are the same for all examples of this subsection. For the chemical and biological problems, we assume homogeneous Neumann boundary conditions. These boundary conditions imply that the problem remains local, such that chemical signals and cell populations do not cross the boundary within the time period of interest. For the mechanical problem, we impose a constant pre-strain of 10% along the x-direction by solving a mechanical equilibrium problem with the appropriate displacement boundary condition before the onset of injury. This boundary condition resembles the pre-stretched state of skin in vivo. In addition, we apply symmetric boundary conditions to reflect the two planes of symmetry. In addition to the material parameters for the source terms described in detail for the homogeneous problem in Section 5.1, we now need to specify the material parameters for the flux terms. Table 1 summarizes the diffusion coefficients for the chemical and biological fields. For the mechanical problem, the Lame´ constants are l ¼ 0.385 MPa and m ¼ 0.254 MPa, and the Holzapfel parameters are c1 ¼ 0.15 MPa, c2 ¼ 0.0418, and k ¼ 0.05 as calibrated from experiments in pig skin (Jor et al. 2011). The collagen fiber orientation is n ¼ ½1; 0; 0
t . The initial conditions for the chemical, biological, and mechanical fields are heterogeneous, with similar values as in Section 5.1 inside an elliptical wounded region and baseline values outside. We choose the center of the wound at ½xc ; yc ; zc
, and parameterize the injured as ½x 2 xc
2 =r 2x þ ½y 2 yc
2 =r 2y , 1 and z 2 zc , 0:5. The injured region initially has an elevated inflammatory signal, c ¼ 1, and is completely depleted of fibroblasts, r ¼ 0, and collagen, w ¼ 0. The healthy tissue outside the wound is free of inflammation, c ¼ 0, and has a baseline fibroblast density, r ¼ r0 ¼ 0:5, and collagen content, w ¼ 1.

12

A. Buganza Tepole and E. Kuhl ρ : fibroblasts

Normalized concentration

c : inflammation

w : collagen Eq. (4) Ref. (15) Ref. (42)

Eq. (1) Ref. (15) Ref. (42) Ref. (51) Ref. (44)

Eq. (7) Ref. (15) Ref. (42)

Normalized time

Normalized time

Normalized time

Downloaded by [University of Otago] at 20:21 01 November 2015

Figure 7. Comparison of temporal evolution of inflammatory signal, fibroblast density, and collagen content based on local versions of Equations (1), (4), and (7) with existing models and experimental data reported in the literature. Our model captures the overall trends: an exponential decay in the inflammatory signal, a rapid increase in the cell population with a possible overshoot, and a monotonic increase in the collagen concentration that approaches the normalized homeostatic concentration within a few days after inflammation.

Figure 8 shows the spatio-temporal and tempo-spatial evolution of the inflammatory signal c, the fibroblast density r, the Green –Lagrange strain Exx, and the collagen content w for a circular wound with a radius of 1 cm. The overall behavior is similar to that of the homogeneous wound depicted in Figure 3: an elevated inflammatory signal c increases the fibroblast density r. This increases the collagen content w and the tissue stiffness, which gradually reduces the strain Exx, see Figure 8(a). The differences between the individual curves in each graph reflect the regional variation across the wound. These differences disappear over time as the injured region gradually recovers its healthy state, see Figure 8(b). The most distinguishing feature of our model is the inclusion of common mechanical features such as deformation, stress, and strain. The third row of Figure 8(a) displays the spatio-temporal evolution of the Green – Lagrange strain Exx in the direction of the collagen fibers. The initial pre-strain of 10%, applied at the edges of the wound, generates an initially heterogeneous strain profile. The strain distribution results from the heterogeneous tissue stiffness introduced through the regionally varying collagen content, with no collagen in the wounded region, w ¼ 0, and baseline values around the wound, w ¼ 1. As healing progresses, the distribution of the strains becomes more and more homogenous as the collagen content in the wound gradually returns to its baseline value: the tissue gradually recovers its healthy material properties. Here we assume that the tissue is naturally prestrained prior to injury (Buganza Tepole et al. 2014). Including the active stress of fibroblasts would create an additional heterogeneity in the strain profiles as healing progresses, an effect which we have neglected here. Moreover, we assume that the collagen fibers maintain their orientation throughout the healing process. Next, we perform a systematic sensitivity analysis to explore the impact of the wound size and shape on the healing process. First, we vary the size of the injured region while maintaining its circular shape. We study a

larger wound with a radius of 1.5 cm and a smaller wound with a radius of 0.5 cm. Then, we vary the wound shape by changing the degree of ellipticity. We study a moderately elliptic wound with an aspect ratio of 3:2 and an elongated wound with an aspect ratio of 3:1. Figure 9 displays the spatio-temporal and tempospatial evolution of the inflammatory signal c, the fibroblast density r, the Green – Lagrange strain Exx, and the collagen content w for a large and small circular wound with a radii of 1.5 cm and 0.5 cm. The large wound displays larger strain variations Exx than the small wound indicating that it takes longer to heal, see Figure 9(a). The elevated inflammatory signal c in the large wound takes longer to decay than in the small wound, which confirms this trend, see Figure 9(b). For the chosen set of material parameters, the time course of healing is only marginally affected by the wound size. Mathematically, this implies that the evolution equations are dominated by local source rather than global flux terms, and diffusion plays a minor role. Biologically, this implies that, for small wounds on the order of one centimeter, the wound size does not affect the recovery time of the wound as a whole. Figure 10 displays the spatio-temporal and tempospatial evolution of the inflammatory signal c, the fibroblast density r, the Green – Lagrange strain Exx , and the collagen content w for a moderately elliptical and elongated wound with aspect ratios of 3:2 and 3:1. The moderately elliptical wound displays larger strain variations Exx than the elongated wound indicating that it takes longer to heal, see Figure 10(a). The elevated inflammatory signal c in the moderately elliptical wound takes longer to decay than in the elongated wound, which confirms this trend, see Figure 10(b). For the chosen set of material parameters, the time course of the healing process is only marginally affected by the wound shape. Similar to the previous example of varying wound sizes, in which the evolution equations are dominated by local source rather than global flux terms, and diffusion plays a minor role. Due to the lack of experimental data, the fully three-

Computer Methods in Biomechanics and Biomedical Engineering

13 Z

(a)

t = 0 [days]

t = 1 [days]

t = 2 [days]

t = 3 [days]

Y

X

t = 4 [days]

0.9 0.545889 0.331106 0.20083 0.121812 0.0738843 0.044814 0.0271817 0.0164869 0.01

c

0.52 0.477273 0.434545 0.391818 0.349091 0.306364 0.263636 0.220909 0.178182 0.135455 0.0927273 0.05

Downloaded by [University of Otago] at 20:21 01 November 2015

ρ

0.11 0.108 0.106 0.104 0.102 0.1 0.098 0.096 0.094 0.092 0.09

Exx

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

w

(b)

A:(0.00,0,1) B:(0.25,0,1) C:(0.50,0,1) D:(0.75,0,1) E:(1.00,0,1)

ρ

c 1.0

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

0.0

0.0 0

1

2

3

4

5

6

0

t [days]

1

2

3

4

t [days]

5

6

w

Exx

0.118 0.116 0.114 0.112 0.110 0.108 0.106 0

1

2

3

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 4

t [days]

5

6

0

1

2

3

4

5

6

t [days]

Figure 8. Wound healing across the spatio-temporal scales for a circular wound. (a) Spatio-temporal evolution and (b) Tempo-spatial evolution of inflammatory signal c fibroblast density r, Green – Lagrange strain Exx in collagen fiber direction, and collagen content w. An elevated inflammatory signal c increases the fibroblast density r. This increases the collagen content w and the tissue stiffness, which gradually reduces the strain Exx .

dimensional calibration of our model remains challenging. Longitudinal measurements of the wound area would be critical to truly calibrate our model (Roy et al. 2009). However, the evolution of the wound area is a combined result of re-epithelialization and contraction, which we do not consider here. In addition to these area measurements, recent experiments have highlighted the direct correlation between strains beyond the physiological regime and the degree of fibrosis (Gurtner et al. 2011). The proposed framework allows us to capture a spatial variation of strains and a logical next step will be to calibrate the collagen deposition based on these experimental observations.

6.

Discussion

Hypertrophic scarring is a cutaneous condition characterized by the excessive deposition of collagen, which gives

rise to red, thick, stiff, and sometimes painful scar tissue (Chambert et al. 2012). In physiological wound healing, the production of new collagen and breakdown of old collagen balance one another and the overall collagen content remains constant. In pathological wound healing, however, collagen production dominates collagen breakdown and the overall amount of collagen increases. Fortunately, hypertrophic scars do not extend beyond the initial wounded region, but they may continue to grow for weeks or even months (Gurtner et al. 2008). Mechanics has long been neglected to play a crucial role in scar formation (Agha et al. 2011). It is well-known that an increased collagen deposition increases tissue stiffness, and might result in impaired motion when the wound is located close to a joint (Verhaegen et al. 2009). In regenerative medicine, surgeons now manipulate the mechanical microenvironment to their advantage: they minimize scaring through negative pressure wound

14

A. Buganza Tepole and E. Kuhl Z

t = 0 [days]

(a)

t = 1 [days]

t = 2 [days]

t = 3 [days]

t = 4 [days]

Y

X

Exx r = 1.5

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

w

0.11 0.108182 0.106364 0.104545 0.102727 0.100909 0.0990909 0.0972727 0.0954545 0.0936364 0.0918182 0.09

Exx

0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45

Downloaded by [University of Otago] at 20:21 01 November 2015

r = 0.5 w

(b)

r = 1.5

A:(0.00,0,1) B:(0.50,0,1) F:(0,0.50,1) C:(1.00,0,1) G:(0,1.00,1) D:(1.50,0,1) H:(0,1.50,1) E:(2.00,0,1) I:(0,2.00,1)

ρ

c

1

2

3

4

5

6

0.115 0.110 0.105 0.100 0.095 0.090 0.085 0.080 0

1

t [days] A:(0.00,0,1) B:(0.50,0,1) F:(0,0.50,1) C:(1.00,0,1) G:(0,1.00,1) D:(1.50,0,1) H:(0,1.50,1) E:(2.00,0,1) I:(0,2.00,1)

2

2

3

5

6

0

1

4

5

6

0.6 0.5 0.4 0.3 0.2 0.1 0.0

0

1

t [days]

2

3

2

3

4

5

6

0

1

t [days]

ρ

c

1

4

1.2 1.0 0.8 0.6 0.4 0.2 0.0

t [days]

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0

3

2

4

5

6

3

4

5

6

5

6

t [days]

w

Exx

t [days]

0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

w

Exx

0.6 0.5 0.4 0.3 0.2 0.1 0.0

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0

r = 0.5

0.11 0.108182 0.106364 0.104545 0.102727 0.100909 0.0990909 0.0972727 0.0954545 0.0936364 0.0918182 0.09

0.125 0.120 0.115 0.110 0.105 0.100 0.095

1.2 1.0 0.8 0.6 0.4 0.2 0.0 0

1

2

3

4

t [days]

5

6

0

1

2

3

4

t [days]

Figure 9. Wound healing across the spatio-temporal scales for varying wound sizes with radii of r ¼ 1:5 and r ¼ 0:5. (a) Spatiotemporal evolution and (b) Tempo-spatial evolution of inflammatory signal c fibroblast density r, Green – Lagrange strain Exx in collagen fiber direction, and collagen content w. The large wound displays larger strain variations Exx than the small wound indicating that it takes longer to heal. The elevated inflammatory signal c in the large wound takes longer to decay than in the small wound, which confirms this trend.

therapy using vacuum-assisted closure devices (Agha et al. 2011) or through controlled stress shielding using prestrained polymeric patches (Gurtner et al. 2011). A computational model could help identify optimal pressure or pre-strain ranges to accelerate wound healing and reduce scarring. Here, we have established a novel computational framework for the chemo-bio-mechanics of wound healing to understand the fundamental mechanisms of scar formation. Our novel approach toward simulating wound healing is unconditionally stable, geometrically flexible, and conceptually modular. Unconditional stability is guaranteed by the use of an implicit backward Euler scheme to discretize the evolution equations in time, both globally and locally. Using implicit

time integration schemes is algorithmically robust and allows for larger time steps than explicit schemes. For the solution of the resulting nonlinear system of equations, we suggest an incremental iterative Newton – Raphson scheme, again both globally and locally. While the generic equations of wound healing can be bi-directionally coupled, here we have focused on a unidirectionally coupled model problem. For this specific case, we could have used a sequential solution algorithm. However, we are currently in the process of introducing bi-directional coupling. To advance all fields simultaneously in time, it proves convenient to adopt a Newton –Raphson-based solution strategy. The conceptual advantage of Newton– Raphson schemes is that they are not only computationally

Computer Methods in Biomechanics and Biomedical Engineering

15 Z

t=0 [days]

(a)

t=1 [days]

t=2 [days]

t=3 [days]

t=4 [days] Y

X

0.11 0.107778 0.105556 0.103333 0.101111 0.0988889 0.0966667 0.0944444 0.0922222 0.09

Exx

3:2

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

w

0.11 0.107778 0.105556 0.103333 0.101111 0.0988889 0.0966667 0.0944444 0.0922222 0.09

Exx

Downloaded by [University of Otago] at 20:21 01 November 2015

3:1

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

w

(b)

A:(0.00,0,1) F:(0,0.50,1) C:(1.00,0,1) G:(0,1.00,1) D:(1.50,0,1) H:(0,1.50,1) E:(2.00,0,1) I:(0,2.00,1)

3:2

ρ

c

B:(0.50,0,1)

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

0

1

2

3

4

5

6

0.6 0.5 0.4 0.3 0.2 0.1 0.0

0

1

t [days]

2

3

4

5

6

0.115 0.110 0.105 0.100 0.095 0.090 0.085 0.080

w

Exx

0

1

t [days]

2

3

1.2 1.0 0.8 0.6 0.4 0.2 0.0 4

5

6

0

1

t [days]

2

3

4

5

6

5

6

t [days]

A:(0.00,0,1) F:(0,0.50,1) C:(1.00,0,1) G:(0,1.00,1) D:(1.50,0,1)

3:1

H:(0,1.50,1) E:(2.00,0,1) I:(0,2.00,1)

ρ

c

B:(0.50,0,1)

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

0.6 0.5 0.4 0.3 0.2 0.1 0.0 0

1

2

3

4

5

6

0

1

t [days]

2

3

4

t [days]

5

6

0.115 0.110 0.105 0.100 0.095 0.090 0.085 0.080

w

Exx

0

1

2

3

4

t [days]

5

6

1.2 1.0 0.8 0.6 0.4 0.2 0.0

0

1

2

3

4

t [days]

Figure 10. Wound healing across the spatio-temporal scales for varying wound ellipticities with aspect ratios of 3 : 2 and 3 : 1. (a) Spatio-temporal evolution and (b) Tempo-spatial evolution of inflammatory signal c fibroblast density r, Green – Lagrange strain Exx in collagen fiber direction, and collagen content w. The moderately elliptical wound displays larger strain variations Exx than the elongated wound indicating that it takes longer to heal. The elevated inflammatory signal c in the moderately elliptical wound takes longer to decay than in the elongated wound, which confirms this trend.

efficient, but they can easily be supplemented by ad hoc time adaptive schemes, which simply adjust the times step based on the number of required Newton iterations. Geometrical flexibility is a crucial novelty of the proposed model. Existing models have mainly been restricted to zero-, one-, and two-dimensional approximations (Schugart et al. 2008). Our general threedimensional setting allows us to move forward in the spatial complexity. It is a pivotal step toward the simulation of healing in patient-specific geometries (Bo¨l et al. 2011). We achieve this flexibility by using a finite element discretization (Javierre et al. 2009). As opposed to conventional finite volume or finite difference techniques,

finite elements, by design, allow for arbitrary geometries (Zo¨llner et al. 2012; Buganza Tepole et al. 2014). For the first time, we have simulated the healing process in an arbitrary three-dimensional domain. For the sake of illustration, we have used idealized circular and elliptical wound geometries (Wyczalkowski et al. 2013). The extension to more realistic geometries is, of course, straightforward and part of our current work. Conceptual modularity allows us to adjust our approach to other existing models (Cumming et al. 2009) or to expand on the particular model proposed here (Sun et al. 2009; Xue et al. 2009). We have systematically divided the problem of wound healing into three building

Downloaded by [University of Otago] at 20:21 01 November 2015

16

A. Buganza Tepole and E. Kuhl

blocks: chemical, biological, and mechanical (Buganza Tepole and Kuhl 2013). The chemical fields obey a system of partial differential equations common to all reaction – diffusion systems. The biological fields follow a more complex system of partial differential equations that can be specialized for the individual cell populations involved in the healing process. The mechanical fields fall into two categories, global and local, characterized through systems of partial and ordinary differential equations wellestablished for the continuum mechanics of soft biological tissues. We have highlighted the constitutive coupling between these three different fields for general chemo-biomechanical problems. Within this generic setup, we have specified particular constitutive equations to modelspecific aspects of wound healing (Sherratt and Murray 1990). For conceptual simplicity, here we have focused only on studying the impact of biology on mechanics through the collagen deposition by fibroblast. In future, we will also include active mechanical effects such as tissue contraction. We will also include the impact of mechanics on biology. This will allow us to simulate effects of mechanotaxis through an additional flux term for mechanically guided diffusion and of mechanotransduction through an additional source term for mechanically induced mitosis and apoptosis (Wong et al. 2011; Zo¨llner et al. 2013). Along these lines, we have considered tissues to deform quasi statically; yet, we could also include the effect of passive convection through deformation-dependent chemical and biological flux terms (Javierre et al. 2009; Vermolen and Javierre 2012). Here, our elemental model problem depends on only a few parameters and thus allows for systematic parametric experiments. By identifying the critical players in the healing process, it is our goal to iteratively enhance the constitutive equations to investigate the pathways of pathological scarring. We have seen that the collagen deposition rate can have an impact well beyond the transient phase of the healing process. Concomitantly, hypertrophic scaring is characterized by increased collagen concentrations, which also prevail for months after injury. Refining the constitutive equation for the collagen deposition by incorporating growth factors and mechanical cues is a logical next step to model hypertrophic scaring. Our current approach not only explores relevant healing scenarios for a particular model, but effectively creates a generic framework that can be easily expanded to incorporate other features such as the impact of mechanical cues on cell mitosis or apoptosis. As such, it is not only applicable to explore chemo-bio-mechanical interaction during wound healing in skin, but also in other inflammation-based systems, for example in healing infarcts in cardiac muscle (Rouillard and Holmes 2012). Ultimately, a better understanding of the healing mechanisms in living systems can inspire the design of

novel, self-healing engineering systems (Mergheim and Steinmann 2013). In addition to these algorithmic aspects, our model accounts for a state-of-the-art mechanical characterization of skin (Limbert and Simms 2013) within a continuum mechanics approach (Jor et al. 2011; Buganza Tepole et al. 2014). Recently there has been significant development in the theory of the mechanics of living soft collagenous tissues (Holzapfel et al. 2000; Buganza Tepole et al. 2011; Saez et al. 2014). Unfortunately, these advances have been almost entirely disconnected from recent trends in systems biology, which have been confined to either rigid geometriesor viscoelastic fluids (Xue et al. 2009). These simplifications impose great limitations toward understanding the role of mechanical cues during wound healing. A rigorous, accurate mechanical characterization is a fundamental knowledge gap in existing models for wound healing. Here, we characterize skin using a hyperelastic strain energy function parameterized in terms of a set of microstructure variables such as collagen orientation (Kuhl and Holzapfel 2007) and collagen content (Saez et al. 2013). The importance of time-varying material properties has recently been identified as a critical aspect in wound healing (Ben Amar and Wu 2014). By allowing our microstructural variables to evolve in time, we establish clear relations between the action of the different cell populations and tissue remodeling (Menzel and Kuhl 2012). Recent studies also provide evidence that pre-strain and tissue tension have a significant effect on the healing characteristics of circular and elliptical wounds (Wyczalkowski et al. 2013; Buganza Tepole et al. 2014). Our model allows us to impose physiological boundary conditions such as pre-strain (Rausch and Kuhl 2013) and predict the spatio-temporal evolution of tissue tension across arbitrarily shaped wounds throughout the healing process. In conclusion, the proposed framework introduces a new generation of wound healing models that may provide fundamental insight into the role of mechanics in scar formation. A unified monolithic finite element treatment of the underlying chemical, biological, and mechanical fields is a first step toward the smooth incorporation of realistic environmental conditions and personalized individual geometries. Our model has the potential to significantly improve effective wound management and optimize treatment options on a patient-specific basis. Funding This work was supported by the Consejo Nacional de Ciencia y Tecnologia CONACyT Fellowship and the Stanford Graduate Fellowship to Adria´n Buganza Tepole and by the National Science Foundation CAREER award CMMI 0952021 and INSPIRE award 1233054 and the National Institutes of Health grant [grant number U01 HL119578] to Ellen Kuhl.

Computer Methods in Biomechanics and Biomedical Engineering Conflict of interest disclosure statement No potential conflict of interest was reported by the author(s).

Downloaded by [University of Otago] at 20:21 01 November 2015

References Aarabi S, Bhatt KA, Shi Y, Paterno J, Chang EI, Loh SA, Holmes JW, Longaker MT, Yee H, Gurtner GC. 2007. Mechanical load initiates hypertrophic scar formation through decreased cellular apoptosis. FASEB J. 21:3250 – 3261. Agha R, Ogawa R, Pietramaggiori G, Orgill DP. 2011. A review of the role of mechanical forces in cutaneous wound healing. J Surg Res. 171:700 – 708. Bayat A, McGrouther DA, Ferguson MW. 2003. Skin scarring. BMJ. 326:88– 92. Ben Amar M, Wu M. 2014. Re-epithelialization: advancing epithelium frontier during wound healing. J R Soc Interface., doi:10.1098/rsif.2013.1038. Bo¨l M, Sturmat M, Weichert C, Kober C. 2011. A new approach for the validation of skeletal muscle modeling using MRI data. Comp Mech. 47:591– 601. Bowes LE, Jimenez MC, Hiester ED, Sacks MS, Brahmatewari J, Mertz PE, William H. 1999. Collagen fiber orientation as quantified by small angle light scattering in wounds treated with transforming growth factor-beta2 and its neutalizing antibody. Wound Repair Regen. 7:179– 186. Brown BC, McKenna SP, Siddhi K, McGrouther DA, Bayat A. 2008. The hidden cost of skin scars: quality of life after skin scarring. J Plast Reconstr Aesthet Surg. 61:1049 – 1058. Buganza Tepole A, Gart M, Gosain AK, Kuhl E. 2014. Characterization of living skin using multi view stereo and isogeometric analysis. Acta Biomater. 10:4822 – 4831. Buganza Tepole A, Gosain AK, Kuhl E. 2012. Stretching skin: the physiological limit and beyond. Int J Nonlin Mech. 47:938 – 949. Buganza Tepole A, Gosain AK, Kuhl E. 2014. Computational modeling of skin: using stress profiles as predictor for tissue necrosis in reconstructive surgery. Comput Struct. 143:32 – 39. Buganza Tepole A, Kuhl E. 2013. Review: systems-based approaches towards wound healing. Pediatr Res. 73:553 –563. Buganza Tepole A, Ploch PC, Wong J, Gosain AK, Kuhl E. 2011. Growing skin: a computational model for skin expansion in reconstructive surgery. J Mech Phys Solid. 59:2177 – 2190. Buganza Tepole A, Steinberg JP, Kuhl E, Gosain AK. 2014. Application of finite element modeling to optimize flap design with tissue expansion. Plast Reconstr Surg. 134:785 – 792. Chambert J, Zhao L, Remache D, Jacquet E. 2012. Numerical analysis of keloid scar in the presternal area. Comput Method Biomech Biomed Eng. 15:23 –24. Chary SR, Jain RK. 1989. Direct measurement of interstitial convection and diffusion of albumin in normal and neoplastic tissues by fluorescence photobleaching. Proc Natl Acad Sci USA. 86:5385– 5389. Chiquet M, Gelman L, Lutz R, Maier S. 2009. From mechanotransduction to extracellular matrix gene expression in fibroblasts. Biochim Biophys Acta. 1793:911– 920. Cumming BD, McElwain DLS, Upton Z. 2009. A mathematical model of wound healing and subsequent scarring. J R Soc Interface. 7:19 – 34. Dallon JC, Sherratt JA, Maini PK. 2001. Modeling the effects of transforming growth factor-beta on extracellular matrix alignment in dermal wound repair. Wound Repair Regen. 9:278 – 286.

17

Diegelmann RF, Evans MC. 2004. Wound healing: an overview of acute, fibrotic and delayed healing. Front Biosci. 9:283 –289. Flynn C, Taberner A, Nielsen P. 2011. Mechanical characterisation of in vivo human skin using a 3D force-sensitive micro-robot and finite element analysis. Biomech Model Mechanobiol. 10:27– 38. Garikipati K, Olberding JE, Narayanan H, Arruda EM, Grosh K, Calve S. 2006. Biological remodeling: stationary energy, configurational change, internal variables and dissipation. J Mech Phys Solid. 54:1493– 1515. Garzon-Alvarado DA, Sandoval RPC, Acosta RCV. 2012. A mathematical model of medial collateral ligament repair: migration, fibroblast proliferation and collagen formation. Comput Meth Biomech Biomed Eng. 15:571– 583. Gurtner GC, Dauskardt RH, Wong VW, Bhatt KA, Wu K, Vial IN, Padois K, Korman JM, Longaker MT. 2011. Improving cutaneous scar formation by controlling the mechanical environment. Ann Surg. 254:217– 225. Gurtner GC, Werner S, Barrandon Y, Longaker MT. 2008. Wound repair and regeneration. Nature. 453:314– 321. Herbert SP, Stainier DY. 2011. Molecular control of endothelial cell behaviour during blood vessel morphogenesis. Nat Rev Mol Cell Biol. 12:551 – 564. Himpel G, Menzel A, Kuhl E, Steinmann P. 2008. Timedependent fiber reorientation of transversely isotropic continua: finite element formulation and consistent linearization. Int J Number Method Eng. 73:1413 –1433. Holzapfel GA, Gasser TC, Ogden RW. 2000. A new constitutive framework for anterial wall mechanics and a comparative study of material models. J Elast. 61:1 –48. Hunter PJ, Borg TK. 2003. Integration from proteins to organs: the physiome project. Nat Rev Mol Cell Biol. 4:237– 243. Javierre E, Moreno P, Doblare M, Garcia-Aznar JM. 2009. Numerical modeling of a mechano-chemical theory for wound contraction analysis. Int J Solid Struct. 46:3597– 3606. Javierre E, Vermolen FJ, Vuik C, Zwaag S. 2009. A mathematical analysis of physiological and morphological aspects of wound closure. J Math Biol. 59:605 – 630. Jor JWY, Nash MP, Nielsen PMF, Hunter PJ. 2011. Estimating material parameters of a structurally based constitutive relation for skin mechanics. Biomech Model Mechanobiol. 10:767 – 778. Kitano H. 2002. Computational systems biology. Nature. 420:206 – 210. Kuhl E, Holzapfel GA. 2007. A continuum model for remodeling in living structures. J Mater Sci. 42:8811 – 8823. Kuhl E, Steinmann P. 2004. Computational modeling of healing: an application of the material force method. Biomech Model Mechanobiol. 2:187– 203. Lanir Y. 1983. Constitutive equations for fibrous connective tissues. J Biomech. 16:1 – 12. Laurent GJ. 1987. Dynamic state of collagen: pathways of collagen degradation in vivo and their possible role in regulation of collagen mass. Am J Physiol Cell Physiol. 252:1 –9. Limbert G, Simms C. 2013. Special issue on skin mechanobiology. J Mech Behav Biomed Mater. 28:395 – 396. Maini PK, McElwain DL, Leavesley DI. 2004. Traveling wave model to interpret a wound-healing cell migration assay for human peritoneal mesothelial cells. Tissue Eng. 10:475– 482. Martin P. 1997. Wound healing-aiming for perfect skin regeneration. Science. 276:75 – 81. McDougall S, Dallon J, Sherratt J, Maini P. 2006. Fibroblast migration and collagen deposition during dermal wound

Downloaded by [University of Otago] at 20:21 01 November 2015

18

A. Buganza Tepole and E. Kuhl

healing: mathematical modelling and clinical implications. Philos Trans A Math Phys Eng Sci. 364:1385 – 1405. Menzel A. 2007. A fibre reorientation model for orthotropic multiplicative growth. Biomech Model Mechanobiol. 6:303 – 320. Menzel A, Kuhl E. 2012. Frontiers in growth and remodeling. Mech Res Commun. 42:1– 14. Mergheim J, Steinmann P. 2013. Phenomenological modelling of self-healing polymers based on integrated healing agents. Comput Mech. 52:681 –892. Mi Q, Riviere B, Clermont G, Steed DL, Vodovotz Y. 2007. Agentbased model of inflammation and wound healing: insights into diabetic foot ulcer pathology and the role of transforming growth factor-beta1. Wound Repair Regen. 15:671–682. Mutsaers SE, Bishop JE, McGrouther G, Laurent GJ. 1997. Mechanisms of tissue repair: from wound healing to fibrosis. Int J Biochem Cell Biol. 29:5– 17. Olsen L, Sherratt JA, Maini P. 1995. A mechanochemical model for adult dermal wound contraction: on the permanence of the contracted tissue displacement profile. J Theor Biol. 177:113 – 128. Olutoye OO, Zhu X, Cass DL, Smith CW. 2005. Neutrophil recruitment by fetal porcine endothelial cells: implications in scarless fetal wound healing. Pediatr Res. 58: 1290– 1294. Paterno J, Vial IN, Wong VW, Rustad KC, Sorkin M, Shi Y, Bhatt KA, Thangarajah H, Glotzbach JP, Gurtner GC. 2011. Akt-mediated mechanotransduction in murine fibroblasts during hypertrophic scar formation. Wound Repair Regen. 19:49 – 58. Pettet GJ, Byrne HM, McElwain DL, Norbury J. 1996. A model of wound-healing angiogenesis in soft tissue. Math Biosci. 136:35 – 63. Qutub AA, Mac Gabhann F, Karagiannis ED, Vempati P, Popel AS. 2009. Multiscale models of angiogenesis. IEEE Eng Med Biol Mag. 28:14 – 31. Rausch MK, Kuhl E. 2013. On the effect of prestrain and residual stress in thin biological membranes. J Mech Phys Solids. 61:1955 – 1969. Rouillard AD, Holmes JW. 2012. Mechanical regulation of fibroblast migration and collagen remodeling in healing myocardial infarcts. J Phyisol. 590:4585 – 4602. Roy S, Biswas S, Khanna S, Gordillo G, Bergdall V, Green J, Marsh CB, Gould LJ, Sen CK. 2009. Characterization of a preclinical model of chronic ischemic wound. Physiol Genomics. 37:211– 224. Saez P, Pena E, Martinez MA, Kuhl E. 2013. Mathematical modeling of collagen turnover in biological tissue. J Math Biol. 67:1765– 1793. Saez P, Pena E, Martinez MA, Kuhl E. 2014. Computational modeling of hypertensive growth in the human carotid artery. Comput Mech. 53:1183 – 1196. Schugart RC, Friedman A, Zhao R, Sen CK. 2008. Wound angiogenesis as a function of tissue oxygen tension: a mathematical model. Proc Nat Acad Sci. 105:2628– 2633. Sherratt JA, Murray JD. 1990. Models of epidermal wound healing. Proc R Soc B Bio Sci. 241:29– 36. Simpson CL, Patel DM, Green KJ. 2011. Deconstructing the skin: cytoarchitectural determinants of epidermal morphogenesis. Nat Rev Mol Cell Biol. 12:565– 580.

Southern J, Pitt-Francis J, Whiteley J, Stokeley D, Kobashi H, Nobes R, Kadooka Y, Gavaghan D. 2008. Multi-scale computational modelling in biology and physiology. Prog Biophys Mol Biol. 96:60 –89. Sun T, Adra S, Smallwood R, Holcombe M, MacNeil S. 2009. Exploring hypotheses of the actions of TGF-beta1 in epidermal wound healing using a 3D computational multiscale model of the human epidermis. PLoS ONE. 4: e8515. Tranquillo RT. 1987. Stochastic model of leukocyte chemosensory movement. J Math Biol. 25:229 – 262. Tranquillo R, Murray J. 1992. Continuum model of fibroblastdriven wound contraction – inflammation-mediation. J Theor Biol. 158(1992):135 – 172. Tranquillo RT, Murray JD. 2007. Continuum model of fibroblastdriven wound contraction: inflammation-mediation. Biomech Model Mechanobiol. 158:361 – 371. Velnar T, Bailey T, Smrkolj V. 2009. The wound healing process: an overview of the cellular and molecular mechanisms. J Int Med Res. 37:1528– 1542. Verhaegen PDHM, van Zuijlen PPM, Pennings NM, van Marle J, Niessen FB, van der Horst CMAM, Middelkoop E. 2009. Differences in collagen architecture between keloid, hypertrophic scar, normotrophic scar, and normal skin: an objective histopathological analysis. Wound Repair Regen. 17:649–656. Vermolen FJ, Javierre E. 2012. A finite-element model for healing of cutaneous wounds combining contraction, angiogenesis and closure. J Math Biol. 65:967 – 996. Vodovotz Y. 2010. Translational systems biology of inflammation and healing. Wound Repair Regen. 18:3– 7. Vodovotz Y, Csete M, Bartels J, Chang S, An G. 2008. Translational systems biology of inflammation. PLoS Comput Biol. 4:e1000014. Werner S, Grose R. 2003. Regulation of wound healing by growth factors and cytokines. Physiol Rev. 83:835– 870. Wong VW, Akaishi S, Longaker MT, Gurtner GC. 2011. Pushing back: wound mechanotransduction in repair and regeneration. J Invest Dermatol. 131:2186 – 2196. Wong VW, Levi K, Akaishi S, Schultz G, Dauskardt RH. 2012. Scar zones: region specific differences in skin tension may determine incisional scar formation. Plast Reconstr Surg. 129:1272 – 1276. Wyczalkowski MA, Varner VD, Taber LA. 2013. Computational and experimental study of the mechanics of embryonic wound healing. J Mech Behav Biomed Mater. 28:125 –146. Xue C, Friedman A, Sen CK. 2009. A mathematical model of ischemic cutaneous wounds. Proc Natl Acad Sci. 106:16782 – 16787. Zo¨llner AM, Buganza Tepole A, Gosain AK, Kuhl E. 2012. Growing skin: tissue expansion in pediatric forehead reconstruction. Biomech Model Mechanobiol. 11:855– 867. Zo¨llner AM, Buganza Tepole A, Kuhl E. 2012. On the biomechanics and mechanobiology of growing skin. J Theor Biol. 297:166 – 175. Zo¨llner AM, Holland MA, Honda KS, Gosain AK, Kuhl E. 2013. Growth on demand: reviewing the mechanobiology of stretched skin. J Mech Behav Biomed Mater. 28:495 – 509.