|Home | About | Journals | Submit | Contact Us | Français|
The pathogenesis of vocal fold scarring is complex and remains to be deciphered. The current study is part of research endeavors aimed at applying systems biology approaches to address the complex biological processes involved in the pathogenesis of vocal fold scarring and other lesions affecting the larynx.
We developed a computational agent-based model (ABM) to quantitatively characterize multiple cellular and molecular interactions involved in inflammation and healing in vocal fold mucosa after surgical trauma. The ABM was calibrated with empirical data on inflammatory mediators (eg, tumor necrosis factor) and extracellular matrix components (eg, hyaluronan) from published studies on surgical vocal fold injury in the rat population.
The simulation results reproduced and predicted trajectories seen in the empirical data from the animals. Moreover, the ABM studies suggested that hyaluronan fragments might be the clinical surrogate of tissue damage, a key variable that in these simulations both is enhanced by and further induces inflammation.
A relatively simple ABM such as the one reported in this study can provide new understanding of laryngeal wound healing and generate working hypotheses for further wet-lab studies.
Treatment outcomes for patients with vocal fold scarring remain grim, despite remedial efforts that have been undertaken to date.1–3 A possible factor in this poor prognosis is the wide individual variability displayed by vocal fold scarring,4 which hampers efforts to predict disease progression and treatment response across patients.5 Surgically induced trauma, like all other forms of trauma, causes a highly complex inflammation and healing process that is modulated by numerous cells and their products such as cytokines, growth factors, and extracellular matrix (ECM) substances. The complexity of the disorder complicates the decision-making task for voice management professionals, and the task becomes increasingly daunting in the context of an era emphasizing personalized, preemptive, and predictive medicine.6 Systems biology is a fast-advancing field that has significant potential utility to address this challenge. A systems-oriented (instead of reductionist) approach may capture the complex dynamics in disease progression and treatment response by using an integrated in vitro/in vivo/in silico (computational) approach.7 The present report is part of a series of systems biology–driven studies aimed at developing a computational platform with which to aid researchers and clinicians in 1) investigating the complex processes involved in the pathogenesis of vocal fold injury and inflammation and 2) testing the effects of behavioral and pharmaceutical treatments on stressed or traumatized vocal folds.8–10 In this article, we report the first generation of a computational model of vocal fold surgical trauma that was primarily built on published animal data from experimental vocal fold injury.
Extensive experimental and clinical data on vocal fold scarring in animals have been accumulated across several laboratories.2,4,11–23 These studies have undeniably advanced our understanding of the microscopic and molecular characteristics of the disease. However, only a small fraction of the system behavior in vocal fold inflammation and healing has to date been represented in a clinically useful way. The aim of the present study was to integrate some of these time-varying data via a mechanistic computational model, in order to construct a general picture of the pathogenesis of vocal fold scarring. The framework of choice for modeling the processes of injury, inflammation, and repair was the agent-based model (ABM).24–26 This class of models represents a relatively new approach for defining macroscopic emergent properties via “bottom-up” simulation of molecular and cellular processes. This type of “bottom-up” modeling represents a complex system as a collection of entities called agents. Governed by a set of predetermined rules, each agent can individually execute a series of rule-based operations. Note that an agent can represent cell(s), protein(s), or gene(s) as an entity. The rules can involve mathematical equations or “if…then” conditional statements. The relative importance of various rules is dictated by the values of model parameters. This type of model is unique, because it can produce stochastic behavior, which may account for the issue of variability in population dynamics as observed in the real world. Also, the programming languages used to create an ABM are relatively intuitive and concrete. The biological behavior identified in basic science is easier to translate into the rules in an ABM than are the mathematical equations in equation-based modeling. Agent-based modeling has been applied for human acute phonotrauma with satisfactory simulation accuracy.9 However, the magnitude of injury in phonotrauma is remarkably smaller than the injury in surgical trauma, and thus, differentiated tissue responses would be expected between phonotrauma and surgical trauma. In the current study, the human phonotrauma model was modified and recalibrated to specify the model to the surgical injury of interest.
In the present study, an ABM simulating the response to surgical trauma in animals was designed to augment a preexisting human phonotrauma ABM,9 by use of published rat messenger RNA (mRNA) data. Studies of rat vocal folds were used for model calibration and validation because these data were the most comprehensive among the animal species in terms of 1) the wide spread of time points following injury and 2) the relatively complete profiles of the changes in inflammatory mediators and ECM substances following injury.11,14,27–29
Original empirical mRNA tissue data that provided the basis for 2 published articles on rat vocal fold injury11,27 were used for model calibration and validation in this experiment (see below; individual data points were graciously provided by the authors). The animal surgical protocols were identical in these 2 studies. In brief, Sprague-Dawley male rats (4 to 6 months old) were used, and injuries were induced with a 25-gauge needle and microforceps to strip the vocal folds until the thyroarytenoid muscle was exposed. All laryngeal specimens were harvested and stored in the same manner after injury. Real-time reverse transcription–polymerase chain reactions were used to measure in vivo mRNA for the expression of inflammatory mediators and ECM substances.
Messenger RNA levels were expressed as the ratio of the concentration of target gene to that of the housekeeping gene β-2 microglobulin in a natural logarithmic (ln) scale. Mathematically, the ln scale can only be defined for positive real numbers or nonzero complex numbers. However, from the practical consideration of modeling, we could not exclude the case that zero values would be predicted by the ABM, ie, that no mRNA expression would be present for a particular marker. In that case, an error output would be returned if an ln scale was used in the model. Accordingly, nontransformed data were used for modeling purposes.
Next, data were inspected by use of the SPSS 15.0 statistical program (SPSS Inc, Chicago, Illinois) for each marker at each time point. Individual data showing more than 3 times the interquartile range (ie, the difference between the 75th percentile and the 25th percentile) were regarded as “extremes” and were excluded from the data pool for subsequent model calibration and validation.
The freeware Netlogo 4.0.3 (Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, Illinois) was used as the platform for model building and simulation. The previously reported ABM of human phonotrauma9 was expanded and modified for the purpose of simulating vocal fold surgical trauma in this study. Originally, the model9 was specified to the setting of phonotraumatic injury through iterative verification and calibration procedures with vocal fold secretion data from a human phonotrauma study (Verdolini Abbott et al, unpublished observations). The model was composed of 1) platelets; 2) cells (neutrophils, macrophages, and fibroblasts); 3) a growth factor (transforming growth factor [TGF]–β1) and 3 cytokines (interleukin [IL]-1β, tumor necrosis factor [TNF]–α, IL-10) involved in inflammation and wound healing; 4) a matrix substance (collagen type I); and 5) a tissue damage function analogous to alarm or danger signals constituting positive feedback to induce additional inflammation.30 In the present study, additional mediators, growth factors, and ECM substances were added to the previous phonotrauma ABM.9 Specifically, the new model variables were 2 cytokines (IL-6 and IL-8), a collagenase (matrix metalloproteinase [MMP]–8), a growth factor (basic fibroblast growth factor [bFGF]), and 2 ECM components (elastin and HA). MMP-8 and bFGF have been reported for their presence in vocal folds and their effects on vocal fold ECM components, such as collagen, elastin, and HA.31–36 Also, IL-6 and IL-8 have been suggested as active mediators in inflammation and repair within the respiratory system.37–53 Thus, the expansion of the model was intended to yield a more complete simulation of the interplay among cells, inflammatory mediators, growth factors, and ECM substances in the progression of vocal fold healing following acute injury. In the current ABM framework, agent variables were used to represent platelets, cells, and ECM, whereas patch variables were used to represent inflammatory mediators. Because of space limitations, a fuller description of ABMs is not given here; it can be found elsewhere.8,9
The model was then specified for the animal surgical trauma setting by calibrating the values of model parameters by use of the respective empirical data. Detailed literature on inflammation and healing was reviewed to identify rules for the further development of the ABMs in the present context.54–77 Further, relevant literature on the vocal folds17,31,34 was used to adapt the model to the setting of vocal fold injury. The cell source and biological functions of the existing and augmented (in italics) models are summarized in Table 1.
The overall flowchart of the new ABM is displayed in Fig 1. In the model, the assumption was that 1 step of simulated time represents 0.5 hours. The changes in temporal concentration of platelets, cells, mediators, tissue damage, and ECMs were plotted and refolded into the model at each time step.
In the initial setting of the model, circulating neutrophils, macrophages, and platelets were in the blood capillary region, whereas resident macrophages and fibroblasts were present with a random distribution within the tissue region. An initial mucosal injury (with magnitude set to an arbitrary value of 40, calibrated for the mean inflammatory response to surgical trauma; see below) was simulated as an initial traumatic stimulus to the mucosal tissue in the middle of the region, triggering platelet degranulation and native matrix substance degradation. Shortly afterward, a chemoattractant gradient was created that stimulated the infiltration and activation of neutrophils and macrophages from blood capillaries. Shortly thereafter, resident macrophages and fibroblasts were activated by mediators and “alarm or danger signals” (also known as damage-associated molecular patterns [DAMPs]) composed of ECM fragments.78,79 Activated neutrophils and macrophages secreted collagenase (MMP-8) and a proinflammatory mediator (TNF-α) that further degraded the matrices, causing secondary tissue damage. Fibroblasts secreted ECM molecules to repair both the initial and the inflammation-induced damage. The complete set of ABM rules governing these variables is detailed in a dissertation by Li.8
Pattern-oriented analysis80,81 was used to estimate the model’s parameter values through an iterative model calibration process. By this approach, the patterns of simulation-generated data curves were compared with 1) the patterns of wound healing reported in the wound healing literature across a roughly 2-week period (qualitative model calibration) and 2) the empirical data of inflammatory mediators and matrix substance in laryngeal tissue from rat vocal fold studies following surgical injury (quantitative model calibration). If the model-predicted and empirical curves failed to match, the model would be calibrated to minimize differences between the curves. Of note, not the structure of the model (ie, components and rules), but only the values of the parameters, were adjusted during the calibration process.
First, a qualitative parameter estimation was carried out to test whether the model reproduced the generally accepted patterns of cellular and molecular responses according to the literature on surgical skin wound healing,57,62,65,67,82–84 as well as on surgical vocal fold wound healing14,28,29 (Table 2). The user-defined initial magnitude of mucosal injury was first set at a value of 40 (range, 0 to 40 in arbitrary units of damage), because that setting represented realistic predictions of massive mucosal damage and healing when compared with the general consensus around surgical wound healing documented in the literature. The pretraumatic values of inflammatory markers (IL-1β, IL-6, IL-8, IL-10, TNF-α, and MMP-8) were set to zero. Simulations were then carried out to determine model outcomes and compare the model’s outputs with prespecified patterns reported in the skin and vocal fold wound healing literature14,28,29,57,62,65,67,82–84 (Table 2). If the model-generated and empirical curves failed to match according to the gross qualitative criteria for that literature, the model’s parameters were adjusted iteratively to produce a better qualitative match.
When the qualitative behavior of the simulation appeared satisfactory, quantitative calibration of the model was carried out by comparing the model outputs with specific experimental data.11,27 First, the model’s behavior was calibrated by adjusting parameter values to fit the quantity and time course of measured mRNA levels of vocal fold mediators and ECM products in the surgery-traumatized tissue.11,27 Simulations were then carried out for the rat population up to 16 hours after surgical injury by inputting the average baseline mRNA levels of mediators (IL-1β, TNF-α, and TGF-β1) and matrices (procollagen type I, elastin synthase, a surrogate marker for elastin, and hyaluronic acid synthase [HAS]-2, a surrogate marker for hyaluronan) in rat laryngeal tissue11,27 and then adding a surgical trauma event.
The initial magnitude of mucosal injury, which denoted the surgical trauma event, was again set at a value of 40 as in the qualitative calibration. Simulation outputs for each inflammatory marker (IL-1β, TNF-α, and TGF-β1) and matrix marker (procollagen type I, elastin synthase, and HAS-2) in laryngeal tissue for the rat population were compared with the empirical data across 4 time points: 1 hour, 4 hours, 8 hours, and 16 hours following surgery.11,27 The model parameter values were iteratively adjusted to achieve optimal fit to the empirical laryngeal mRNA data. The quantitative model calibration iterative process was continued until the model eventually yielded a satisfactory match between the simulation data and the empirical data.
After the quantitative model calibration, the calibrated ABM was tested for its accuracy in predicting mRNA levels of mediators and matrices at the 24-hour and 72-hour time points.11 After manual input of the population’s baseline levels of IL-1β, TNF-α, TGF-β1, procollagen type I, elastin synthase, and HAS-2, a surgical trauma event was added. (The initial magnitude of mucosal injury was 40.) The ABM simulation was run 100 times for up to 5 simulated days after surgery, in order to generate a representative data pool for subsequent statistical analysis (Kevin Kim, personal communication, 2009). Subsequently, the ABM was statistically evaluated (SPSS 15.0 statistical software) by comparing the predicted levels of each inflammatory marker and ECM marker with the corresponding marker levels at 24 hours and 72 hours for the rat population as a whole. The α level was set as 0.05. Because multiple comparisons of 6 markers were involved, a Bonferroni adjustment (α/6) was used to control for the family-wise type I error. Thus, a 99% confidence interval was computed for each marker; ie, 6 confidence intervals were computed in total (for the 6 markers IL-1β, TNF-α, TGF-β1, procollagen type I, elastin synthase, and HAS-2) at the 24-hour and 72-hour time points from the simulation runs. If the empirical result for a given marker fell within the 99% confidence interval of the simulation runs, the model was considered adequate to predict the levels of markers seen in the empirical experiment.
Although our previously described ABM was calibrated and validated against protein levels of inflammatory mediators in laryngeal secretions from human subjects, the current ABM of surgery-induced inflammation was calibrated against mRNA levels. We hypothesized that this approach was appropriate because prior studies85,86 generally suggest a fairly good temporal and quantitative correspondence between mRNA and protein levels for the ECM markers included in the ABM. The surgical injury ABM generally reproduced and predicted population-trend inflammatory mediator and ECM marker trajectories in the rat population following trauma (Fig 2). Specifically, the ABM predicted 24-hour and 72-hour mediator values in 83% of cases (ie, 10 of 12 cases were within the 99% confidence interval).
Both empirical and simulation results showed time-dependent changes in mRNA gene expression of mediators and ECM markers from rat vocal folds after surgical injury. The IL-1β expression was elevated to a great degree immediately after injury, remained elevated up to 8 hours later, decreased notably by 16 hours after injury, and remained at a low level until the end of the simulation (day 5 of simulated time). The TNF-α expression showed relative fluctuation across the simulation period. The TGF-β1 expression increased progressively after injury up to day 3, and was predicted to drop progressively until the end of the simulation.
The ECM marker elastin synthase showed early expression immediately after injury and persisted until the end of simulation. The model further predicted that the expression of the ECM markers procollagen type I and HAS-2 would be elevated starting on day 1 after injury and would peak at day 3.
The model predicted a large inflammatory response in the form of massive inflammatory cell infiltration triggered by the high magnitude of tissue damage induced by surgical trauma. In particular, extensive platelet and neutrophil infiltration was predicted for the initial 12 hours after injury (Fig 3A). Macrophages and fibroblasts were predicted to start accumulating in the wound area slightly after neutrophils — between 12 hours and 1 day after injury. Neutrophil and macrophage counts declined from day 3 after injury in the simulation, whereas fibroblasts remained at high levels until day 5. All ECM products, especially collagen type I, were predicted to be secreted in great quantities by activated fibroblasts after surgical trauma (Fig 3B). Collagen was predicted to be initially elevated at approximately days 1 to 3 after injury. In fact, collagen type I accumulation has been well described in the literature on scarred vocal folds following vocal fold “stripping” in different animal models.2
The trajectory of tissue damage was predicted to be high (Fig 3B). Damage slowly subsided until day 3, temporally corresponding to the decline of inflammatory cell counts (neutrophils and macrophages). A small rebound “spike” was seen between days 0.5 and 1 after injury. Of note, the timing of the “spike” in the damage curve was not programmed into the ABM code, but rather, emerged as a property of the system. This spike might be attributed to a positive feedback loop involved in the function of tissue damage. Inflammatory cells produce positive feedback to induce further inflammation because of their collateral ECM damage (Fig 1), thereby amplifying the proinflammatory response and delaying the tissue healing process.
In this study, the ABM was calibrated with experimental data across 2 biological organization levels, ie, molecules (inflammatory mediators) and tissues (ECM substances), with the long-term goal of reconstructing the link between molecular factors and their tissue-level manifestation. There was an overall good concordance between the predicted dynamics of these variables and the empirical data on mRNA expression of inflammatory and ECM markers in rats. Overall, compared to the acute phonotrauma ABMs,8,9 the current animal surgical ABM showed better predictive value, possibly because of the calibration of more time points (4 calibration data points in the current surgical ABM versus 2 calibration data points in the phonotrauma ABMs8,9). Another possible reason is that this animal data set is relatively cleaner than the human phonotrauma data, mainly because the experimental injury is comparatively consistent and controllable in animals. Also, the raw data used in the animal ABM were derived directly from tissues rather than estimated from secretions, as occurred for the human ABM. As a consequence, the healing outcomes became more predictable in the case of animals than in their human counterparts. If the animal data set ultimately can be extrapolated to the human case, this ABM could serve to predict a likely course of healing in patients after vocal fold surgery, if the degree of initial tissue injury can be established, as discussed below.
To date, the research of vocal fold pathology has mainly addressed damage at the tissue level. Vocal fold scarring due to surgical trauma has been characterized in terms of altered ECM structure and distribution,2 rather than by a profile of genetic or protein markers. Clinically, instrumental measures, such as phonation threshold pressure and an array of other aerodynamic and acoustic measures, are used to capture physiological signs of vocal fold tissue changes.87–91 These measures may capture accumulated damage but fail to effectively capture the early signs of tissue stress or damage, which may be needed in order to optimize treatment before overt histopathologic changes are detected in vocal fold tissues.
The trajectory of tissue damage generated from the surgical ABM was not prospectively calibrated with any empirical data in this study. Surprisingly, the predicted dynamics of this key variable (a surrogate for tissue health status) was found to be in good correspondence with recent empirical reports by other authors. A rat vocal fold study showed that epithelialization and fibroblast proliferation were observed on day 1 after mucosal stripping surgery, and fibroblasts were at their peak number on day 3 after surgery.29 These observations suggest that the degree of ultimate tissue damage might be controlled by the initiation of a fibroblast-dominant repair program starting on day 1 after injury. Also, the damaged tissues were expected to be replaced by neomatrix from day 3 after injury, at the time when fibroblasts were the most abundant in the wound site. In the current work, the predicted trajectory of tissue damage corresponded fairly well with the aforementioned empirical observations; ie, simulated damage decreased from day 1 after injury onward, nearly resolving by day 3 (Fig 3B).
The molecular correlate of the simulated vocal fold tissue damage has yet to be empirically verified. Indeed, this variable may represent a group of alarm or danger signals that are produced in response to inflammation and that further propagate inflammation. These signals are endogenous inducers of inflammation, which can be constitutive or inducible, intracellular or secreted, or even a part of the ECM.79,92,93 Some alarm or danger signals have been identified empirically, including heat-shock proteins, chromatin-associated protein high-mobility group box 1, adenosine-5′-triphosphate, uric acid, free DNA, IL-1α, IL-18, and degraded ECM components.78,83,94–102
Ideally, we would like to identify a marker that is sensitive to various stress levels (from low phonatory stress to high surgical stress) and indicates the earliest sign of tissue damage, well before conventional indicators or overt histopathologic changes of vocal fold damage can be seen. We focused on 3 ECM degradation products from collagen type I, elastin, and hyaluronan (HA).78,79 Collagen and elastin are both structural proteins in the vocal folds and are sparsely found in the vocal fold superficial lamina propria. These structural proteins may be more resistant to destruction than other ECM proteins in the vocal folds, at least within the range of physiological stresses typical of phonatory stresses examined in the present series. On the other hand, HA is abundant in the superficial lamina propria and has high turnover rates.103 Thus, we speculate that HA would not be as “inert” as its ECM counterparts and that it may be more promptly degraded even under the low levels of mechanical stress that occur during phonation. Hyaluronan fragments are generated during the moment of tissue damage and the subsequent inflammatory response and then are removed by activated macrophages in the wound site during the transition from inflammation to healing. In fact, HA fragments, because of their hydrophilic properties, have been reported to accumulate in inflamed tissues and contribute to edema formation, one of the cardinal signs of inflammation and wound repair.82,83 Failure to remove HA fragments was suggested to lead to persistent inflammation and possibly scar formation.104 Thus, the level of HA fragments seems to be sensitively associated with the state of the tissues. We hypothesize that the level of HA fragments could be an indicator of very early vocal fold stress or damage and might be in good concordance with the injury or trauma status of mucosal tissue. Currently, when a patient comes to the clinic with an existing vocal fold lesion, it is impossible for clinicians to accurately estimate the degree of tissue damage by examining the appearance of the vocal folds. Our future work will advocate defining the molecular correlate for the ABM-simulated tissue damage. Once the biomarker of tissue stress or damage is identified, the data will be used to quantitatively calibrate and validate the parameter of simulated tissue damage in ABMs and improve the model’s accuracy in predicting tissue status. Further, we assume that the ideal prescription type and dose for treatments will vary with tissue status. Once an in vivo measuring tool of tissue damage is available, clinicians will be able to increase their diagnostic accuracy and prescribe a tailored intervention for patients based on current tissue status.
Last, the model is intended to serve as a prototype for further developments, including in silico clinical trials and personalized diagnostics. In the latter case, when human data on inflammation and healing following vocal fold surgery become available, the model will be recalibrated and personalized to individuals with the ultimate goal of assisting clinicians to prescribe patient-specific treatment regimens that will optimize tissue healing after vocal fold scarring.
This study represents an attempt to advance systems-driven medicine in voice care. Computational modeling and the simulation of disease processes has become increasingly important for knowledge integration and hypothesis generation and testing. Future work will involve modeling various treatment effects on vocal fold injury, as well as attempts to define the molecular correlate for our ABM’s damage variable. It is hoped that these advances will lead to better diagnosis and treatment of vocal fold injuries.
The authors especially thank Dr Nathan Welham and Xinhong Lim for their generous help in sharing their published animal data for this project. The authors also thank Dr Kevin Kim for statistical advice. Last, the authors thank Dr John Durrant and Dr Susan Shaiman for providing valuable advice to improve this study, conducted as part of the doctoral dissertation project of Dr Li.
This work was funded by NIH grants R01-DC-008290and P50-GM-53789.