|Home | About | Journals | Submit | Contact Us | Français|
Wound healing in the pediatric patient is of utmost clinical and social importance, since hypertrophic scarring can have aesthetic and psychological sequelae, from early childhood to late adolescence. Wound healing is a well-orchestrated reparative response affecting the damaged tissue at the cellular, tissue, organ, and system scales. While tremendous progress has been made towards understanding wound healing at the individual temporal and spatial scales, its effects across the scales remain severely understudied and poorly understood. Here we discuss the critical need for systems-based computational modeling of wound healing across the scales, from short-term to long-term and from small to large. We illustrate the state of the art in systems modeling by means of three key signaling mechanisms: oxygen tension regulating angiogenesis and revascularization; TGF-β kinetics controlling collagen deposition; and mechanical stretch stimulating cellular mitosis and extracellular matrix remodeling. The complex network of biochemical and biomechanical signaling mechanisms and the multi-scale character of the healing process make systems modeling an integral tool in exploring personalized strategies for wound repair. A better mechanistic understanding of wound healing in the pediatric patient could open new avenues in treating children with skin disorders such as birth defects, skin cancer, wounds, and burn injuries.
Dermal wound healing in the pediatric patient is a symphony of events, precisely synchronized to repair the damaged tissue, restore its protective barrier function, and safely return it to its homeostatic equilibrium state1. Although the underlying processes, cell-matrix interaction, cell-cell cross-talk, and cellular mechanotransduction, involve a complex cascade of events, dermal wound healing is robust and rarely diverges to malignant transformation2. Yet, it is not always perfect. While pre-natal skin usually heals smoothly to seamlessly restore the state prior to injury3, post-natal skin is incapable of healing wounds tracelessly, leaving scar behind. Post-natal skin can easily restore its protective barrier function; however, the resulting scar rarely has the same microstructure, collagen content, and mechanical properties as the native tissue4. In extreme cases, pronounced fibrotic activity might even initiate hypertrophic scaring, characterized through an excessive collagen deposition5. In the pediatric patient, excessive scarring has consequences throughout early childhood and adolescence, and can lead to low self-esteem or even stigmatization6. The prevalence of hypertrophic scarring in the pediatric population is overwhelming: Of the total cases of burn scars and keloids, 70% occur in children7.
The underlying mechanisms of scar formation are now better understood than ever before, and we have made tremendous progress towards improving and accelerating healing mechanisms8. We have come to appreciate that the healing process spans various temporal and spatial scales and that is affected by both chemical and mechanical cues. However, even with the detailed insight that traditional approaches have provided on the individual scales, the behavior of the system as a whole remains elusive. Computational modeling is increasingly recognized as a powerful tool to provide insight into the dynamics of wound healing and the interaction of biochemical and biomechanical phenomena across the different scales9 .
Fortunately, dermal wound healing, like all inflammation-based processes in the human body, is based on various redundant signals and cross talk between different signaling networks10. While redundancy is hugely beneficial for the biological system itself, it complicates the overall understanding of the healing process: Even if individual elements of the signaling network are well understood in isolation, the coupling of these elements is hugely complex, and it is virtually impossible to gain basic insights based on sparse experimental data. Computational systems biology is of acknowledged importance to advance our holistic understanding of pediatric wound healing: Short term, computational systems biology allows us to systematically explore controlled what-if-scenarios and virtually probe various hypothesis to better understand the healing process as a whole11. Long term, given the incredible variability of healing responses between different individuals, computational systems biology is an integral ingredient to shape the future of personalized medicine12.
From the moment of injury until the tissue reaches its final configuration in the form of a mature scar, months or even years may elapse. However, the protective function of skin has to resume immediately to avoid dehydration, infection, and loss of tissue integrity. Accordingly, the initial phase of healing only takes a few minutes. The process of wound healing is commonly divided into four overlapping phases: hemostasis, inflammation, proliferation, and remodeling, see Figure 1, rows 1 through 4.
Hemostasis occurs in the order of minutes. During hemostasis, the wounded space rapidly fills with a clot to stop blood loss and to reestablish a barrier from the outside world. In the later phases of healing, this clot will serve as a temporary matrix for the cells that migrate towards the wound to reconstruct the dermal tissue (Figure 1.1c)13. The clot is composed mainly of fibrin fibers and platelets. The cytoplasm of the platelets carries α-granules from which a cocktail of growth factors and cytokines is released when the cells degranulate (Figure 1.1b)14. Pro-inflammatory signals from the platelets are directly reflected in molecular changes in the endothelial cells of the blood vessels near the injury site. A chemotactic response attracts leukocytes, such as neutrophils and monocytes (Figure 1.1a), which sense these molecular changes and adhere to the endothelial cells. The ongoing chemotactic pathway and the interaction between leukocytes and endothelial cells lead to the capture and transmigration of the neutrophils and monocytes in a process called diapedesis15.
Neutrophil diapedesis marks the beginning of the inflammatory phase. Neutrophil transmigration occurs also very rapidly; in fact, the first neutrophils reach the wound within a few minutes, but they keep being recruited for the next two or three days. The primary function of neutrophils is to dispose pathogens entrapped in the clot at the moment of injury. Their secondary function is to amplify the pro-inflammatory indicators. Monocytes arrive at the wound within two days of the initial insult and differentiate into macrophages (Figure 1.2a). The role of monocytes is twofold, they phagocyte the remaining matrix debris, pathogens and neutrophils, and they produce soluble chemical mediators such as TGF-β and VEGF. Both TGF-β and VEGF are crucial to guide the recruitment of fibroblasts and endothelial cells later in the healing process (Figure 1.2c).
During the proliferative phase, two processes occur simultaneously: re-epithelialization and replacement of the fibrin clot by granulation tissue (Figure 1.3c). These steps in the repair process continue throughout two to three weeks. In this phase, the sharp chemoattractant gradients generated through the inflammatory phase attract different cell types to manufacture new tissue (Figure 1.3a). At the epidermal level, keratinocytes in the wound edge dissolve their adhesions to the basal lamina in order to crawl and reproduce on top of the fibrin matrix, rebuilding the epidermis over the injured region (Figure 1.3b). On the dermal level, fibroblasts are in charge of depositing collagen, the main load-bearing constituent of skin (Figure 1.3b). The metabolic needs resulting from increased cellular proliferation and migration require the formation of a vast network of capillary tubes to provide nutrients and oxygen. Endothelial cells rapidly form a new vasculature in a process referred to as angiogenesis (Figure 1.3b). Towards the end of this phase, some fibroblasts transform into myofibroblasts, which actively pull on the wound edges to contract the injured tissue.
Once the fibrin matrix has been replaced by granulation tissue made of fibroblasts, the remaining macrophages, and a vascular network, the remodeling phase begins. At this point, the wound is fully closed, yet the tissue has very poor quality. It consists primarily of thick, aligned collagen bundles, instead of the interwoven collagen networks found in native skin (Figure 1.4b). Dynamic changes continue, but on a much slower time scale. Remodeling can go on for months or even years. Finally, the vascular network retracts and most of the cells undergo apoptosis or migrate out of the affected region. The remaining fibroblasts (Figure 1.4a) gradually restore skin integrity and mechanical homeostasis (Figure 1.4c). However, unlike pre-natal skin, post-natal skin usually never fully recovers its original native state.
Similar to virtually all biological phenomena, in wound healing, temporally and spatially interacting events at the molecular, subcellular, cellular, tissue and organ levels ultimately converge to a well-defined system level response16, see Figure 1, columns a through d. This macroscopic response is inherently rooted in the hierarchical structure, which is established through a precisely defined biological organization17. Although it is in principle possible to zoom in to the subcellular or even molecular scales, the cellular level is the preferred starting point in the context of wound healing modeling, mainly because cells are the smallest autonomous building blocks in the biological hierarchy18. One level above in the organization, at the tissue level, we find interacting cell aggregates embedded in the extra-cellular matrix (ECM). The next larger scale corresponds to the organ level, the level of the wound itself. In principle, we could even contemplate higher levels of organization, for example entire organ systems. Here, however, we do not consider larger spatial scales, because all relevant phenomena take place within the local wound environment.
The smallest scale we explore here is the cellular level, a scale in the order of 10um. At this scale, four cell-types are critical to wound healing: endothelial cells, fibroblasts, macrophages, and keratinocytes (Figures 1.1a to 1.4a). For each cell type, there are specific regulatory mechanisms, which are key to a successful repair process. The most important control mechanism is chemotaxis, the migration of cells in the direction of increasing chemical gradients19. Endothelial cells, fibroblasts, and leukocytes migrate into the wound site attracted by the strong concentration of different growth factors and other chemical mediators released during the inflammatory phase20-22. Chemotaxis in eukaryotic cells is complex and, unlike for prokaryotic systems, there are no established models to predict how eukaryotic cells migrate as a function of chemical gradients. Prokaryotic cell migration resembles a biased random walk, which can be easily represented by simple mathematical models. Eukaryotic cells, however, move in continuous paths with smooth turns towards ascending signal concentrations, but display an overall stochastic behavior23-25.
In addition to chemical cues, mechanical cues play a significant role in wound healing26,27. The ultimate goal of skin repair is to reestablish the mechanical load bearing capacity to restore the homeostatic equilibrium state. Key to this restoration is the translation of mechanical signals into chemical reactions, a process known as mechanotransduction, which is the landmark characteristic of fibroblasts in the repair sequence (Figure 1.4a)28. Fibroblasts anchor to the ECM through focal adhesions. These adhesion sites are linked to stretch-activated ion channels across the cell membrane and to the cytoskeleton inside the cell. Altering the mechanical scenario of the ECM directly creates a flow of charged ions via these ion channels and indirectly governs cellular behavior through perturbed cytoskeletal dynamics29. Fibroblasts respond to mechanical changes in the ECM in several ways: by depositing collagen, by decreasing apoptosis rates, by increasing inflammatory signals, and by transforming into myofibroblasts to actively pull on the wound edges30-32.
For endothelial cells and keratinocytes (Figure 1.3a), cell-cell and cell-matrix adhesions are of primary importance throughout the entire repair process. Keratinocytes critically depend on cell-cell and cell-matrix crosstalk to define their position in the epidermal lattice. The polarization of keratinocytes directs lateral migration and proliferation during re-epithelialization33. Endothelial cells also depend on these types of interaction during angiogenesis34. Some cells in the tip of sprouting vessels are activated and protrude filopodial extensions, which actively interact with their immediate microenvironment to guide angiogenesis. The remaining endothelial cells support the leading tip by maintaining the connectivity with the parent vessels, and by helping in the maturation of the newly formed capillary tubes35.
On the next level of the hierarchy, we are interested in the bulk mechanical properties of the ECM, the diffusion of growth factors, and the overall behavior of cellular aggregates (Figures 1.1b to 1.4.b). The properties of the ECM and the diffusion of growth factors are crucial at this scale, because they contribute a key aspect of the healing process: They pass information across distances larger than the characteristic cell size, coordinating the action of multiple cells without using direct cell-cell communication. Soluble mediators diffuse chemical cues and the ECM transmits mechanical cues across the tissue (Figures 1.1b and 1.2b). To characterize chemical diffusion in an isotropic medium, only a single constant needs to be specified. To characterize mechanical signaling of the ECM, however, we can select from a huge variety of constitutive models, and two or even more constants need to be determined (Figure 1.4b)36-38.
Two other phenomena that become important at the tissue level are pattern formation and collective cell migration. An illustrative example is re-epithelialization, associated with the formation of an advancing front that has distinct collective features entirely different from the individual cellular response. The leading edge of keratinocytes that invade the wound adopts a wave propagation profile, which can be characterized using Fisher's equation (Figure 1.3b)39. Another even more intriguing proof of complex organization is angiogenesis, a highly coordinated process that cannot be understood solely from the study of a single cell. During angiogenesis, a new vasculature forms from existing blood vessels that surround the injured region40. This process consists of combined tip chemotactic diffusion and sprouting to create a new network of capillary tubes (Figure 1.3b)41.
The final level of interest in wound healing is the spatial region that contains the injured skin and its immediate surrounding healthy tissue (Figures 1.1c to 1.4c). At this level, we witness the healing capacities of skin as the result of a perfectly organized interaction of all phenomena described above: from microstructural decisions of the individual cells, via mesostructural properties of the cells embedded in the ECM, to macrostructural reaction-diffusion systems of the soluble mediators and the patterning capabilities of cellular aggregates.
Mathematical representations can foster our understanding of biological systems. An abstraction of a biological system consists of characteristic variables, such as the concentration of a chemical substance, or the cell density in a region of tissue, and the interaction rules between the different components of the system. Mathematical models are representations of the real world and, as such, they can be created in many different ways depending on the features of the problem itself, on the information at hand, and on the mathematical tools preferred by the modeler. Creating a mathematical model is only the first step of a simulation process. The second step consists of the converting the mathematical model into a computational model that allows solving complex patient-specific problems.
In the context of wound healing, there are two conceptually different mathematical approaches to represent a biological system, continuum models and discrete models. In the physical world, substances can vary smoothly over a given domain. In this case we represent them through a function called a continuous field. Continuous fields are the basic representation of continuum models. Typical continuous fields in wound healing are concentrations of chemical mediators, such as TGF-β. Some components of the system may be large enough to be considered as discrete entities. To model individual entities, we typically use discrete models that represent the individual components explicitly and characterize their interaction through simple mathematical equations. Typical discrete entities in wound healing are cell populations, such as densities of fibroblasts and macrophages. The choice between continuum and discrete models strongly depends the particular scale of interest. Components that can be considered as discrete entities at the molecular and cellular level become smeared out at the tissue and organ level, where they can be approximated as continuous fields. In general, the scale of interest will determine which modeling approach is more appropriate.
In systems biology, different types of interacting components have to be represented through the mathematical model, either as continuous fields or as discrete entities. A multi-field model is one that consists of several interacting fields. Typical multi-field phenomena in wound healing models are various interacting chemical signaling pathways, such as TGF-β and VEGF. In some cases, it might be advantageous to represent some variables as continuous fields while others might be represented as discrete entities. Accordingly there are three basic types of representations: purely continuum, purely discrete and hybrid continuum/discrete. In addition, a mathematical model for each of these representations can span various spatial scales, from the cellular level, via the tissue and organ levels, to the system level. We characterize these types of models as multi-scale models.
From a modeling point of view, biological systems consist of three basic ingredients: chemical substances, cell populations, and mechanical factors. Chemical substances are generally treated as continuous fields. Cell populations can either be modeled individually or as fields, depending on the scale of interest. Here, we use the notion biological field for fields that represent cell populations. Mechanical factors include forces and deformations and are almost always treated as continuous fields.
In continuum models, the evolution of different fields in space and time is governed by partial differential equations. A partial differential equation establishes: i) how a field interacts with other fields through coupling terms; ii) how a field evolves in time; and iii) how a field varies in space over the domain of interest. One partial differential equation is needed for each field. Continuum models represented through partial differential equations are an elegant and compact way to characterize spatial variations of chemical concentrations, mechanical strains, or mechanical stresses.
In discrete models, the behavior of each entity is governed by one or more ordinary differential equations. An ordinary differential equation establishes: i) how a characteristic variable of a single entity interacts with other variables through coupling terms; and ii) how a characteristic variable evolves in time.
In pediatric wound healing, the mathematical model is typically too complex to be solved analytically. Through a process called discretization, we can transfer the mathematical model, the set of equations, into a computational model, an algorithmic solution of these equations, that can easily simulate phenomena of various degrees of complexity.
Continuum models are typically represented through coupled partial differential equations. To solve the underlying set of equations, it is necessary to subdivide or discretize both space and time. Time is almost always discretized into small intervals by means of the finite difference methods. Discretizing space can be preformed in several ways: The simplest approach to represent the domain of interest is a regular mesh, as done in the finite difference method and the finite volume method. To represent arbitrarily shaped domains such as particular parts of the body, the preferred approach is the use the finite element method. This method approximates the domain by subdividing it into little pieces called elements. While finite elements are geometrically more versatile then finite differences and finite volumes, their algorithms are more sophisticated and require special care of implementation.
Once time and space are discretized, we obtain a system of algebraic equations. The solution of this system characterizes the evolution of all continuous fields, both in space and time. Depending of the required level of accuracy, the solution procedure can become a challenging and computationally expensive task. Consequently, we can select between various useful simplifications: In the real world, the domain of interest is three-dimensional. However, we can often extract useful information about the biological system by approximating it as one or two-dimensional. Another simplification could be a particular symmetry of the system. Typical examples are cylindrical symmetries, which can be approximated by one-dimensional formulation referred to as axisymmetric.
Discrete models are typically represented through ordinary differential equations. These equations do not need to be solved over the entire domain; they are solved locally to predict the evolution of each individual entity in time. Since only the time needs to be partitioned into small intervals, and not the space, it is common to use the finite difference method. While discrete models are conceptually simpler than continuum models, they are computationally more expensive, since every single entity has to be represented through its own equation. When dealing with entire cell populations, their computational cost becomes overwhelming.
From the above discussion on the basic approaches towards systems-based mathematical representations, and to establish a framework for a guided discussion of the most recent wound healing models, we propose a few useful guidelines for effective computational modeling of wound healing.
One-dimensional and axisymmetric models provide a quick first insight into the interplay of individual signaling pathways during wound healing. They are of great value to quickly test therapeutic hypotheses. Because of their conceptual simplicity, one-dimensional and axisymmetric models should be viewed as a validated and well-calibrated starting point for further model refinement. For detailed analyses that focus on the spatio-temporal interaction of signaling pathways in wound healing, however, we recommend using two- or three-dimensional models.
We recommend using discrete models to study the biochemistry of individual cells. Discrete models are typically based on simple relations and straightforward to implement. Since they are computationally expensive, however, they are only feasible to answer fundamental scientific questions at small scales. We recommend fully avoiding discrete models when diffusion or mechanical cues play an important role.
We recommend using continuum models to study the interaction of biochemical and biomechanical fields. Continuum models are computationally efficient, but require special care when translating the collective action of individual cells to a continuous field. When embedded in efficient finite element solvers, continuum models allow for a patient-specific analysis based on individualized, realistic three-dimensional geometries. Continuum models are ideal to answer translational clinical questions at larger scales.
Although models focusing on a single individual aspect of wound healing can provide specific insight, holistic models are currently becoming the gold standard in wound healing simulations. In contrast to single-field models, holistic multi-field models can reliably represent the cross-signaling networks during wound healing. When creating a multi-field model, we recommend to include a representation of at least five interacting fields: i) fibroblast population, ii) endothelial cell population and/or nutrient concentration, iii) indicator of ECM restoration such as collagen content, iv) chemoattractant concentration, and v) inflammatory cells population. This baseline five-field model can, of course, be further refined if specific aspects of the healing process are of interest.
Over the past two decades, mathematical and computational modeling have advanced as key players in the quest for understanding the complex multi-scale phenomena involved in wound healing. Ideally, once calibrated and validated, these models can serve as tools to quantify the impact of selected perturbations to the baseline system, and, ultimately, to predict the effect of different therapeutic treatment options.
In an effort to gain basic scientific insight into the process of wound healing, researchers in the early nineties started to develop mathematical models of the repair response of skin. Sherrat and Murray are often considered as the fathers of modern mathematical modeling of wound healing42. Since their first model was published more than two decades ago, a multitude of modifications, improvements and enhancements have been proposed. The characteristics, the focus, and the simplifications of each model are different, based on the fundamental question it seeks to answer. This makes it challenging to define a unique benchmark for a scientist or a clinician that would like to make use of the existing models or propose a new one.
It is becoming clear, however, that the current trend in computational systems biology converges towards multi-scale, multi-field modeling. In the time domain, most analyses focus on short time scales, mainly on the inflammatory and proliferative phases, and only a few more recent models integrate larger times scales including the remodeling period43. In the space domain, modelers are now giving priority to multi-scale integration, from a basic mathematical representation of the cell to fully three-dimensional models of the wound site.
Here, we suggest a classification of existing models based on the aspect of wound healing they seek to address and on the mathematical tools they use to derive their governing equations, see Table 1. In this classification, the three major focus areas are: i) re-epithelialization and cell migration, ii) angiogenesis, and iii) wound contraction and collagen deposition. The four degrees of modeling complexity are i) one-dimensional or axisymmetric continuum models; ii) two-dimensional continuum models; iii) two- or three-dimensional discrete models; and iv) two-dimensional hybrid discrete/ continuum approaches.
We hope this classification alongside with the guidelines discussed above will provide new computational scientists and clinicians with the necessary framework to make use of the existing models or to create improved versions attending to the specific question they want to answer.
To better illustrate which modeling approaches have been used and how they have enhanced the understanding of wound healing from the perspective of systems biology, the following subsections present, for each major focus area defined in Table 1, the modeling and solution methods that have proved to be more effective followed by an example of how computational simulations have tested virtual hypothesis by altering the baseline model and analyzing the outcome of the computations.
The simplest mathematical models in wound healing focus on re-epithelialization and cell migration. These models typically characterize a cell population over a domain in terms of point-wise cell densities rather than considering each individual cell. Fisher's equation is the preferred representation of the population behavior. It consists of a diffusion term to capture cell migration and a logistic term to capture cell proliferation39. This abstraction is well-suited for the study of re-epithelialization, and is also considered useful for fibroblast and endothelial cell migration in wound healing39,44. During re-epithelialization, keratinocytes remain in stretch contact with one another, maintaining tissue continuity throughout the epidermis. A continuum assumption is therefore well justified. The solution to Fisher's equation is a wave propagation pattern where the advancing front of cells moves with an almost constant velocity44. Although the model is conceptually two dimensional45,46, most numerical implementations are based on the additional assumption of axisymmetry to further simplify the modeling. A more sophisticated alternative to these continuum models are discrete models, which represent each individual cell explicitly47. Discrete models typically describe the motion of each cell as a reinforced random walk in the direction of a stimulus and regulate cell division through an internal clock48.
During wound healing, keratinocytes, initially located at the edge of the lesion, crawl over the lesion to restore the protective barrier. While cell-cell cross-talk is key to this process, external signaling strongly influences the behavior of the leading edge. It overwrites the default program to initiate horizontal growth of skin, temporally bypassing the upward motion of keratinocytes from the basal membrane to the stratum corneum. In this early response of keratinocytes, TGF-β1 has been recognized as a crucial signaling pathway controlling two major events: It increases keratinocyte migration and decreases keratinocyte proliferation. A balanced action of TGF-β1 is therefore critical to avoid pathological re-epithelialization. With the help of agent-based computational models, Sun et al.49 explored the role of TGF-β1 in epidermal wound healing and identified spatio-temporal sequences of events in normal and pathological wound healing. In their baseline model, increased concentration of TGF-β1 at the edge of the wound induced a population of keratinocytes to migrate inwards and blocked their proliferation. As the leading edge moves into the lesion, the adjacent population of keratinocytes was exposed to normal levels of TGF-β1 but to an increased presence of other growth factors that stimulate cell proliferation. Normal re-epithelialization can thus be viewed as a non-proliferative, highly motile leading edge followed by a proliferative population of keratinocytes. In these computational experiments, disrupting the balance of TGF-β1 created either chronic wounds or hypertrophic scars49 .
The aspect of wound healing that has received the most attention in modeling is angiogenesis. The reason for the disproportional interest in new vessel formation is that it is not only crucial to wound healing but also to tumor growth. The basic assumptions of angiogenesis models are valid for both circumstances, and combined efforts between these two scientific communities have markedly pushed the frontiers in mathematical modeling of angiogenesis within the past decade40,50. Researchers in wound healing and tumor growth have collectively used various approaches to study the creation of new vasculature and some of their most remarkable models are extremely elaborate. They include most features of the repair sequence: inflammatory cells, one or more chemoattractants, fibroblasts, and the ECM41,51. In continuum models, new vasculature is represented as the density of endothelial cells in a point-wise fashion or as a combination of capillary tip and sprout densities52. A partial differential equation dictates the evolution of the endothelial cell density field, in analogy to Fisher's equation, but with additional coupling terms between the endothelial cells and the other fields in their immediate environment53,54. In discrete models, fractal-like approaches or cellular automata have been proposed to account for the contact between individual endothelial cells as the capillary tubes branch out55,56. Hybrid models have also successfully reproduced the angiogenesis mechanism57. The greatest value of angiogenesis models, continuous or discrete, is that they allow us to explore the effects of altered nutrient and oxygen supplies on the healing process.
During wound healing, elevated metabolic needs in the proliferative phase increase oxygen demand. In fact, the lack of oxygen or hypoxia is prevalent in the wound during the inflammatory phase and is an important cue for macrophages to release cytokines that recruit other cell types for the subsequent phases58-60. If hypoxia is severe, for example in chronic wounds, the tissue is incapable of creating new vasculature, and the healing process is significantly impaired61. Thus, oxygen therapy has been proposed to accelerate the healing process62. However, excess of oxygen can also be harmful, because it can lead to intoxication63. Hypoxia is an essential cue during the inflammatory phase and is linked indirectly to fibroblast recruitment62. The interplay of oxygen and cellular recruitment has successfully been studied with the aid of computational models64.
For example, Schugart et al. proposed a model of wound healing, in which oxygen tension across the wound is considered as an additional variable65. Their baseline model reproduced a normal healing response throughout a period of 10 days. Then, different degrees of hypoxia and hyperpoxia were simulated. Their simulations revealed that severe hypoxia cannot sustain angiogenesis, and that extreme hyperpoxia reduces the proliferation of endothelial cells, see Figure 2, top. Finally, they studied different hyperbaric oxygen therapies and concluded that 90 minutes of hyperbaric oxygen per day, enhance the healing process, see Figure 2, bottom.
Similar to angiogenesis, collagen deposition and wound contraction have mainly been studied holistically66,67. Unique features of their mathematical models are the consideration of fibroblasts and myofibroblasts68,69. Fully discrete models are not appealing to represent the ECM, and therefore modelers have either turned to hybrid or fully continuous frameworks. In the former, the chemical species and the ECM are modeled as continuous fields, whereas the cells are modeled discretely as individual entities70. The major focus of hybrid discrete/continuum models has been on collagen deposition. If we decide to also represent cell populations as continuous fields, we can select plain continuum models and more advanced mechanical theories to answer questions such as myofibroblasts-based active wound contraction71.
While the biochemical aspects of wound healing have received significant attention, mathematical modeling of the mechanical aspects of wound healing remains largely unexplored. The role of mechanical cues is currently gaining importance though since recent experimental data suggest that releasing mechanical stresses in the wound may shorten the inflammatory phase and reduce scarring8,72. Mechanical stress is transmitted across the ECM and directly affects fibroblast behavior73,74. For example, the local environment of a fibroblast can induce its transformation into a myofibroblast75. Myofibroblasts are responsible for contracting the tissue and bringing the edges of the wound together69. Although there are a few models that have incorporated wound contraction by myofibroblasts, these models have not yet been used to optimize the mechanical environment and enhance therapeutic outcomes. Nonetheless, since pathologic scarring is a major concern in wound healing, researchers have focused on modeling the flow of chemical information, which regulates collagen deposition and creation of the new ECM75.
For example, Cumming et al. studied the response of fibroblasts to TGF-β, a cytokine released by macrophages during inflammation, throughout a period of 14 days76. Using a predictive mathematical model, they were able to show that altered TGF-β kinetics in the wound have a significant impact on the healing process. They found that in their baseline model, fibroblasts tend to cluster around the zones with highest concentration of chemoattractants, where they gradually replace the fibrin clot with collagen as they reach the damaged region, see Figure 3, left. According to the simulations, reduced TGF-β diffusion causes a clustering of fibroblasts, reduced collagen synthesis, and significantly altered healing kinetics, see Figure 3, right. Another remarkable example of the use of agent based models is the study conducted by Mi et al.77 , which focuses on TGF-β1, a specific isoform of TGF-β. Altered kinetics of this chemokine are an important factor in diabetic ulcer pathology. According to their simulations, increased tumor necrosis factor-α (TNF) and decreased TGF-β1 lead to an impaired healing response of diabetic patients. The model further showed that altering these chemoattractants increased the concentration of TGF-β1, increased collagen deposition, reduced concentration of TNF-α, and reduced necrosis.
Despite intense research over the past two decades, most existing models for wound healing are still one or two-dimensional and focus on the acute rather than on the chronic response, see Table 1. The majority of these models use finite difference or finite volume methods to discretize their governing equations in space, limiting the geometries to idealized settings. An elegant way to incorporate realistic three-dimensional geometries of the wound and skin anatomy is finite element modeling. Recent trends in computational biology focus on predicting chronic soft tissue adaptation using mechanistic finite element models78,79. In an attempt to quantify how elevated mechanical stretch can alter collagen deposition and fibroblast mitosis, several groups have recently modeled chronic skin growth in response to changes in the mechanical environment80,81. These models have immediate clinical applications in skin expansion in plastic and reconstructive surgery. Predicting stress, strain, and skin area gain, skin growth models have the potential to enhance treatment for patients with birth defects, burn injuries, or breast tumor removal82.
For example, Zöllner et al. simulated skin expansion in pediatric forehead reconstruction of a one-year old girl throughout a period of 12 weeks83. Their model incorporated a thermo-mechanically consistent representation of the dermis and a phenomenological scalar field that quantified the amount of newly grown skin82. A conceptually simple and elegant abstraction of the mechanotransduction pathways and the corresponding cellular response defined the evolution of this scalar field in terms of mechanical stimuli such as tissue stretch83,84, see Figure 4. Refining this framework to incorporate the true mechanotransduction pathways during mechanical overstretch or during wound healing appears to be a next logical step towards predicting and improving of wound healing therapies in realistic three-dimensional geometries of pediatric patients.
Systems-based mathematical modeling of wound healing has achieved remarkable sophistication and is on its way to becoming an irreplaceable tool in personalized medicine. The first mathematical models of wound healing were proposed two decades ago and considered only specific aspects of the healing process. Current models include most of the key components that interact throughout the repair sequence. They have opened the floor to advanced hypotheses testing and enhanced wound management therapies. The state of the art in computational modeling of wound healing is the temporal and spatial integration of different cell types, chemoattractants, nutrients, and the ECM, interacting jointly to restore tissue integrity. Current trends in computational modeling indicate that this knowledge gained on the cellular level should be integrated in holistic multi-scale multi-fields continuum models through a bottom-up approach. The ultimate goal would be to create high-resolution system-level models with parameters calibrated at their generative level of resolution. Rather than fitting physiologically meaningless phenomenological parameters to system-level measurements, which has been the standard for the past decades, all parameters would then have a clear physiological interpretation.
Greatest attention has been paid to the biochemical features of healing such as diffusion of chemoattractants and oxygen tension. The first models have now advanced far enough to reliably predict how systematic manipulations of baseline parameters in chemical signaling can change the healing process. The biomechanical features of healing are currently receiving increasing attention. Mechanical models focus primarily on collagen deposition and active wound contraction. Yet, to date, it still remains unclear, how exactly these mechanisms are regulated when the injured skin tries to recover its homeostatic equilibrium state.
One of the major challenges in the mathematical and computational modeling of wound healing is the consideration of complex, physiological geometries in two and three dimensions. Current models have been tested primarily in axisymmetric conditions, which limits their application to a clinical setting and reduces their translational potential. Another current roadblock is the incorporation of a detailed mechanical representation of skin with temporally, spatially, and directionally varying material properties. Fortunately, recent advances in molecular biology, mechanotransduction, soft tissue mechanics, and computational modeling of soft tissues might soon allow us to bridge these gaps. The first patient-specific computational model of pediatric wound healing is likely to appear within this decade and it will constitute a major breakthrough in the progress of systems biology towards a better care for the pediatric patient.
STATEMENT OF FINANCIAL SUPPORT: This work was supported by the Consejo Nacional de Ciencia y Tecnologia (CONACyT) Fellowship and the Stanford Graduate Fellowship to Adrian Buganza Tepole; by the National Science Foundation CAREER award CMMI-0952021 and INSPIRE award 1233054; and by the National Institutes of Health grant U54 GM072970 to Ellen Kuhl.
The authors declare that they have no conflict of interest and no relationships with industry relevant to this work.