## Abstract

Wound fibrosis (i.e., excessive scar formation) is a medical problem of increasing prevalence, with poorly understood mechanistic triggers and limited therapeutic options. In this study, we employed an integrated approach that combines computational predictions with new experimental studies in mice to identify plausible mechanistic triggers of pathological scarring in skin wounds. We developed a computational model that predicts the time courses for six essential cell types, 18 essential molecular mediators, and collagen, which are involved in inflammation and proliferation during wound healing. By performing global sensitivity analyses using thousands of model-simulated wound-healing scenarios, we identified five key processes (among the 90 modeled processes) whose dysregulation may lead to pathological scarring in wounds. By modulating a subset of these key processes, we simulated fibrosis in wounds. Moreover, among the 18 modeled molecular mediators, we identified TGF-β and the matrix metalloproteinases as therapeutic targets whose modulation may reduce fibrosis. The model predicted that simultaneous modulation of TGF-β and matrix metalloproteinases would be more effective in treating excessive scarring than modulation of either therapeutic target alone. Our model was validated with previously published and newly generated experimental data, and suggested new in vivo experiments.

## Introduction

Traumatic skin injuries are often prone to exaggerated skin scar formation (1). Although all skin scars are inferior to original skin in terms of integrity and function, normally they regain the natural properties of skin after years of maturation (2). However, in many wounds, disruptions in cellular and molecular signaling lead to excessive scarring, i.e., fibrosis. Fibrotic scars (e.g., hypertrophic scars and keloids) are thicker and functionally weaker than normal scars, which makes the tissue under the scar vulnerable to subsequent injuries (3). Pathological scarring is typically associated with excessive collagen deposition and disorganized orientation of collagen fibers (4). Despite significant progress in wound care practices, pathological scarring continues to be a significant medical and economic burden (5).

Current scar management therapies include surgery, silicone sheeting, anti-inflammatory medications, laser/radiation, corticosteroids, and a multitude of topical agents (3, 6). However, the clinical efficacy of these therapies is modest (7, 8). The processes involved in scarring are complex and depend on the coordinated signaling of several molecular mediators and environmental factors (9, 10). Altered fibroblast apoptotic signaling, mechanical loading of the extracellular matrix, and altered TGF-β1 expression in fibroblasts have recently been implicated as potential initiators and propagators of wound fibrosis (11–14). Based on these findings, the targeting of mechanical signaling pathways (15) and TGF-β inhibition (16) have been investigated as possible therapeutic wound-healing strategies; however, these approaches have only had limited clinical success (17).

The absence of effective therapeutic strategies for fibrosis is due primarily to the poor understanding of the pathogenesis of excessive scarring. New approaches, complementary to traditional experimentation, may allow systematic identification of promising molecular targets and the design of optimal therapeutic interventions to treat pathological scarring. We hypothesized that pathological scarring is predominantly triggered by the dysregulation of a few key processes whose identification can be assisted by computational modeling. Computational models representing wound scarring have been developed previously (18–26). However, existing models are limited in their ability to reflect scar formation initiation and molecular signaling. The main goal of this study is 2-fold: 1) to develop a quantitative kinetic model of scar formation that describes both the inflammatory and proliferative phases of wound healing; and 2) to use the model to predict potential mechanistic drivers of pathological scarring. To achieve this, we extended our computational model of wound inflammation (27) to represent the proliferative phase of wound healing (Fig. 1). Using this extended model, we simulated normal and pathological (i.e., excessive) scarring conditions. We specifically focused on pathological scarring scenarios that involve excessive collagen accumulation.

We validated the model predictions for normal scar formation, using newly generated data from a mouse skin wound model. By analyzing 40,000 model-simulated wound-healing scenarios, we generated testable predictions regarding plausible mechanistic triggers of pathological wound scarring and promising molecular targets for its reduction. First, among the 90 molecular and cellular processes represented in the model, we identified five processes whose dysregulation may be the strongest potential triggers of pathological scarring: macrophage crowding, TGF-β degradation, fibroblast apoptosis, collagen production by fibroblasts, and collagen degradation by fibroblast enzymes. Although abnormal collagen deposition has been shown to occur in hypertrophic scars (13, 28), a systematic investigation of generally prevalent mechanistic factors affecting collagen formation (and possibly causing abnormal scarring) has not been previously undertaken. Second, among the 18 molecular mediators included in the model, we identified two classes of mediators (i.e., TGF-β and matrix metalloproteinases [MMPs]) as potential targets whose modulation may reduce excessive scarring. Third, the model predicted that simultaneous modulation of these two molecular targets may be more efficacious in reducing excessive scarring than modulation of either target alone. TGF-β and MMPs have been individually targeted to treat abnormal healing with limited success (16, 29), but combined therapies simultaneously targeting TGF-β and MMPs have not been used for treating excessive scarring. Our analysis suggests that therapeutic interventions targeting these proteins may be efficacious for treating human wounds. Our results are corroborated by existing experimental and clinical data as well as our own.

## Materials and Methods

### Animals and wound models

Ten-week-old female C57BL/6J mice purchased from The Jackson Laboratory (Bar Harbor, ME) were used. Six 3 mm full-thickness wounds were made using a biopsy punch on shaved dorsal skin under ketamine (100 mg/kg) and xylazine (5 mg/kg) anesthesia. Wound tissues harvested at designated time points were frozen at −80°C and embedded in OCT compound (Fisher Scientific, Pittsburgh, PA), or fixed in 3.7% formalin (Fisher Scientific). This animal study was approved by the University of Illinois at Chicago Institutional Animal Care and Use Committee and the U.S. Army Animal Care and Use Review Office (Fort Detrick, MD).

### Collagen measurement

Wounds from C57BL/6J mice at days 7, 14, 21, and 28 after injury were excised, weighed, and stored at −80°C until analysis. Each wound and normal skin tissue sample was hydrolyzed in 1 ml of 6 N HCl overnight at 95°C for 20 h. The hydroxyproline content was then analyzed using a hydroxyproline assay kit (QuickZyme Biosciences, Leiden, Netherlands).

### Immunofluorescence for neutrophil, macrophage, and fibroblast measurement

Frozen 8 μm sections were fixed in cold acetone and blocked with 10% goat serum for 30 min. Sections were incubated with rat anti-mouse CD68 (Abcam, Cambridge, MA), rat anti-mouse Gr-1 (BD Bioscience, San Jose, CA), and chicken anti-mouse vimentin (Abcam) for 45 min for macrophage, neutrophil, and fibroblast staining, respectively. The secondary Abs used were Alexa Fluor 594 goat anti-rat IgG, Alexa Fluor 488 goat anti-rat IgG (Invitrogen, Carlsbad, CA), and Texas Red goat anti-chicken IgY H/L (Abcam), respectively. Stained sections were observed and imaged using a digital camera. The number of neutrophils in the wound margin and wound bed per 20× field was counted. Because of overwhelming staining of CD68 and vimentin at certain time points, accurate counting of the positively stained cells was not feasible. Instead, the cell densities were quantified using ImageJ software (30).

### ELISA

Wound samples and normal skin tissue were homogenized in ice-cold PBS containing protease inhibitors for mammalian cells (Sigma Aldrich) followed by sonication. The supernatants were collected after centrifugation at 16,000 × *g* for 15 min. Protein concentrations of TNF-α, IL-6, chemokine CXCL1, CXCL2 (macrophage inflammatory protein-2), and CCL2 (MCP-1) were measured by a multiplex ELISA assay kit (ProcartaPlex; eBioscience, San Diego, CA). Protein concentrations of MMP-9 and fibronectin were determined by ELISA kits purchased from LifeSpan BioSciences (Seattle, WA).

### Computational analysis

The computational model presented in this study is an extension of our previously developed model of injury-induced local wound inflammation (27). The wound inflammation model described the kinetics of platelets, four inflammatory cell types, 11 molecular mediators, and their essential interactions during normal and chronic inflammation in response to injury. To that model, we added mathematical descriptions of the kinetics of two proliferative cell types (namely, fibroblasts and myofibroblasts), seven molecular mediators (fibronectin, basic fibroblast growth factor, MMP-1, 2, and 9, TIMP-1, and MCP-1), and three forms of collagen (tropocollagen, collagen fibrils, and collagen fibers) (Fig. 1). We modeled these components because they are widely regarded as essential cell types and molecular mediators involved in proliferation (2, 14). We modeled 90 different processes in total [e.g., cellular chemotaxis, cellular phenotype conversion, cellular apoptosis, collagen production and polymerization, molecular mediator production/degradation, and mechanical stress effects (see Supplemental Table I)] that govern the kinetics of both the inflammatory and proliferative responses via 108 model parameters (Supplemental Table I). The default parameter set represents the cellular and molecular kinetics during a normal scarring scenario. Our model is a coupled system of 27 ordinary differential equations and one delay differential equation (Supplemental Table II). We performed all computations in the software suite MATLAB R2015b (MathWorks, Natick, MA) and solved the model equations using the MATLAB solver DDE23 with default tolerance levels. We computed time courses for each of the 28 model variables for 40 d after wounding. We additionally considered three variables representing the total concentrations of neutrophils, macrophages, and collagen. We calculated the values for these variables from the model variables representing two neutrophil phenotypes, two macrophage phenotypes, and three collagen forms. Thus, we computed the time courses for 31 variables in total. The MATLAB code for our computational analyses is provided as supplemental data in the online version of this article.

### Goodness of fit

To assess the goodness of fit of our model with respect to our experimental data, we calculated the root mean squared errors (RMSEs) between the model predictions and the data shown in Figs. 2 and 3. Additionally, we quantified the likelihood that our model predictions would fall within the ranges of our corresponding experimental data. To this end, we first used our experimental data to calculate a tolerance interval (TI) for each measured cell type and molecular mediator at each time point. For a variable *Y*, we defined the TI to be the interval whose midpoint is the mean value, *Ŷ*, of *Y* and whose width, *w*, is evaluated by the formula given below (assuming that *Y* is normally distributed, which was true for the vast majority of our experimental measurements):where σ is the SD of *Y* and *k* is a factor that determines the width of the TI. The value of *k* is calculated as described previously (31); it depends on the specified fraction of the population (we used a value of 95% in our calculations) that falls within the interval, the specified confidence level (we used a value of 0.05 in our calculations), and the sample size. Therefore, the TI is expected to contain 95% of the values of a random variable. Next, for each measured (and model-predicted) molecular mediator or cell type, we calculated the percentage of the model-predicted time-course points that fell within the TIs of the corresponding experimental measurements (i.e., the measurements at the corresponding time point, when such measurements were available in our data set). Thus, the higher this percentage, the greater the likelihood that a model prediction for a given model variable will fall within the same range as the experimental measurements. This is because a new measurement at a given time point is expected, with a 95% probability, to lie within the TI of the data points measured at the same time point.

### Sensitivity analysis

First, we performed a local sensitivity analysis to assess the model’s robustness and remove any non-essential interactions, as previously described (27). In this analysis, the model parameters were varied near their default values (±1%). Next, we performed two distinct types of global sensitivity analysis (GSA): partial rank correlation coefficient (PRCC) analysis and extended Fourier amplitude sensitivity testing (eFAST) analysis. For these analyses, we used a slightly modified version of the MATLAB code provided in a publication by Marino et al. (32). In both PRCC and eFAST analyses, we analyzed 40,000 simulated wound-healing scenarios, wherein the model parameter values were randomly selected from a 4-fold interval (2-fold in each direction) around their default values. PRCCs provide a measure of the strength of monotonic dependence between a model parameter and a model variable, whereas the eFAST sensitivities reflect the strength of propagation of model parameter variation to the model variables. Although these two methods measure different quantities, they both provide a measure of the strength of influence of the model parameters on the model variables, when the model parameters are varied simultaneously.

### PRCC

For this analysis, we used Latin hypercube sampling to create 40,000 unique parameter sets (33). We accomplished this using the MATLAB function LHSDESIGN. Next, we simulated the time courses for each model variable in 40,000 wound-healing scenarios using the 40,000 generated parameter sets. We then calculated PRCCs (with their associated *p* values) between each of the 31 model variables and each of the 108 model parameters at 40 simulation time points, each representing 1 d after wounding. PRCC values varied between −1 and +1. Higher absolute values of the PRCCs indicated a stronger influence of a particular model parameter on a specific model variable. A PRCC with a *p* value <0.01 indicated that it was significantly different from zero. The sign of the PRCC values indicated the positive or negative directionality of the correlation between model parameter and a model variable. At the end of our PRCC analysis, we obtained 108 PRCCs for each model variable.

### eFAST

In this analysis, the model parameters were varied sinusoidally at different frequencies, encoding the identity of parameters in the frequency of the parameters’ variation. We performed ∼40,000 simulations, in which each model parameter was sampled 75 times and resampled five times for a given frequency of variation. Once the time courses were generated by solving the model equations (Supplemental Table II) for the generated parameter sets, the variance in each model variable was partitioned to determine the fraction of variance contributed by each model parameter. Fourier analysis was used to calculate the strength of each parameter’s frequency in the variations of the model’s variables. The eFAST sensitivities obtained at the end of the analysis reflect the strength of propagation of a model parameter’s variation to a given model variable. The values of these sensitivities varied between 0 and 1. For each of the model’s 31 variables, we obtained 108 individual eFAST sensitivities.

### Identification of strongly influential processes

To identify plausible mechanistic triggers for pathological scarring, we focused on the PRCC and eFAST analysis results for two specific model variables, namely, the total collagen concentration and the fibroblast concentration. For both variables, we divided the respective 108 PRCCs and 108 eFAST sensitivities (one PRCC and one eFAST sensitivity per each model parameter; the PRCCs and eFAST sensitivities were analyzed separately) into three groups using a *k*-means clustering algorithm (34). We used the MATLAB function KMEANS to perform the clustering. Group 1 for PRCCs contained the highest absolute PRCCs and group 1 for eFAST sensitivities contained the highest eFAST sensitivities. We used the minimal values in group 1 from both the PRCC and eFAST clustering analysis as a threshold for identifying the model parameters with the strongest influence on scarring. For a given model variable (i.e., the total collagen concentration or the fibroblast concentration), the model parameters for which the PRCC and eFAST sensitivities were both above their respective threshold values and the *p* values associated with the PRCCs were statistically significant (i.e., *p* < 0.01) were regarded as the parameters exerting the strongest influence on that model variable. The dysregulation of the biological processes represented by these most influential model parameters was regarded as a plausible mechanistic trigger of pathological scarring.

### Estimation of mechanical stress effect functions

We modeled the effect of mechanical stress on the production of certain molecular mediators and on the production of tropocollagen, as described below. First, we assumed that the mechanical stress, *M*_{stress}, generated by the extracellular matrix (ECM) would be directly proportional to the concentration of collagen fibers present in the wound, [*Coll _{fib}*]. We calculated the coefficient of proportionality by performing linear regression of the experimental data from Roeder et al. [see Figure 10 and Table II in Ref. (35)]. The linear dependence of mechanical stress on collagen concentrations between 0.3 and 3.0 mg·ml

^{−1}was defined by the following equation:

*M*

_{stress}= 3 × 10

^{9}[

*Coll*]. Next, we calculated the strain (i.e., percentage of displacement in the ECM, designated as

_{fib}*M*

_{strain}) induced by a given level of mechanical stress,

*M*

_{stress}. We performed linear regression of the stress-strain curve for collagen type I from Aarabi et al. [see Figure 1B in Ref. (13)]. The linear dependence of the induced strain on applied mechanical stress was described by the following equation:

*M*

_{strain}= 0.203

*M*

_{stress}. This allowed us to calculate the strain induced in the ECM as a function of the mechanical stress applied by the collagen fibers in the ECM. Based on the level of strain in the ECM, the production of molecular mediators and collagen by fibroblasts is altered. The calculated strain value (i.e.,

*M*

_{strain}) was used as the input to five different functions describing the effect of the applied strain on the rates of IL-6, CXCL8, TGF-β, and tropocollagen production by fibroblasts, as well as on the fibroblast proliferation rate (see Supplemental Table II). These functions are shown in Supplemental Table I.

## Results

### Prediction and validation of normal inflammation and scarring time courses

We measured the levels of different inflammatory and proliferative cell types and molecular mediators in excisional wounds of wild-type C57BL/6J mice (Fig. 1). We then used these data, as well as published data from other animal wound models, to validate our extended computational model (developed independently of these validation data sets; see Supplemental Tables I, II). To establish that the model predictions are not exclusive to just one specific mouse strain (i.e., C57BL/6J mice), we compared our model predictions with previously published data on traumatic skin injuries in BALB/c mice from six independent studies (Fig. 2, black symbols and black lines). Moreover, to rule out the possibility of gender bias in our data, we compared the time courses of four major inflammatory cytokines (TNF-α, IL-1β, IL-6, and CXCL1) from our data (generated using female C57BL/6J mice) with those of previously published wound data in both male and female C57BL/6J mice (Supplemental Fig. 2). In these comparisons (Fig. 2, Supplemental Fig. 2), the model-predicted time courses for the inflammatory cells and molecular mediators showed reasonably good agreement with experimental data from BALB/c and C57BL/6J mouse wounds. Interestingly, the inflammation kinetics in our extended computational model (which represented both the inflammation and proliferation phases of wound healing) was only marginally different from the predictions of our previously developed, experimentally validated computational model that represented only the inflammation phase [see Figure 3 in Ref. (27)]. Yet, the direct comparison of the extended model predictions with the wound inflammation data in Fig. 2 was necessary, because our original model was validated using inflammation data obtained for conditions different from traumatic skin injuries (36, 37).

In our model predictions for the proliferative phase, we detected a peak in fibroblast and myofibroblast concentrations on days ∼6–7 postwounding, respectively (Fig. 3A, 3B). These predictions matched our cell density measurements in the mouse, as well as those from guinea pig (38), and pig skin (39) wound models (Fig. 3A, 3B). In our experiments, the total collagen concentration (which represents the sum of the tropocollagen, collagen fibril, and collagen fiber concentrations) peaked around day 7 (Fig. 3C). Our model predictions for the total collagen concentration matched these data and those from pig wounds (39) (Fig. 3C). Moreover, our model predictions regarding the kinetics of key proliferative mediators, such as fibronectin, MMP-9, and TIMP-1, showed reasonable qualitative agreement with available experimental data from different animal wound models (Fig. 3D, 3F).

Quantitative comparisons were hampered by differences between the units used to report experimental measurements and the units in our model predictions. Nonetheless, to test the goodness of fit of our model, we calculated the RMSEs between the normalized model predictions and our normalized experimental data for 13 model variables representing different inflammatory cells and molecular mediators. The RMSEs for the model variables ranged from 0.28 to 0.48 (for individual RMSE values, see the legends for Figs. 2, 3). Moreover, we computed the percentages of model-predicted values that fell within the TIs of the corresponding experimental data points. With the exception of three model variables (namely, the total macrophage concentration and the concentrations of fibronectin and MMP-9), more than 70% of the model predictions fell within the TIs of their respective experimental data points. Interestingly, for three model variables (i.e., fibroblast, TNF-α, and MCP-1 concentrations), 100% of the model-predicted values fell within the TIs (for individual likelihood percentage values, see the legends for Figs. 2, 3). Overall, the differences between our model predictions and our experimental data were comparable to the differences between distinct experimental data sets characterizing the same cell type or molecular mediator. This indicates that our computational model captured typical experimentally observed, injury-induced inflammatory and proliferative responses.

### Identification of plausible mechanistic triggers of pathological scarring

We identified plausible mechanistic factors triggering pathological scarring by applying two distinct GSA methods commonly used in computational biology ‒ the calculation of PRCCs and eFAST sensitivities (see *Materials and Methods* section). For each model variable, we calculated 108 PRCCs (Supplemental Fig. 1) along with their associated *p* values, and 108 eFAST sensitivities corresponding to the 108 main parameters in our model. Both the PRCCs and the eFAST sensitivities reflect the strength of influence that a given model parameter exerts on a specific model variable. Higher absolute values of statistically significant PRCCs (i.e., *p* < 0.01) and eFAST sensitivities indicate a stronger influence, but because PRCCs and eFAST sensitivities are calculated differently, the two GSA methods complement each other. For this analysis, we focused on the model variables representing total collagen and fibroblast concentrations in the wound, because they are often altered in pathological scarring (13, 40). On the last day reflected in our model predictions (i.e., day 40 postwounding, which shows the final scarring outcome in our study), the majority of the modeled processes did not show a significant influence on the total collagen (Fig. 4A) or the fibroblast concentration (Fig. 4B) in both the PRCC and the eFAST analyses (Fig. 4). Note that the bars representing the eFAST sensitivities are superimposed on the bars representing the PRCCs. After thresholding the PRCCs and the eFAST sensitivities (Fig. 4, horizontal lines; see also *Materials and Methods* section), we identified the rate of collagen production by fibroblasts, the rate of collagen degradation by fibroblast enzymes, and fibroblast apoptosis rate as the strongest modulators of the total collagen concentration. The strongest modulator of the fibroblast concentration was the fibroblast apoptosis rate.

Using both the PRCC and eFAST results, we identified the most influential processes for the collagen and fibroblast variables in our computational model across the simulated 40 d postinjury. Interestingly, macrophage crowding showed the highest influence on the fibroblast concentration at earlier time points (days 1–5 postinjury) (Fig. 5A, 5B), whereas after day 5 this concentration was controlled primarily by the fibroblast apoptosis rate (Fig. 5A, 5B). Similar to the fibroblasts, the total collagen was strongly influenced by macrophage crowding between day 1 and day 5 (Fig. 5C, 5D). Moreover, the TGF-β degradation rate appeared to strongly influence collagen on day 1 through day 10 (Fig. 5C, 5D). During the proliferative phase (i.e., day 10 through day 40), collagen was strongly influenced by the fibroblast-dependent collagen production rate, fibroblast apoptosis rate, and the rate of collagen degradation by fibroblast-released enzymes (Fig. 5C, 5D). Importantly, although our computational model was developed using in vitro data, two of the five processes identified as the most influential (i.e., fibroblast apoptosis and TGF-β signaling) had earlier been recognized as potential initiators of pathological scarring in vivo (28, 41). To sum up, our GSA predicted five processes whose dysregulation may trigger pathological scarring: macrophage crowding, fibroblast apoptosis, collagen production by fibroblasts, collagen degradation by fibroblast enzymes, and TGF-β degradation.

### Restoration of normal collagen level during pathological scarring

Considering the evidence that fibroblast apoptosis and fibroblast collagen expression are altered in hypertrophic scars (13, 40), we modified these two model parameters to generate predictions for pathological scarring (Fig. 6). Although the shape of the model-predicted time courses for total collagen and fibroblasts during pathological scarring (data not shown) was qualitatively similar to their normal scarring kinetics, the collagen concentrations during model-predicted pathological scarring were ∼2-fold higher (Fig. 6A), which is consistent with experimental data (42–44).

We used this strategy to investigate whether the modulation of TGF-β degradation and the collagen degradation by fibroblast enzymes may be promising intervention strategies to reduce pathological scarring. We simulated three intervention scenarios for pathological scarring. In the first intervention scenario, we increased the TGF-β degradation rate. This scenario represents the removal of TGF-β from the wound or the introduction of a TGF-β inhibitor into the wound. In the second intervention scenario, we increased the fibroblast-mediated collagen degradation rate. This scenario represents the addition of collagen-degrading fibroblast-derived enzymes, such as MMPs. In the third intervention scenario, we combined the first two scenarios. We looked at the intervention effects on the time course of total collagen. Although the strategies involving either TGF-β degradation modulation or enzymatic collagen degradation modulation were each successful in reducing the pathological-scarring collagen levels (Fig. 6A), the best restoration of collagen levels to a normal scarring scenario was achieved when the two intervention strategies were implemented simultaneously (Fig. 6A).

We validated our model predictions regarding the TGF-β–targeting intervention (Fig. 6B) by comparing them with experimental data from the mouse model of bleomycin-induced scleroderma (45). By comparing our model predictions for total collagen concentration on day 28 and day 42 with these experimental data (Fig. 6B), we showed that our approach to modeling TGF-β degradation–rate modulation could capture the effect of TGF-β inhibition reasonably well. Moreover, this analysis showed that TGF-β and collagen-degrading enzymes, such as MMPs, could be potential molecular targets whose combined regulation may reduce pathological scarring in wounds.

## Discussion

The multitude of complex interactions during wound healing obfuscates the prognosis and treatment of wound-healing pathologies. Our computational model of scarring could successfully capture temporal dynamics of cell types, molecular signaling mediators, and collagen during normal scarring, which was validated by comparisons with experimental data (Figs. 2, 3). Among the 90 modeled processes, we identified five key processes (macrophage crowding, fibroblast apoptosis, collagen production by fibroblasts, collagen degradation by fibroblast-derived enzymes, and TGF-β degradation) whose dysregulation plausibly triggers pathological scarring. Moreover, we predicted potential efficacious interventions to reduce pathological scarring, which involved individual and simultaneous modulation of two classes of molecular targets (i.e., TGF-β and MMPs).

Hypertrophic scars and keloids are conceivably the most widely encountered forms of pathological scarring in skin wounds (40). These scars have a markedly different collagen content and collagen architecture compared with normal scars (40, 42, 44, 46). Recently, abnormal fibroblast apoptotic signaling has been recognized as a potential initiator of fibrosis (28). In fact, hypertrophic scar–derived fibroblasts exhibit disrupted apoptotic signaling (13, 28, 47) and upregulated collagen and TGF-β expression (41), which is consistent with our model predictions regarding fibrosis triggers. Although the data in experimental studies are, by necessity, specific to the chosen conditions and animals, our model predictions characterized the effects of fibroblast apoptosis and collagen production modulation in 40,000 simulated wound-healing scenarios, which accounted for intersubject variability and variations in wound initiation. Our analysis thus suggests that the triggering of pathological scarring by altered fibroblast apoptosis signaling and increased collagen production may be widely spread and, in fact, typical for mammals.

Among the molecular mediators involved in wound healing, TGF-β has the broadest spectrum of effects during all wound-healing phases (48). Because of its pleiotropic nature, targeting TGF-β to treat pathological scarring has naturally been the focus of extensive investigations (16, 17, 49). Recent approaches to antagonize TGF-β–stimulated fibrosis include the use of TGF-β–neutralizing Abs, an antagonist of TGF-β activation, and TGF-β receptor blockers (4, 47, 49). In accordance with these results, our computational analysis indicated that, among the 18 modeled mediators, TGF-β may be the most promising target for reducing excessive scarring. Although different TGF-β inhibitors have shown varying degrees of success in animal and preclinical studies, the major problem lies in their modest performance in clinical trials (17). It is conceivable that the identification of viable molecular targets needs to be reinforced by the identification of optimal strategies for their targeted regulation, which include determining the optimal intensity and timing (41). This notion, combined with our results, highlights the necessity of such investigations for TGF-β inhibition. Importantly, comprehensive searches for optimal dosage and intervention time points for therapeutic strategies can efficiently be performed using computational models (50, 51).

Besides TGF-β, we identified the collagen degradation by fibroblast-produced enzymes as a strongly influential process for collagen regulation (Figs. 5, 6). Collagen remodeling in wounds is primarily regulated by the levels of MMPs and TIMPs (52, 53). In fact, scarless wounds have a higher ratio of MMP to TIMP levels, promoting collagen remodeling and reducing collagen accumulation (54). Although the effect of MMPs on collagen remodeling is known, it has not been explored thoroughly as a therapeutic modality for pathological scar treatment. One of the reasons may be the known propensity of MMPs to trigger wound chronicity (55). Our model predictions suggest that MMP modulation during the later phases of healing can strongly influence collagen accumulation (Figs. 5A, 5B, 6A). Therefore, a thorough evaluation of MMP regulation effects can determine safe dosage ranges and optimal intervention time points for therapies involving MMP administration to treat pathological scarring. Moreover, our results suggest that simultaneous modulation of both TGF-β and MMPs may provide a more efficacious means to control fibrosis compared with modulating these molecular components individually. This finding, which is consistent with the rapidly evolving paradigm centered on the use of combinatorial therapies for proliferative diseases (including cancer) (56), warrants further investigations into combinatorial therapeutic approaches for wound healing.

Our computational approach has several limitations arising from the simplifying assumptions required to describe the complex wound-healing biology. First, we assumed that the stress and the corresponding strain experienced by the proliferative cells depend solely on the wound collagen content. Second, we describe the polymerization of tropocollagen into collagen fibrils and the polymerization of the fibrils into collagen fibers using two constant rates. In reality, these rates may vary depending on the tropocollagen and collagen fibril concentrations (57). Lastly, a fibrotic wound may have different mechanistic triggers depending on a wide variety of factors, including the wound origin (e.g., burn wounds, ulcers, etc.), size, and infection level. In fact, even the load of commensal skin bacteria in a wound is known to affect its healing (58). As another example, neutrophil extracellular traps (which attack microbial agents in the wound) and their release from apoptotic neutrophils (i.e., NETosis) have recently been shown to influence the signaling of different wound macrophage phenotypes, thereby impairing wound healing (59, 60). Because our work is focused on injury-induced (rather than infection-induced) inflammation and sufficient mechanistic data to model processes such as NETosis are not yet available, we do not account for such additional layers of complexity in our (rather parsimonious) computational model of wound healing.

Computational models offer a systematic way to study complex biological systems. Not only can they provide mechanistic insights (61, 62), but they can also act as surrogate systems to design therapeutic approaches and develop diagnostic tools for various pathological conditions (63, 64). Our computational results illustrate the utility of integrated computational/experimental approaches in the studies of wound fibrosis.

## Disclosures

The authors have no financial conflicts of interest.

## Acknowledgments

The authors are grateful to Dr. Sridhar Ramakrishnan for comments regarding the statistical analysis presented in the paper.

## Footnotes

This work was supported by the Clinical and Rehabilitative Medicine Research Program of the U.S. Army Medical Research and Materiel Command, Fort Detrick, MD.

The opinions and assertions contained in this study are the private views of the authors and are not to be construed as official or as reflecting the views of the U.S. Army or of the U.S. Department of Defense.

The online version of this article contains supplemental material.

Abbreviations used in this article:

- ECM
- extracellular matrix
- eFAST
- extended Fourier amplitude sensitivity testing
- GSA
- global sensitivity analysis
- MMP
- matrix metalloproteinase
- PRCC
- partial rank correlation coefficient
- RMSE
- root mean squared error
- TI
- tolerance interval.

- Received July 22, 2016.
- Accepted November 15, 2016.

- Copyright © 2017 by The American Association of Immunologists, Inc.