Next Article in Journal
Assessment of Deep Learning Methodology for Self-Organizing 5G Networks
Next Article in Special Issue
Toward Development of a Vocal Fold Contact Pressure Probe: Sensor Characterization and Validation Using Synthetic Vocal Fold Models
Previous Article in Journal
A Four-Stage Method for Active Control with Online Feedback Path Modelling Using Control Signal
Previous Article in Special Issue
Geometry of the Vocal Tract and Properties of Phonation near Threshold: Calculations and Measurements
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Towards a Physiological Scale of Vocal Fold Agent-Based Models of Surgical Injury and Repair: Sensitivity Analysis, Calibration and Verification

1
Department of Biological and Biomedical Engineering, McGill University, Montreal, QC H3A 0G4, Canada
2
School of Communication Sciences and Disorders, McGill University, Montreal, QC H3A 1G1, Canada
3
Department of Electrical and Computer Engineering, University of Maryland, College Park, MD 20742, USA
4
Virginia Tech Carilion Research Institute, Roanoke, VA 24016, USA
5
Department of Biology, McGill University, Montreal, QC H3A 1G1, Canada
6
Department of Mechanical Engineering, McGill University, Montreal, QC H3A 0G4, Canada
7
Department of Otolaryngology–Head and Neck Surgery, McGill University, Montreal, QC H3A 1G1, Canada
*
Author to whom correspondence should be addressed.
This author currently has no affiliation.
Appl. Sci. 2019, 9(15), 2974; https://doi.org/10.3390/app9152974
Submission received: 12 May 2019 / Revised: 3 July 2019 / Accepted: 15 July 2019 / Published: 25 July 2019
(This article belongs to the Special Issue Computational Methods and Engineering Solutions to Voice)

Abstract

:
Agent based models (ABM) were developed to numerically simulate the biological response to surgical vocal fold injury and repair at the physiological level. This study aimed to improve the representation of existing ABM through a combination of empirical and computational experiments. Empirical data of vocal fold cell populations including neutrophils, macrophages and fibroblasts were obtained using flow cytometry up to four weeks following surgical injury. Random Forests were used as a sensitivity analysis method to identify model parameters that were most influential to ABM outputs. Statistical Parameter Optimization Tool for Python was used to calibrate those parameter values to match the ABM-simulation data with the corresponding empirical data from Day 1 to Day 5 following surgery. Model performance was evaluated by verifying if the empirical data fell within the 95% confidence intervals of ABM outputs of cell quantities at Day 7, Week 2 and Week 4. For Day 7, all empirical data were within the ABM output ranges. The trends of ABM-simulated cell populations were also qualitatively comparable to those of the empirical data beyond Day 7. Exact values, however, fell outside of the 95% statistical confidence intervals. Parameters related to fibroblast proliferation were indicative to the ABM-simulation of fibroblast dynamics in final stages of wound healing.

1. Introduction

Voice disorders involve impairments with varying pitch, quality and loudness in daily conversation as well as in occupational voice use such as during voice performance and teaching. Among all voice disorders, vocal fold scarring is considered as one of the most perplexing clinical problems [1,2,3]. Current treatment options to fully repair the fibrotic tissue are deemed to be limited [2,3,4]. Vocal fold scars can develop within the lamina propria after surgical removal of benign or malignant vocal fold lesions. The scarred tissue alters the microarchitecture and vibratory functions of the vocal folds, resulting in a debilitating condition of the human voice, namely, dysphonia [2,3,4,5,6,7,8]. Iatrogenic vocal fold scarring triggers a cascade of cellular and molecular events associated with inflammation and wound repair, whereas the degree of response varies across patients depending on myriad personal and lifestyle factors. Computer simulation has become an appealing approach to integrate and analyze massive patient data towards more personalized and precise medical treatments.
Agent-based computational models (ABM) are often used to simulate complex system dynamics. The basic framework of ABM consists of agents and agent-rules. Agents are decision-making units that interact with other agents and the environment. Agent-rules are formulated to control the action and decision of agents in the virtual world [9]. ABM have been widely applied to ecology, economics, geography, business, political sciences and social sciences [10,11]. For biomedical applications, ABM are used to understand complex acute and chronic diseases such as diabetes, traumatic brain injury, spinal cord injuries, acute liver failure and sepsis [12,13,14,15,16]. ABM can be further developed to numerically simulate individual treatment response of a disease as a function of patient profile [17]. In silico experiments allow the exploration of a much wider parameter space and a much longer time scale than what would be very costly or sometimes impossible with animal and human models. In the era of personalized/ precision medicine, computational models have become an indispensable tool to integrate different forms of patient data (e.g., demographics, clinical, lifestyle etc.) and generate specific disease phenotypes or cohorts for in silico clinical trials [18,19].
Our laboratory has developed vocal fold agent-based models (VF-ABM) to numerically simulate the inflammatory and healing process following a surgical injury [20,21,22,23,24,25,26,27]. VF-ABM were composed of: (1) cells (platelets, macrophages, neutrophils and fibroblasts), (2) extracellular matrix (ECM) substances (collagen type I, elastin and hyaluronic acid (HA)) and (3) chemical mediators (tumor necrosis factor-α (TNF-α), interleukin-1β (IL-1β), interleukin-10 (IL-10), interleukin-6 (IL-6) interleukin-8 (IL-8), basic fibroblast growth factor (bFGF), and transforming growth factor (TGF-β1)). The model used initial wound size as the inputs to numerically simulate the post-injury trajectories of inflammatory cytokines and ECM markers in surgical vocal fold injury [24]. The dimensions of the VF-ABM correspond to the physiological size of rat vocal folds, since most of the empirical data were available for this animal model [28,29,30,31,32,33,34,35,36,37]. Partial calibration and verification were performed with vocal fold biological data [36,37,38,39,40]. However as quantitative data on vocal fold cellularity were limited, the cellular data have not been fully calibrated and verified in the current VF-ABM.
Model calibration is an important step to ensure sufficient authenticity of model representation to the real world. Model calibration involves systematically altering the values of parameters in the model in an iterative fashion until the simulated outputs and the observable behavior of the system are as close as possible. Some common parameter estimation methods for ABM calibration have been reported to reduce the uncertainty of the model (Table 1).
The pattern-oriented approach is used to estimate parameters by comparing observed trends and patterns in the model with empirical data [42,49,50]. For complex models having a large number of parameters, this approach has been limited to the identification of relevant patterns for discovering hidden information [50]. The particle swarm optimization (PSO) algorithm is a heuristic optimization method that is more suitable for computational models with a large number of parameters [45,51]. The PSO investigates the mutual collaboration between particles for sharing the internal information and reducing parameter space by discarding implausible inputs [45,51]. Since PSO has a low convergence rate in the iterative process, it is not considered as ideal for ABM, which would require thousands of iterations for each simulation [51]. Genetic algorithms (GA) are another heuristic optimization method that use a natural selection process to estimate as many solutions as possible simultaneously [46,52]. In GA, the model parameter values are progressively modified using genetic operators namely mutation and crossover, to optimize the model’s fitness in predicting the empirical data [46,52]. One limitation of GA is parameter overfitting and large variation in the results for each run [52]. Other existing tools for parameter estimation in ABM include parameter sweeping, Bayesian approaches, greedy algorithms and regressions, hybrid approaches and nonlinear multi-grid/finite difference methods. These methods are often computationally expensive and are not ideal for biological models at large scales [41,43,44,47,48].
When the model has a large number of parameters with unknown values, model calibration can also require expensive computing resources. The parameter search space scales up with the number of parameters in the model [53,54]. Lengthy simulation times, ambiguities in model design and a large number of unknown parameters can also limit the efficiency of model calibration [55]. Sensitivity analysis is thus required to reduce the number of parameters before model calibration. Sensitivity analysis is a statistical technique to explore the effects of systematically varying input parameters on the variability of the model outputs. Multiple sensitivity analysis methods have been reported for ABM [56,57]. The most common techniques are one parameter at a time (OPAT), Fourier amplitude sensitivity test (FAST) regression-based methods and variance-based methods [56,58,59]. These methods, however, are often computationally expensive or unreliable for large numbers of agents or they require excessive numerical simulation times [56,58,59,60,61]. Alternative schemes for sensitivity analysis are therefore worth exploring for running large-scale ABM with long-time scale simulations. For instance, Random forests, a popular machine learning technique, has been applied for the quantification of parameter importance and parameter selection [62,63,64]. Random forest was proposed to handle high dimensional data, thousands of input variables, complex parameter interaction, and missing data in complex models [65].
The primary goal of this study was to improve the accuracy of existing VF-ABM in simulating cellular and molecular activities associated with surgical vocal fold injury and repair. Experimental and computational studies were conducted. This paper was structured to present: (1) The experimental study of identifying and phenotyping immune and repair cells in surgically injured rat vocal folds up to four weeks after surgical injury; as well as (2) the computational study of optimizing the biological representation of VF-ABM through sensitivity analysis, model calibration and verification using flow cytometry data.

2. Materials and Methods

2.1. Vocal Fold Cell Phenotyping

The animal study was approved by the Institutional Animal Care and Use Committee of the University of Maryland-College Park (protocol number: R-12-85). The study design was a one-way between-subjects design with time [0 (uninjured control), one, two, three, or five days, and one, two, or four weeks post-surgery] as the independent variable. Dependent variables were percentages of neutrophils, macrophages, fibroblasts, and endothelial cells in the overall cell population based on multiparametric flow cytometry analyses.

2.1.1. Vocal Fold Surgical Model

A total of 160 male Sprague-Dawley rats (four to six months old, 450-500 g) were used in the study. Vocal fold injuries were performed in 140 rats using an established protocol. Twenty rats were used as uninjured controls [35]. Regarding the surgical protocol, animals were anesthetized, and their vocal folds were bilaterally injured using a 25-gauge needle to strip the vocal fold mucosa until the thyroarytenoid muscle was exposed [35]. At each end-point of the study, animals were euthanized via CO2 asphyxiation, and larynges were removed immediately following euthanasia [35]. For anesthesia induction, the rats were injected with isoflurane (2% to 3% delivered at 0.8 to 1.5 L/min) followed by maintenance with an intraperitoneal injection of ketamine hydrochloride (90 mg/ kg) and xylazine hydrochloride (9 mg/kg). Rats with injured vocal folds were euthanized at each of seven post-surgery time points: Days 1, 2, 3, 5 and Weeks 1, 2, 4. Due to the unexpected death of the animals, between 16 to 19 animals survived at each assigned time point for laryngeal harvest and the subsequent flow cytometry analysis.

2.1.2. Vocal Fold Cell Isolation

The mucosae of both sides of vocal folds were dissected and separated from the underlying thyroarytenoid muscles under a stereoscope. Vocal fold mucosal samples within an experimental group (e.g., 1 day after injury) and cytometric panel (see below) were pooled and subjected to single cell isolation. The strategy of sample pooling has been used for rat vocal fold cell and protein analysis considering the small size of rat vocal folds [66,67]. Samples were dispersed into single cell suspensions using digestion, centrifugation and filtration steps [66]. First, dissected mucosa were placed in conical panels covered in aluminum foil and then incubated in Ca/Mg free DPBS solution (Mediatech, USA, Cat#: 21-031-CV) with 0.05% collagenase (Gibco, USA, Cat#: 17018-029) /0.001% DNAse I (Sigma Aldrich, USA, Cat#: D4513) for 20 min at 37 °C to dissolve the extracellular matrix for cell release. The samples were thoroughly mixed and incubated again in 37 °C water bath. Once samples were visually confirmed to contain no large chunks of tissue, the solution was filtered with 40 nm cell strainers. Filtered samples were centrifuged at 290g for five minutes at 4 °C. After aspirating the supernatant, remaining cells were suspended in HBSS (Corning, USA, Cat#: 21-022-CV) from 1M HEPES solution with 2% fetal calf serum (HyClone, Canada, Cat#: SH30071.03) and 10nM HEPES (diluted from 1M HEPES solution, Corning, USA, Cat#: 25-060-CI). Resulted single cell suspension samples were transferred to flow cytometry polystyrene tubes and centrifuged at 201g for 10 min at 4 °C. After aspiration and suspension, the total cell number was determined with trypan blue staining and a hemacytometer. The number of isolated vocal fold cells ranged between 2.15 × 105 cells and 2.02 × 106 cells across time points.

2.1.3. Sample Preparation for Flow Cytometry

Experimental Samples

Two independent flow cytometry panels (Panels A and B) were designed to cross-validate the results for cell identification. Cell surface markers specific to the cells of interests were selected based on the literature of rat immunology. Multiple markers were used for cell identification since no single surface marker is specific to the cells of interest (Table 2 and Table 3). Panel A consisted of 11 parameters including FSC (cell size), SSC (cell granularity), one cell viability marker (AmCyan) and eight fluorescent cell surface markers (CD11b/c, CD29, CD44H, CD45, CD68, CD105, CD106 and His48). Panel B was composed of eight parameters including FSC, SSC, one cell viability marker (AmCyan) and five fluorescent cell surface markers (CD31, CD45, CD90, CD163 and His48) (see Table S1 and Table S2 for the corresponding biological functions of each marker in Supplementary Materials).
Pre-conjugated primary antibodies were used to facilitate specific affinity to surface antigens (see Table S3 and Table S4 for catalogue number and the fluorescence conjugate in Supplementary Materials). Isolated cells were first incubated in 1:100 diluted purified mouse anti-rat CD32 - FcγII blocker (Monoclonal D34-485, 0.5 mg/ml, BD) in staining buffer for 20 min at 4 °C. This blocking step was required to prevent non-specific antibody reactions to the Fc-receptor on cells such as monocytes and macrophages. Samples were then stained in 1:100 dilution of either eight (Panel A) or five (Panel B) preconjugated antibodies with staining buffer for 30 min at 4 °C in a dark room. At the end, we collected cell type information across a total of 16 datasets (8 time points × 2 flow panels).

Control Samples

Unstained samples were used as negative controls that contained 5% FBS and DPBS without Ca/Mg ions. Samples were washed with DPBS without Ca/Mg ions twice and then stained with 1:1000 diluted fixable viability dye in DPBS without Ca/Mg ions for 30 min at 4 °C in a dark room. A fixable viability dye was used to exclude dead cells from the analysable population. UltraComp eBeads, which are beads conjugated with individual fluorochromes, were used as single-color controls for compensation setup to correct the spectral overlap. For the fixable viability dye and UltraComp beads, procedures were performed according to manufacturers. Samples were washed with staining buffer twice and then transferred to clear polystyrene flow cytometry tubes until analysis. Samples were read using BD FACSAria II (BD Bio sciences, San Jose, CA, USA) in the Maryland Pathogen Research Institute. Compensation steps were performed to correct spectral overlap before every run of the flow analysis [80,81].

2.1.4. Flow Cytometry Data Processing and Analysis

Bivariate gating, an approach commonly used in flow cytometric analyses, was used primarily to analyze the multi-dimensional flow cytometry data [82,83,84,85,86,87,88]. FlowJo software (version 10.0.7, LLC, Ashland, OR, USA) was used in this study. Cell viability and compensation were applied during the setup of the experiment. For all datasets, gating analysis was performed by one individual (AG) who was blind to experimental conditions. The gating strategy for Panel A first used FSC-A vs SSC-A to represent the distribution of cells based on size, granularity and intracellular composition. These two scattering parameters led to the exclusion of debris, other non-cellular particles and lymphocytes [89,90]. To exclude doublets that were considered as single cells, three tests were implemented in the following order: (1) FSC-W vs FSC-H (W=width, H=height) to select low FSC-W cells, (2) FSC-W histogram to check the threshold and (3) FSC-H vs FSC-A (A=area) to select cells that were clustered diagonally [91,92].
A positive gating approach was used in which, instead of excluding negative cells (showing negative expression of a marker), priority was given to cells showing positive expression of a marker. The CD29 vs CD45 plot was used to distinguish CD29+CD45+ hematopoietic cells (neutrophils and macrophages) from CD29+CD45- non-hematopoietic cells (endothelial cells and fibroblasts) [68,71,72,73,74,75]. Here CD29 not only separated hematopoietic cells from non-hematopoietic cells, but it was also used as a conventional marker to gate out erythrocytes. For the gating of hematopoietic cells, CD45+His48+ cells were selected as neutrophils and were further verified using a His48+CD11b/c+ gates [77,78]. From the same gate of hematopoietic cells, macrophages were identified by using CD106-CD44H+ and His48-CD68+ gates [68,71,74,76]. In the case of non-hematopoietic cells (CD29+CD45-), CD106 was used to distinguish CD29+CD106+ endothelial cells from CD29+CD106- fibroblasts [68,73]. Endothelial cells and fibroblasts were then confirmed by using CD44H+CD106+ and CD29+CD105+ gates, respectively [68,71,73,75].
Regarding Panel B, FSC-A vs SSC-A and doublet exclusion tests were used in a similar way as for Panel A. CD90 was used with FSC-A to separate CD90-FSC-A+ hematopoietic cells from CD90+FSC-A+ non-hematopoietic populations [68,71,73,75]. CD31 was used to confirm the population of hematopoietic cells and to separate non-hematopoietic cells into endothelial cells and fibroblasts [68,71,73]. His48 was used to separate His48+CD31+ neutrophils from His48-CD31+ macrophages by using a His48 vs CD31 gate [77,78]. CD45 and CD163 were used to confirm populations of neutrophils and macrophages, respectively [68,71,72,73,74,75].
Multiple plots and parameters were used to decide the threshold limit in the generation of gates and to verify gate positions (see Figure S1 for illustration in Supplementary Materials). Contour and density plots (Figure S1A (iii) and (iv)), which show expression levels and relative density of data, are commonly used in placing the gates [86,93]. Backgating was used to verify the gate positions (Figure S1B). Backgating not only helped to analyse cells identified in a gate on dot plots with different parameters but also showed the final gated population within the population of its ancestors. The same gating strategy was applied across all eight time points within the same panel to compute the percentage of each cell type in the total cell population.

2.1.5. Flow Cytometry Statistical Analysis

Statistical analyses were performed using JMP software (SAS Institute, Cary, NC). Mixed effects models and correlation analyses were used to evaluate the degree of similarity in results from Panels A and B. Distinct mixed effects models were used to assess variation between panels for cell types (neutrophils, macrophages, endothelial cells and fibroblasts). Density data for a particular cell type (i.e., % of cells in total cell population) at each time point was obtained from one flow cytometry run. To analyze the effect of flow panel on estimates of cell density, we used the data for each of the eight time points as independent data points (i.e., time not included in the model, n = 8 for each panel). We additionally ran a full-factorial model with flow panel (A vs. B) and cell type as the independent variables, which allowed us to analyze the effect of the panel simultaneously across cell types. Because each flow cytometry run generated data for each cell type, the identity of each flow cytometry run was included as a random variable for the mixed effects model. Pearson’s correlations were used further to analyze the relationship between the results from Panel A and B, and to complement the results from mixed effects models.

2.2. Vocal Fold Agent-based Models

2.2.1. VF-ABM Implementation

The overview of the VF-ABM is illustrated in Algorithm 1. In our implementation, the ABM world represents a rat VF tissue and is spatially discretized into rectangular volumes called 3D patches. Each mobile agent represents a cell that moves from one patch to an adjacent patch and makes decisions to perform specific biological functions at discrete time steps. Agents make decisions based on the state of the patches. For example, chemokines and ECM proteins are associated with the states of the patches. Cells would perform their actions as a function of the concentrations of chemokines and ECM on their current and nearby patches. A detailed description of VF-ABM can be found in [25]. The overview of the VF-ABM is illustrated in Algorithm 1.
Algorithm 1: Overview of 3D rat vocal fold agent-based models (VF-ABM)
Procedure VFABM
 Initialization of patches
 Initialization of chemicals
 Initialization of cells
 Initialization of ECM
 Initialization of damage

 For each tick
   /* Model Computation */
   For each Patch
      Seed Cell Function
      ECM Function
      ECM Fragmentation
   For each Cell
      Cell Function
   For each Chemical Type
      Diffuse Chemical

   /* Model Update */
   For each Patch
      Update ECM
      Update Patch
      Update Chemicals
   For each Cell
      Update Cell
The initial model size and configuration of VF-ABM are summarized in Table 4. Agent-rules of VF-ABM were formulated based on reported mechanisms in the vocal fold wound healing literature [36,37,38,39,40]. Each time-step (or tick) corresponds to approximately 30 min of “real time”. VF-ABM were implemented on a computer node with two NVIDIA Tesla P100 12 GB Graphics Processing Units (GPUs) and 32 Intel E5-2683 v4 computer processing units (CPUs). GPUs were mainly used for running the chemical diffusion component of the model and implementing the visualizations of the simulations Open Graphics Library (OpenGL). The model was implemented in the object-oriented programming language C++ and Open Multi-Processing (OpenMP) for parallel computing [25,26,27]. Each simulation took about 22 min to complete for 28 days of “real time”. The source code of the VF-ABM can be found at https://github.com/VF-ABM/hpc-abm-vf-version_0_6.

2.2.2. Sensitivity Analysis

The goal of sensitivity analysis was to determine the influence of model parameters on the ABM cell count outputs in order to reduce the number of parameters for calibration by selecting the most important parameters. Random Forests were chosen as the sensitivity analysis method herein because it does not require many samples for the implementation and thus makes the computational cost feasible for the notable scale of VF-ABM.
Random Forests use an ensemble method with decision trees that are constructed using bootstrapping [94] (see Figure S2 for the workflow of Random Forests in Supplementary Materials). To split the node on any variable (or parameter) within a tree, the GINI impurity criterion (GINI Index) is used [95,96]. The GINI Index for two descendent nodes should be less than that of the parent node. A node is split if the change in GINI Index is significant [95], i.e.,
Δ i ( t ) = i ( t )   N L N i ( t L )   N R N i ( t R )
where, tL is the node on the left, and tR the node on the right, NL is the number of samples on left node and NR is the number of samples on right node.
The importance of any variable (or parameter) is determined by Mean decrease GINI, given by [95]
I ( P ) =   T N Δ i P   ( N ,   T )  
The mean decrease GINI (I) for a parameter (P) is evaluated by adding up the weighted GINI indices (i) for all nodes (N) where parameter P is used (averaged over all trees T in the forest) [95]. This quantity indicates how often a parameter P was selected for a split and provides a relative ranking of the parameters.
To find the suitable number of iterations, n for our sensitivity analysis, pilot tests were implemented with n = 3000, 4000, 4500, 5000, and 5500. It was observed that after 4500, the top three parameters for each cell type and time-point were identical and thus n = 5000 was chosen for this study. Sensitivity analysis using Random Forest was then implemented independently for the first four time-points (Day 1, Day 2, Day 3 and Day 5) as these time points would be used for the subsequent calibration.
The R package RandomForest was implemented in our computations (Algorithm 2) [97]. The user setting was set as (1) Number of Input parameters = 213, (2) Number of Output Parameters = Neutrophils, Macrophages, or Fibroblasts and (3) Number of Samples/ Number of Iterations of Model = 5000. The ranking of parameter importance was then obtained by sorting the parameters according to their Mean Decrease GINI.
Algorithm 2: Random Forests in R
Library – clusterGeneration
Library – mnormt
Require – randomForest
Library – caret

Number of Trees = 600

X=Samples
Y=Model Output for a time point T and cell type C

Df = data.frame (Y,X)
allX = paste (“X”, 1:ncol(X),sep=““)
names(df) = c(“Y”, allX)
fit= randomForest(factor(Y)~., data=df)
VI_F = importance(fit)
varImp(fit)
  varImpPlot(fit, type=2)

2.2.3. Model Calibration

The Statistical Parameter Optimization Tool for Python (SPOTPY) package was used for parameter calibration [98]. The SPOTPY package has a library of algorithms and objective functions for model calibration and verification. The algorithm has five major steps: (1) Sampling, (2) Simulation, (3) Evaluation, (4) Objective Function and (5) Parameter Estimation. Sampling was implemented for the most important parameters based on the sensitivity analysis using a uniform distribution and a predefined range. After sampling was done, the simulation was executed with the generated parameter set using the sampling function and the output was stored as simulated data. In our case, the empirical flow cytometry data were used as evaluation data. An objective function was then used to estimate the performance of the model for that given parameter set. The root mean square error (RMSE) was then used as the objective function, which determines the fitness of simulated data to evaluated data. RMSE is given by
R M S E = 1 / m i m ( e ( i ) s ( i ) ) 2  
Here, e(i) is evaluated data, s(i) is simulated data and m is the total number of parameters.
Lastly, the Robust Parameter Estimation (ROPE) algorithm was used for parameter optimization based on the concept of data depth [98]. Data depth refers to a statistical method that assigns a numerical value to a data point with respect to a set of points based on its centrality. It is a non-parametric multivariate statistical method that no distributional assumptions are needed [98,99]. ROPE was selected for VF-ABM calibration because it accelerates the calibration process by iteratively using information from previous simulations to estimate the outputs of subsequent ones [98,100]. The principle of ROPE is to identify a set of parameter vectors with the high model performance and subsequently generate a set of parameter vectors with high data depth with respect to the first set [98,100].
Instead of requiring thousands of random parameter sets, ROPE uses the results from the previous simulation as knowledge to define the parameter sets for the subsequent simulations. ROPE is recommended when the model requires parallel computing and the simulations are expensive, as the case of VF-ABM herein [98,100]. According to ROPE, the top 10% results are used to generate the next parameter sets. After estimating the best set of parameters, the algorithm repeats the sampling again. This time, the algorithm uses parameter sets learned from previous runs to determine the parameter set for the next run. The procedure repeats either until the defined number of iterations is reached or the RMSE reaches zero. In the present study, flow cytometry data from the first four time-points (Day 1, Day 2, Day 3 and Day 5) was used and the VF-ABM was run 800 times iteratively for each calibration step (Figure 1).
The top three parameters for each cell type (neutrophils, macrophages and fibroblasts) for each time point were used to proceed for calibration. Hence, a total of 36 top parameters were ranked for each cell type at each time point (3 parameters x 3 cell types x 4 time points); however, 12 of these parameters were found to be the same across certain conditions. As a result, a set of 24 unique parameters was calibrated iteratively with all other parameters being fixed as constants until the model eventually yielded a satisfactory match between simulation data and empirical data.
All four time points (Day 1, Day 2, Day 3 and Day 5) were first calibrated for each cell type individually and calibrated values of parameter sets for each cell type were determined. Subsequently, a new parameter range was reduced to 20% of the size of the original range, targeting the calibrated value, and was then re-entered manually in SPOTPY. The modified range for a particular parameter was estimated using
M i = { C     0.1 × ( M a M i ) ,   i f   [ C     0.1 × ( M a M i ) ]   >   M i   M i ,   i f   [ C     0.1 × ( M a M i ) ]   <   M i
M a = { C +   0.1 × ( M a M i ) ,   i f   [ C +   0.1 × ( M a M i ) ] <   M a M a ,   i f   [ C +   0.1 × ( M a M i ) ] >   M a
where, Mi’ and Ma’ are new minima and maxima, Mi and Ma are old minima and maxima and C is the calibrated value of the parameter. The parameters were then recalibrated again using the modified range with ROPE.

2.2.4. Model Verification

For model verification, the calibrated model was evaluated statistically for its accuracy in cell outputs at the Day 7, Week 2 and Week 4. Given ABM’s stochastic properties, the model was run 100 times up to Day 28 to generate a representative data set for statistical evaluation. To avoid making assumptions on the distributions of the output data, a non-parametric method was implemented to obtain predictive confidence intervals. A two-sided 95% confidence interval was computed for each cell type at each time-point by directly calculating the quantiles of the 100 simulation outputs using ordered statistics. The VF-ABM-simulated outputs would be considered accurate if the empirical results for a given cell type fall within the 95% confidence interval of the ABM-simulation outputs.

3. Results

3.1. Cell Compositions in Injured Rat Vocal Folds

Across both Panels A and B, flow cytometric analyses revealed variations in vocal fold cellular composition following injury (Figure 2A,B). Across all cell types, data from Panels A and B (i.e., two types of analyses) were not significantly different from each other (mixed effects model: F1,14 = 0.4, p = 0.5630) and significantly positively correlated with each other (Pearson’s correlation: p < 0.01 for pairwise correlations for each cell type). This result is consistent with results from the full-factorial mixed effects model (Panel and Cell Type as independent variables, cytometry run as a random variable), with Panel and the interaction between Panel and Cell Type not being significant (p > 0.05 for each). Importantly, the effect of cell type was significant (F3,42 = 59.3, p < 0.0001) and was driven by the fact that fibroblasts were more abundant than any other cell type analyzed (Tukey’s HSD, p < 0.05). No significant differences were found in cell density among the remaining cell types. These analyses suggest convergent estimates of cell types emerge from the distinct types of cytometry panels; consequently, we averaged the results from the two panels for the final characterization of changes to cell types across time (Figure 2C).
In uninjured controls, neutrophils, macrophages, endothelial cells and fibroblasts, respectively, represented 7.70%, 4.70%, 9.54% and 53.49% of the total cell population. In surgically injured vocal folds, fibroblasts also appeared as the prevalent cell type (>50%) among the four cell types at all time-points except Day 2 and Day 3. Neutrophils ranged from 2.31% to 14.37% on average throughout the examined time course. Macrophages appeared as the most dominant cell type (~33%) on Day 2. Shortly after, both neutrophils and macrophages were found in less than 5% of the total cell populations. The population of endothelial cells was below 5% in injured vocal folds on average except a relatively sharp increase during Week 2 (14.10%).

3.2. ABM Parameter Ranking

Random Forests were used to quantify the variance contribution of influential parameters to model outputs. Mean decrease GINI was used to rank the parameters from most influential to least influential. The ranking of parameter importance was then obtained by sorting the parameters according to their Mean Decrease GINI (Figure 3 for the top 25 parameters for Day 1; see Figures S3–S5 for the complete results in Supplementary Materials).

3.3. VF-ABM Calibration with Most Influential Parameters

VF-ABM were calibrated for Day 1, Day 2, Day 3 and Day 5 from the flow cytometry data. The top three parameters for each cell type for each time point were selected for model calibration (Table 5). As described in Section 2.2.3, 12 out of 36 parameters were overlapped across conditions, leading to a total of 24 unique parameters for the calibration. These top parameters were mostly related to cytokine synthesis, cell activation, ECM synthesis and cell recruitment (sprouting amount).

3.4. Verification of ABM-Predicted Cellular Dynamics

VF-ABM were verified for Day 7, Week 2, and Week 4 with the flow cytometry data. The VF-ABM reproduced the dynamics of all cell populations from the empirical data between Day 1 to Week 4 (Figure 4). For neutrophils, the peak was delayed by one day. Macrophage counts reached a maximum concentration on Day 2, in contrast to Day 1 in the data. For fibroblasts, the model reproduced the peak at the correct time point on Day 2 but did not resemble the oscillation pattern as observed in the empirical data.
In addition to the overall patterns, statistical tests of 95% confidence intervals were used to quantitatively evaluate the VF-ABM cell outputs. As the VF-ABM was run 100 times, a 95% confidence interval was obtained by computing the 2.5% and 97.5% quantiles from the ABM simulation output samples for each time point and cell type. The goal here was to evaluate for each time point (i.e., Day 7, Week 2 and Week 4), how many empirical data points (i.e., counts of neutrophils, macrophages and fibroblasts) fell within the two-sided 95% confidence interval of simulated outputs (Table 6). For Day 7, all empirical data were within the 95% confidence intervals (i.e., 100%, 3/3). For Week 2 and Week 4, all empirical data points fell outside the 95% intervals except for neutrophils at Week 4. However, when examining the data closely, the VF-ABM outputs were reasonably close to the empirical data. For instance, the empirical cell counts for neutrophils and macrophages were 180 and 186 at Week 2 respectively; whereas the VF-ABM outputs were from 188.5 to 222.5 and from 103.5 to 155.5 respectively.

4. Discussion

4.1. Time-Evolution of Cell Populations in Surgical Vocal Fold Injury and Repair

Multi-parametric flow cytometry was used to analyze dynamic changes to cell composition of vocal folds up to four weeks following injury. Results are generally consistent with the temporal dynamics of immune cells (neutrophils and macrophages) known in general wound healing literature. In dermal wound healing, neutrophils arrive at the wound site in the first few hours, peak in abundance one day after the injury, and reduce rapidly by the third day after injury [101,102,103,104,105,106,107]. Macrophages peak slightly later, by the second day after the injury, but also demonstrate rapid decreases on the third day after surgery [101,102,103,104]. Given the similarities in the temporal dynamics of neutrophil and macrophage changes following injury in vocal folds and other types of tissue, comparable cellular functions of immune cells underlying the inflammation and repair processes are speculated.
In addition to immune cells, the temporal pattern of endothelial cells following vocal fold damage also resembles the general wound healing literature [105,106]. In this study, endothelial cells contributed about 9% of the total cell population in uninjured controls that were likely originated from the endothelium layer of undisturbed blood capillaries in native vocal fold mucosae. A rise in the endothelial cell population was observed in Week 2 after surgery, representing about 14% of total cell populations. An active angiogenesis might have started between Week 1 and Week 2 post- surgery. To our best knowledge, this study is the first to evaluate the population of endothelial cells in vocal fold injury and healing. Further investigation on the precise role of endothelial cells, their interaction with other cell types as well as their implications in vocal fold scarring is warranted before we could reliably formulate agent-rules of endothelial cells in VF-ABM.
Fibroblasts normally arrive at the wound site around Day 3 and start the proliferative phase [37,103,105,106]. These cells usually reach their maximum concentrations between Day 5 and Day 7, and start to decrease gradually by the end of Week 1 [101,102,103,104]. Interestingly, fibroblasts were predominantly present in injured rat vocal folds spanning from the early acute inflammatory to later remodeling phases of wound healing as shown in this study. Previous studies reported the abundant presence of fibroblasts up to one week after surgical vocal fold injury [36,66,108]. Our data further indicated a notable accumulation of fibroblasts at the wound site up to four weeks after surgery, likely associated with their important roles in ECM remodeling.
Further, fibroblasts remain the dominant cell type at Day 1 (53.81%) and Day 2 (29.37%) despite the infiltration of neutrophils and macrophages into the acute wound environment. Previous work showed that both macrophages and fibroblasts were the major cell source of damage-associated molecular patterns (DAMP), namely high mobility group box-1, in modulating the inflammatory cytokine production in acute vocal fold injury [67,109]. Our data and others collectively suggest immunological and repair functions of fibroblasts that are unique in vocal fold microenvironment [110,111].

4.2. Random Forests and ROPE for VF-ABM Parameter Optimization

Common sensitivity analysis methods (e.g., One Parameter At a Time (OPAT), Fourier Amplitude Sensitivity Test (FAST) etc.) and calibration algorithms (e.g., Genetic Algorithm (GA), Particle Swamp Optimization (PSO), Bayesian approach etc.) are not optimal for large scale ABM as in our case herein. We thus used a new approach utilizing Random Forests and ROPE for sensitivity analysis and model calibration respectively. The execution of Random Forests does not require a large number of samples for running sensitivity analysis. This feature makes Random Forests less computationally expensive and more reliable for simulating long time scale models as compared to other sensitivity analysis methods [62,63,64,65]. For example, the estimated number of samples required for Fourier Amplitude Sensitivity Test (FAST) is 128 × parameter number ^ 2. That is, a total of 5,766,549 samples for the case of VF-ABM. An estimated time of running FAST is seven years with the computer node of two NVIDIA Tesla P100 12 GB GPUs and two Intel E5-2683 v4 CPUs [60,61]. In contrast, Random Forests uses a bootstrapping method for sampling and required only 5,000 samples for VF-ABM, making the sensitivity analysis much feasible [62,63,64,112].
Among the 24 parameter subjected to ABM calibration (Table 5), they were mostly related to cytokine synthesis, cell activation, ECM synthesis and cell recruitment (sprouting amount). In our VF-ABM, many biological activities of neutrophils, macrophages and fibroblasts are controlled by parameters of sprouting frequency and amount as well as cytokine synthesis. Specifically, parameters of sprouting frequency and amount were used to abstract the migration of neutrophils and macrophages to vocal fold vasculatures and transmigration from capillary to mucosal tissue as well as the recruitment and proliferation of fibroblasts in VF-ABM. Two sprouting-related parameters, namely parameter 155 (WhSproutAmount4) and parameter 156 (WhSproutAmount5) were used to control the amount of neutrophils. Three sprouting-related parameters (parameter 159 (WhSproutAmount8), 160 (WhSproutAmount9) and 161 (WhSproutAmount10) were used to control the amount of macrophages. Four sprouting-related parameters (parameter 150 (WhSproutFreq5), 162 (WhSproutAmount11), 163 (WhSproutAmount12) and 164 (WhSproutAmount13) were used to control the amount of fibroblasts. Further, the sprouting-related parameters for fibroblasts were linked to cytokines, growth factors and ECM contents in VF-ABM. Sprouting-related parameters are thus critical to initiate and sustain cell activities of neutrophils, macrophages and fibroblasts during tissue repair. The ABM outputs of cell numbers are mostly determined by the sprouting amount and sprouting frequency of these three cell types. As cell numbers were used as inputs for calibration, it is reasonable that the parameters of sprouting amount and frequency have the strongest influence on outputs of cell numbers in VF-ABM.

4.3. VF-ABM Performance

VF-ABM were calibrated for Day 1, Day 2, Day 3 and Day 5 and verified for Day 7, Week 2 and Week 4 with the flow cytometry data. The overall temporal trajectories of neutrophils, macrophages and fibroblasts were found to be in reasonable agreement between empirical and simulated data. Statistically, the empirical data of all three cell types fell within the VF-ABM output ranges for Day 7, but slightly outside the 95% confidence interval for remote time points (Week 2 and Week 4 except for neutrophils). Of note, the VF-ABM stimulated data showed a discrepancy of fibroblast count from the empirical trend at Week 4. The simulated fibroblast counts were lower than those of empirically observed by about 14 folds at Week 4. Based on the literature, proliferation rates of fibroblasts are 15.4 ± 1.1% after 4 days, 4.1 ± 0.6% after 1 week, and less than 0.5% after 2 weeks from rodent cardiac literature [113]. Since the proliferation rate of fibroblasts greatly depends on the microenvironment such as cytokine levels and ECM contents, the exact proliferation rate for vocal fold fibroblasts was not fully accurate based on the cardiac tissue data. Additional empirical data are thus needed to better estimate the proliferation rates of vocal fold fibroblasts in both homeostatic and injurious conditions. Once the data become available, parameters related to the sprouting amount and sprouting frequency of fibroblasts can be revised to improve the ABM-simulation of fibroblast dynamics in the final stages of wound healing such as wound remodeling.
We also noticed that VF-ABM neutrophil outputs did not have a perfect match with those from the empirical data in terms of timing and magnitude. (Figure 4A) The estimation of initial parameter range for corresponding sprouting frequencies (i.e., the rate of the cell infiltration when there is tissue damage) was difficult to define for calibration. The infiltration rate of neutrophils was only available as a function of blood flow in the human circulatory system [114]. However, specific information on rates of blood flow and injury-induced vasodilation within the vocal fold lamina propria are not available to date. Such data are necessary to better simulate the transmigration of neutrophils from capillary to mucosal tissue in VF-ABM.

5. Conclusions

The overarching goal of this research was to further develop VF-ABM into computer software that could pre-operatively inform clinicians about an individual’s risk of iatrogenic scarring from surgery. Iatrogenic vocal fold scarring is one of the most perplexing complications in phonosurgery. Surgical injury triggers a cascade of cellular and molecular events in inflammation and healing. Due to the individual-specific response to surgical injury, a computational tool could be useful to aid clinicians in tailoring the vocal treatment. For example, if the clinician can obtain patients’ biological profile before surgery, the computer model can help estimate individual post-operative response. In this study, empirical data on major vocal fold cell population were used to optimize existing VF-ABM through a scheme of sensitivity analysis, calibration and verification. General agreement was observed between the empirical and simulated data but further refinement is still warranted. The success of this work will contribute to the advancement of computational medicine in voice disorders.

Supplementary Materials

The following are available online at https://www.mdpi.com/2076-3417/9/15/2974/s1. Figure S1: Verification of Gating Strategy. Figure S2: Workflow of Random Forests. Figure S3: Top 25 parameters obtained by sensitivity analysis for Day 2. Figure S4: Top 25 parameters obtained by sensitivity analysis for Day 3. Figure S5: Top 25 parameters obtained by sensitivity analysis for Day 5. Table S1: Cell Surface Markers used in Panel A and Their Corresponding Biological Functions. Table S2: Cell Surface Markers used in Panel B and Their Corresponding Biological Functions. Table S3: Laser configuration of FACSAria II and Preconjugated primary antibody-fluorochrome list for Panel A. Table S4: Laser configuration of FACSAria II and Preconjugated primary antibody-fluorochrome list for Panel B. The flow cytometry datasets in this study are available on FlowRepository.org (ID: FR-FCM-ZYC6). The source code of the VF-ABM is available at https://github.com/VF-ABM/hpc-abm-vf-version_0_6. The sample data and source code of Random Forest and SPOTPY are available at https://github.com/VF-ABM/RF-SPOTPY.

Author Contributions

Conceptualization, A.G., N.L.J., N.S.; methodology, A.G., J.A.C.K., M.P., N.L.J.; software, N.S., J.J.J., S.Y., G.Y.; verification, A.G., S.Y., G.Y.; formal analysis, A.G., J.T.S.; investigation, A.G., J.A.C.K., M.P., N.L.J.; resources, L.M., J.J.J., N.L.J; data curation, A.G., G.Y., J.T.S.; writing—original draft preparation, A.G., N.L.J.; writing—review and editing, J.A.C.K., G.Y., J.T.S., N.L.J.; supervision, L.M., J.J.J., N.L.J.; project administration, N.L.J.; funding acquisition, N.L.J.

Funding

This research was funded by National Institute of Deafness and other Communication Disorder of the National Institutes of Health [R03DC012112], National Sciences and Engineering Research Council of Canada [RGPIN-2018-03843] and the Canadian Institutes of Health Research [388583].

Acknowledgments

We would like to acknowledge Susan Thibeault, Suzanne King, Martin Richer, David Mosser and his student, Rahul Suresh, for technical assistance and advice on flow cytometry, as well as Douglas Powell for providing trainings for the animal work. We would also like to acknowledge Compute Canada for their support of clusters and computational resources, which made the sensitivity analysis and calibration work possible.

Conflicts of Interest

The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Friedrich, G.; Dikkers, F.; Arens, C.; Remacle, M.; Hess, M.; Giovanni, A.; Duflo, S.; Hantzakos, A.; Bachy, V.; Gugatschka, M. Vocal fold scars: Current concepts and future directions. Consensus report of the phonosurgery committee of the European laryngological society. Eur. Arch. Oto Rhino Laryngol. 2013, 270, 2491–2507. [Google Scholar] [CrossRef] [PubMed]
  2. Hirano, S. Current treatment of vocal fold scarring. Curr. Opin. Otolaryngol. Head Neck Surg. 2005, 13, 143–147. [Google Scholar] [CrossRef] [PubMed]
  3. Rosen, C.A. Vocal Fold Scar. In The Otolaryngologica Clinics of North America. Voice Disorders and Phonosurgery; Rosen, C.A., Murry, T., Eds.; W.B. Saunders: Philadelphia, PA, USA, 2000. [Google Scholar]
  4. Benninger, M.S.; Alessi, D.; Archer, S.; Bastian, R.; Ford, C.; Koufman, J.; Sataloff, R.T.; Spiegel, J.R.; Woo, P. Vocal fold scarring: Current concepts and management. Otolaryngol. Neck Surg. 1996, 115, 474–482. [Google Scholar]
  5. Shin, Y.S.; Chang, J.W.; Yang, S.M.; Wu, H.W.; Cho, M.H.; Kim, C.H. Persistent dysphonia after laryngomicrosurgery for benign vocal fold disease. Clin. Exp. Otorhinolaryngol. 2013, 6, 166–170. [Google Scholar] [CrossRef] [PubMed]
  6. Woo, P.; Casper, J.; Colton, R.; Brewer, D. Diagnosis and treatment of persistent dysphonia after laryngeal surgery: A retrospective analysis of 62 patients. Laryngoscope 1994, 104, 1084–1091. [Google Scholar] [CrossRef] [PubMed]
  7. Perouse, A.R.; Coulombeau, B. Iatrogenic scarring of the vocal folds after phonosurgery for benign lesions. A descriptive study of 108 patients. Rev. Laryngol. Otol. Rhinol. (Bord) 2014, 135, 57–61. [Google Scholar]
  8. Hansen, J.K.; Thibeault, S.L. Current understanding and review of the literature: Vocal fold scarring. J. Voice 2006, 20, 110–120. [Google Scholar] [CrossRef]
  9. An, G.; Mi, Q.; Dutta-Moscato, J.; Vodovotz, Y. Agent-based models in translational systems biology. Wiley Interdiscip. Rev. Syst. Boil. Med. 2009, 1, 159–171. [Google Scholar] [CrossRef]
  10. Macal, C.M.; North, M.J. Agent-based modeling and simulation: ABMS examples. In Proceedings of the 2008 Winter Simulation Conference, Miami, FL, USA, 7–10 December 2008; pp. 101–112. [Google Scholar]
  11. Macal, C.M.; North, M.J. Tutorial on agent-based modeling and simulation. In Proceedings of the 2018 Winter Simulation Conference (WSC), Gothenburg, Sweden, 9–12 December 2005; p. 14. [Google Scholar]
  12. Mi, Q.; Li, N.Y.K.; Ziraldo, C.; Ghuma, A.; Mikheev, M.; Squires, R.; Okonkwo, D.; Verdolini Abbott, K.; Constantine, G.; An, G.; et al. Translational systems biology of inflammation: Applications to personalized medicine. Pers. Med. 2010, 7, 549–559. [Google Scholar] [CrossRef]
  13. Vodovotz, Y.; Constantine, G.; Faeder, J.; Mi, Q.; Rubin, J.; Bartels, J.; Sarkar, J.; Squires, R.; Okonkwo, D.; Gerlach, J.; et al. Translational systems approaches to the biology of inflammation and healing. Immunopharmacol. Immunotoxicol. 2010, 32, 181–195. [Google Scholar] [CrossRef]
  14. Clermont, G.; Bartels, J.; Kumar, R.; Constantine, G.; Vodovotz, Y.; Chow, C. In silico design of clinical trials: A method coming of age. Crit. Care Med. 2004, 32, 2061–2070. [Google Scholar] [CrossRef] [PubMed]
  15. Kumar, R.; Clermont, G.; Vodovotz, Y.; Chow, C.C. The dynamics of acute inflammation. J. Theor. Biol. 2004, 230, 145–155. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Vodovotz, Y.; Chow, C.C.; Bartels, J.; Lagoa, C.; Prince, J.M.; Levy, R.M.; Kumar, R.; Day, J.; Rubin, J.; Constantine, G.; et al. In silico models of acute inflammation in animals. Shock 2006, 26, 235–244. [Google Scholar] [CrossRef] [PubMed]
  17. De Castro, L.N.; Timmis, J. Artificial Immune Systems: A New Computational Intelligence Approach; Springer Science & Business Media: New York, NY, USA, 2002. [Google Scholar]
  18. Auffray, C.; Chen, Z.; Hood, L. Systems medicine: The future of medical genomics and healthcare. Genome Med. 2009, 1, 2. [Google Scholar] [CrossRef] [PubMed]
  19. Hood, L.; Flores, M. A personal view on systems medicine and the emergence of proactive P4 medicine: Predictive, preventive, personalized and participatory. New Biotechnol. 2012, 29, 613–624. [Google Scholar] [CrossRef] [PubMed]
  20. Li, N.Y.K. Biosimulation of Vocal Fold Inflammation and Healing. Ph.D. Thesis, University of Pittsburgh, Pittsburgh, PA, USA, 2009. [Google Scholar]
  21. Li, N.Y.K.; Abbott, K.V.; Rosen, C.; An, G.; Hebda, P.A.; Vodovotz, Y. Translational systems biology and voice pathophysiology. Laryngoscope 2010, 120, 511–515. [Google Scholar] [CrossRef] [PubMed]
  22. Li, N.Y.K.; Verdolini, K.; Clermont, G.; Mi, Q.; Rubinstein, E.N.; Hebda, P.A.; Vodovotz, Y. A patient-specific in silico model of inflammation and healing tested in acute vocal fold injury. PLoS ONE 2008, 3, e2789. [Google Scholar] [CrossRef]
  23. Li, N.Y.K.; Vodovotz, Y.; Kim, K.H.; Mi, Q.; Hebda, P.A.; Verdolini-Abbott, K. Biosimulation of acute phonotrauma: An extended model. Laryngoscope 2011, 121, 2418–2428. [Google Scholar] [CrossRef] [Green Version]
  24. Li, N.Y.K.; Vodovotz, Y.; Hebda, P.A.; Verdolini, K.A. Biosimulation of inflammation and healing in surgically injured vocal folds. Ann. Otol. Rhinol. Laryngol. 2010, 119, 412–423. [Google Scholar] [CrossRef]
  25. Seekhao, N.; Shung, C.; JaJa, J.; Mongeau, L.; Li-Jessen, N.Y.K. High-Performance Agent-based Modeling Applied to Vocal Fold Inflammation and Repair. Front. Physiol. 2018, 9, 304. [Google Scholar] [CrossRef]
  26. Seekhao, N.; JaJa, J.; Mongeau, L.; Li-Jessen, N.Y. In situ visualization for 3D agent-based vocal fold inflammation and repair simulation. Supercomput. Front. Innov. 2017, 4, 68. [Google Scholar] [PubMed]
  27. Seekhao, N.; Shung, C.; JaJa, J.; Mongeau, L.; Li-Jessen, N.Y. Real-time agent-based modeling simulation with in-situ visualization of complex biological systems: A case study on vocal fold inflammation and healing. In Proceedings of the 2016 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), Chicago, IL, USA, 23–27 May 2016; pp. 463–472. [Google Scholar]
  28. Cortés-Sol, A.; Lara-Garcia, M.; Alvarado, M.; Hudson, R.; Berbel, P.; Pacheco, P. Inner capillary diameter of hypothalamic paraventricular nucleus of female rat increases during lactation. BMC Neurosci. 2013, 14, 7. [Google Scholar] [CrossRef] [PubMed]
  29. Burns, A.R.; Smith, C.W.; Walker, D.C. Unique structural features that influence neutrophil emigration into the lung. Physiol. Rev. 2003, 83, 309–336. [Google Scholar] [CrossRef] [PubMed]
  30. Hohsfield, L.A.; Ammann, C.G.; Humpel, C. Inflammatory status of transmigrating primary rat monocytes in a novel perfusion model simulating blood flow. J. Neuroimmunol. 2013, 258, 17–26. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  31. Ling, C.; Yamashita, M.; Waselchuk, E.A.; Raasch, J.L.; Bless, D.M.; Welham, N.V. Alteration in cellular morphology, density and distribution in rat vocal fold mucosa following injury. Wound Repair Regen. 2010, 18, 89–97. [Google Scholar] [CrossRef] [Green Version]
  32. Tateya, I.; Tateya, T.; Bless, D.M. Homeostasis of hyaluronic acid in scarred rat vocal folds. In Proceedings of the International Conference on Voice Physiology and Biomechanics, Marseille, France, 18–20 August 2004. [Google Scholar]
  33. Tateya, I.; Tateya, T.; Sohn, J.-H.; Bless, D.M. Histological effect of basic fibroblast growth factor on chronic vocal fold scarring in a rat model. Clin. Exp. Otorhinolaryngol. 2016, 9, 56. [Google Scholar] [CrossRef] [PubMed]
  34. Tateya, T.; Tateya, I.; Munoz-del-Rio, A.; Bless, D.M. Postnatal development of rat vocal folds. Ann. Otol. Rhinol. Laryngol. 2006, 115, 215–224. [Google Scholar] [CrossRef]
  35. Welham, N.V.; Montequin, D.W.; Tateya, I.; Tateya, T.; Choi, S.H.; Bless, D.M. A rat excised larynx model of vocal fold scar. J. Speech Lang. Hear. Res. 2009, 52, 1008–1020. [Google Scholar] [CrossRef]
  36. Tateya, I.; Tateya, T.; Lim, X.; Sohn, J.H.; Bless, D.M. Cell production in injured vocal folds: A rat study. Ann. Otol. Rhinol. Laryngol. 2006, 115, 135–143. [Google Scholar] [CrossRef]
  37. Tateya, T.; Tateya, I.; Sohn, J.H.; Bless, D.M. Histological study of acute vocal fold injury in a rat model. Ann. Otol. Rhinol. Laryngol. 2006, 115, 285–292. [Google Scholar] [CrossRef]
  38. Welham, N.V.; Lim, X.; Tateya, I.; Bless, D.M. Inflammatory factor profiles one hour following vocal fold injury. Ann. Otol. Rhinol. Laryngol. 2008, 117, 145–152. [Google Scholar] [CrossRef] [PubMed]
  39. Lim, X.; Bless, D.M.; Munoz-Del-Rio, A.; Welham, N.V. Changes in cytokine signaling and extracellular matrix production induced by inflammatory factors in cultured vocal fold fibroblasts. Ann. Otol. Rhinol. Laryngol. 2008, 117, 227–238. [Google Scholar] [CrossRef] [PubMed]
  40. Tateya, T.; Tateya, I.; Sohn, J.H.; Bless, D.M. Histologic characterization of rat vocal fold scarring. Ann. Otol. Rhinol. Laryngol. 2005, 114, 183–191. [Google Scholar] [CrossRef] [PubMed]
  41. Folcik, V.A.; An, G.C.; Orosz, C.G. The Basic Immune Simulator: An agent-based model to study the interactions between innate and adaptive immunity. Theor. Biol. Med. Model. 2007, 4, 39. [Google Scholar] [CrossRef] [PubMed]
  42. Grimm, V.; Revilla, E.; Berger, U.; Jeltsch, F.; Mooij, W.M.; Railsback, S.F.; Thulke, H.-H.; Weiner, J.; Wiegand, T.; DeAngelis, D.L. Pattern-oriented modeling of agent-based complex systems: Lessons from ecology. Science 2005, 310, 987–991. [Google Scholar] [CrossRef] [PubMed]
  43. Gallaher, J.; Hawkins-Daarud, A.; Massey, S.C.; Swanson, K.; Anderson, A. Hybrid approach for parameter estimation in agent-based models. bioRxiv 2017, 175661. [Google Scholar] [Green Version]
  44. Hussain, F.; Langmead, C.J.; Mi, Q.; Dutta-Moscato, J.; Vodovotz, Y.; Jha, S.K. Automated parameter estimation for biological models using Bayesian statistical model checking. BMC Bioinform. 2015, 16, S8. [Google Scholar] [CrossRef]
  45. Li, T.; Cheng, Z.; Zhang, L. Developing a Novel Parameter Estimation Method for Agent-Based Model in Immune System Simulation under the Framework of History Matching: A Case Study on Influenza A Virus Infection. Int. J. Mol. Sci. 2017, 18, 2592. [Google Scholar] [CrossRef]
  46. Moedomo, R.L.; Pancoro, A.; Ibrahim, J.; Ahmad, A.S.; Mardiyanto, M.S.; Belatiff, M.B.; Tasman, H. Simulation of influenza pandemic based on genetic algorithm and agent-based modeling: A multi-objective optimization problem solving. J. Matemat. Dan Sains 2010, 15, 47–59. [Google Scholar]
  47. Tong, X.; Chen, J.; Miao, H.; Li, T.; Zhang, L. Development of an agent-based model (ABM) to simulate the immune system and integration of a regression method to estimate the key ABM parameters by fitting the experimental data. PloS ONE 2015, 10, e0141295. [Google Scholar] [CrossRef]
  48. Wise, S.M.; Lowengrub, J.S.; Frieboes, H.B.; Cristini, V. Three-dimensional multispecies nonlinear tumor growth—I: Model and numerical method. J. Theor. Biol. 2008, 253, 524–543. [Google Scholar] [CrossRef] [PubMed]
  49. Wiegand, T.; Jeltsch, F.; Hanski, I.; Grimm, V. Using pattern-oriented modeling for revealing hidden information: A key for reconciling ecological theory and application. Oikos 2003, 100, 209–222. [Google Scholar] [CrossRef]
  50. Grimm, V.; Railsback, S.F. Pattern-oriented modelling: A ‘multi-scope’for predictive systems ecology. Phil. Trans. R. Soc. B 2012, 367, 298–310. [Google Scholar] [CrossRef] [PubMed]
  51. Li, M.; Du, W.; Nian, F. An adaptive particle swarm optimization algorithm based on directed weighted complex network. Math. Probl. Eng. 2014, 2014, 1–7. [Google Scholar] [CrossRef]
  52. Katare, S.; Bhan, A.; Caruthers, J.M.; Delgass, W.N.; Venkatasubramanian, V.; Katare, S.; Bhan, A.; Caruthers, J.M.; Delgass, W.N.; Venkatasubramanian, V. A hybrid genetic algorithm for efficient parameter estimation of large kinetic models. Comput. Chem. Eng. 2004, 28, 2569–2581. [Google Scholar] [CrossRef]
  53. Crooks, A.; Castle, C.; Batty, M. Key challenges in agent-based modelling for geo-spatial simulation. Comput. Environ. Urban Syst. 2008, 32, 417–430. [Google Scholar] [CrossRef]
  54. Windrum, P.; Fagiolo, G.; Moneta, A. Empirical validation of agent-based models: Alternatives and prospects. J. Artif. Soc. Soc. Simul. 2007, 10, 8. [Google Scholar]
  55. Fehler, M.; Klügl, F.; Puppe, F. In Approaches for resolving the dilemma between model structure refinement and parameter calibration in agent-based simulations. In Proceedings of the Fifth International Joint Conference on Autonomous Agents and Multiagent Systems, Hakodate, Japan, 8–12 May 2006; pp. 120–122. [Google Scholar]
  56. Iooss, B.; Lemaître, P. A review on global sensitivity analysis methods. In Uncertainty Management in Simulation-Optimization of Complex Systems; Springer: Boston, MA, USA, 2015; pp. 101–122. [Google Scholar]
  57. Marino, S.; Hogue, I.B.; Ray, C.J.; Kirschner, D.E. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J. Theor. Biol. 2008, 254, 178–196. [Google Scholar] [CrossRef] [Green Version]
  58. Ten Broeke, G.; Van Voorn, G.; Ligtenberg, A. Which sensitivity analysis method should I use for my agent-based model? J. Artif. Soc. Soc. Simul. 2016, 19, 1–5. [Google Scholar] [CrossRef]
  59. Saltelli, A.; Tarantola, S.; Chan, K.-S. A quantitative model-independent method for global sensitivity analysis of model output. Technometrics 1999, 41, 39–56. [Google Scholar] [CrossRef]
  60. McRae, G.J.; Tilden, J.W.; Seinfeld, J.H. Global sensitivity analysis—A computational implementation of the Fourier amplitude sensitivity test (FAST). Comput. Chem. Eng. 1982, 6, 15–25. [Google Scholar] [CrossRef]
  61. Henkel, T.; Wilson, H.; Krug, W. In Global sensitivity analysis of nonlinear mathematical models—An implementation of two complementing variance-based algorithms. In Proceedings of the 2012 Winter Simulation Conference, Berlin, Germany, 9–12 December 2012; p. 154. [Google Scholar]
  62. Strobl, C.; Boulesteix, A.-L.; Kneib, T.; Augustin, T.; Zeileis, A. Conditional variable importance for random forests. BMC Bioinform. 2008, 9, 307. [Google Scholar] [CrossRef] [PubMed]
  63. Criminisi, A.; Shotton, J.; Konukoglu, E. Decision forests: A unified framework for classification, regression, density estimation, manifold learning and semi-supervised learning. Found. Trends® Comput. Gr. Vis. 2012, 7, 81–227. [Google Scholar] [CrossRef]
  64. Strobl, C.; Boulesteix, A.-L.; Zeileis, A.; Hothorn, T. Bias in random forest variable importance measures: Illustrations, sources and a solution. BMC Bioinform. 2007, 8, 25. [Google Scholar] [CrossRef] [PubMed]
  65. Breiman, L.; Cutler, A. Random Forests-Classification Description. Available online: https://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm (accessed on 20 May 2019).
  66. Ling, C.; Yamashita, M.; Zhang, J.; Bless, D.M.; Welham, N.V. Reactive response of fibrocytes to vocal fold mucosal injury in rat. Wound Repair Regen. 2010, 18, 514–523. [Google Scholar] [CrossRef] [PubMed]
  67. Li, N.Y.; Lee, B.J.; Thibeault, S.L. Temporal and spatial expression of high-mobility group box 1 in surgically injured rat vocal folds. Laryngoscope 2012, 122, 364–369. [Google Scholar] [CrossRef] [PubMed]
  68. Pilling, D.; Fan, T.; Huang, D.; Kaul, B.; Gomer, R.H. Identification of markers that distinguish monocyte-derived fibrocytes from monocytes, macrophages, and fibroblasts. PLoS ONE 2009, 4, e7475. [Google Scholar] [CrossRef] [PubMed]
  69. Smith, C.; Li, Z. Role of CD11a and CD11b in corneal wound healing and inflammatory process. Investig. Ophthalmol. Vis. Sci. 2004, 45, 125. [Google Scholar]
  70. Mizgerd, J.P.; Kubo, H.; Kutkoski, G.J.; Bhagwan, S.D.; Scharffetter-Kochanek, K.; Beaudet, A.L.; Doerschuk, C.M. Neutrophil emigration in the skin, lungs, and peritoneum: Different requirements for CD11/CD18 revealed by CD18-deficient mice. J. Exp. Med. 1997, 186, 1357–1364. [Google Scholar] [CrossRef]
  71. Biosciences, B. CD Marker Handbook. BD Biosciences, 2014. Available online: https://www.bdbiosciences.com/documents/cd_marker_handbook.pdf (accessed on 10 July 2019).
  72. Ahmed, N.; Vogel, B.; Rohde, E.; Strunk, D.; Grifka, J.; Schulz, M.B.; Grassel, S. CD45-positive cells of haematopoietic origin enhance chondrogenic marker gene expression in rat marrow stromal cells. Int. J. Mol. Med. 2006, 18, 233. [Google Scholar] [CrossRef]
  73. Kundrotas, G. Surface markers distinguishing mesenchymal stem cells from fibroblasts. Acta Medica Litu. 2012, 19. [Google Scholar] [CrossRef]
  74. Inoue, T.; Plieth, D.; Venkov, C.D.; Xu, C.; Neilson, E.G. Antibodies against macrophages that overlap in specificity with fibroblasts. Kidney Int. 2005, 67, 2488–2493. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  75. Van Landuyt, K.B.; Jones, E.A.; McGonagle, D.; Luyten, F.P.; Lories, R.J. Flow cytometric characterization of freshly isolated and culture expanded human synovial cell populations in patients with chronic arthritis. Arthritis Res. Ther. 2010, 12, R15. [Google Scholar] [CrossRef] [PubMed]
  76. Dijkstra, C.; Döpp, E.; Joling, P.; Kraal, G. The heterogeneity of mononuclear phagocytes in lymphoid organs: Distinct macrophage subpopulations in rat recognized by monoclonal antibodies ED1, ED2 and ED3. In Microenvironments in the Lymphoid System; Springer: London, UK, 1985; pp. 409–419. [Google Scholar]
  77. Fecho, K.; Lysle, D.T. Morphine-induced enhancement in the granulocyte response to thioglycollate administration in the rat. Inflammation 2002, 26, 259–271. [Google Scholar] [CrossRef] [PubMed]
  78. Iqbal, S.; Thomas, A.; Bunyan, K.; Tiidus, P.M. Progesterone and estrogen influence postexercise leukocyte infiltration in overiectomized female rats. Appl. Physiol. Nutr. Metab. 2008, 33, 1207–1212. [Google Scholar] [CrossRef] [PubMed]
  79. Strober, W. Trypan blue exclusion test of cell viability. Curr. Protoc. Immunol. 2001, 21, A.3B.1–A.3B.2. [Google Scholar]
  80. Roederer, M. Spectral compensation for flow cytometry: Visualization artifacts, limitations, and caveats. Cytom. Part A 2001, 45, 194–205. [Google Scholar] [CrossRef]
  81. Roederer, M. Compensation in flow cytometry. Curr. Protoc. Cytom. 2002, 22, 1.14.1–1.14.20. [Google Scholar]
  82. Perfetto, S.P.; Chattopadhyay, P.K.; Roederer, M. Seventeen-colour flow cytometry: Unravelling the immune system. Nat. Rev. Immunol. 2004, 4, 648–655. [Google Scholar] [CrossRef]
  83. Arnoulet, C.; Béné, M.C.; Durrieu, F.; Feuillard, J.; Fossat, C.; Husson, B.; Jouault, H.; Maynadié, M.; Lacombe, F. Four-and five-color flow cytometry analysis of leukocyte differentiation pathways in normal bone marrow: A reference document based on a systematic approach by the GTLLF and GEIL. Cytom. Part B Clin. Cytom. 2010, 78, 4–10. [Google Scholar] [CrossRef]
  84. Lugli, E.; Roederer, M.; Cossarizza, A. Data analysis in flow cytometry: The future just started. Cytom. Part A 2010, 77, 705–713. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  85. Autissier, P.; Soulas, C.; Burdo, T.H.; Williams, K.C. Evaluation of a 12-color flow cytometry panel to study lymphocyte, monocyte, and dendritic cell subsets in humans. Cytom. Part A 2010, 77, 410–419. [Google Scholar] [CrossRef] [PubMed]
  86. Pyne, S.; Hu, X.; Wang, K.; Rossin, E.; Lin, T.-I.; Maier, L.M.; Baecher-Allan, C.; McLachlan, G.J.; Tamayo, P.; Hafler, D.A. Automated high-dimensional flow cytometric data analysis. Proc. Natl. Acad. Sci. USA 2009, 106, 8519–8524. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  87. Lugli, E.; Pinti, M.; Nasi, M.; Troiano, L.; Ferraresi, R.; Mussi, C.; Salvioli, G.; Patsekin, V.; Robinson, J.P.; Durante, C. Subject classification obtained by cluster analysis and principal component analysis applied to flow cytometric data. Cytom. Part A 2007, 71, 334–344. [Google Scholar] [CrossRef] [PubMed]
  88. Steinbrich-Zöllner, M.; Grün, J.R.; Kaiser, T.; Biesen, R.; Raba, K.; Wu, P.; Thiel, A.; Rudwaleit, M.; Sieper, J.; Burmester, G.R. From transcriptome to cytome: Integrating cytometric profiling, multivariate cluster, and prediction analyses for a phenotypical classification of inflammatory diseases. Cytom. Part A 2008, 73, 333–340. [Google Scholar] [CrossRef] [PubMed]
  89. Wakabayashi, Y.; Watanabe, H.; Inoue, J.; Takeda, N.; Sakata, J.; Mishima, Y.; Hitomi, J.; Yamamoto, T.; Utsuyama, M.; Niwa, O. Bcl11b is required for differentiation and survival of αβ T lymphocytes. Nat. Immunol. 2003, 4, 533–539. [Google Scholar] [CrossRef] [PubMed]
  90. Shapiro, H.M. Practical Flow Cytometry; John Wiley & Sons: Hoboken, NJ, USA, 2005. [Google Scholar]
  91. Houtz, B.; Trotter, J.; Sasaki, D. BD FACService™ TECHNOTES. Custom. Focus. Solut. 2004, 9, 1–9. [Google Scholar]
  92. López-Sánchez, N.; Frade, J.M. Genetic evidence for p75NTR-dependent tetraploidy in cortical projection neurons from adult mice. J. Neurosci. 2013, 33, 7488–7500. [Google Scholar] [CrossRef] [PubMed]
  93. Herzenberg, L.A.; Tung, J.; Moore, W.A.; Herzenberg, L.A.; Parks, D.R. Interpreting flow cytometry data: A guide for the perplexed. Nat. Immunol. 2006, 7, 681–685. [Google Scholar] [CrossRef]
  94. Breiman, L. Random forests. Mach. Learn. 2001, 45, 5–32. [Google Scholar] [CrossRef]
  95. Menze, B.H.; Kelm, B.M.; Masuch, R.; Himmelreich, U.; Bachert, P.; Petrich, W.; Hamprecht, F.A. A comparison of random forest and its Gini importance with standard chemometric methods for the feature selection and classification of spectral data. BMC Bioinform. 2009, 10, 213. [Google Scholar] [CrossRef] [PubMed]
  96. Rodriguez-Galiano, V.; Chica-Olmo, M.; Abarca-Hernandez, F.; Atkinson, P.M.; Jeganathan, C. Random Forest classification of Mediterranean land cover using multi-seasonal imagery and multi-seasonal texture. Remote Sens. Environ. 2012, 121, 93–107. [Google Scholar] [CrossRef]
  97. Liaw, A.; Wiener, M. Classification and regression by randomForest. R News 2002, 2, 18–22. [Google Scholar]
  98. Houska, T.; Kraft, P.; Chamorro-Chavez, A.; Breuer, L. SPOTting model parameters using a ready-made python package. PLoS ONE 2015, 10, e0145180. [Google Scholar] [CrossRef] [PubMed]
  99. Krauße, T.; Cullmann, J. Towards a more representative parametrisation of hydrologic models via synthesizing the strengths of Particle Swarm Optimisation and Robust Parameter Estimation. Hydrol. Earth Syst. Sci. 2012, 16, 603–629. [Google Scholar] [CrossRef]
  100. Bárdossy, A.; Singh, S. Robust estimation of hydrological model parameters. Hydrol. Earth Syst. Sci. 2008, 12, 1273–1283. [Google Scholar] [CrossRef] [Green Version]
  101. Kirsner, R.S.; Eaglstein, W. The wound healing process. Dermatol. Clin. 1993, 11, 629–640. [Google Scholar] [CrossRef]
  102. Martin, P. Wound healing--aiming for perfect skin regeneration. Science 1997, 276, 75–81. [Google Scholar] [CrossRef] [PubMed]
  103. Witte, M.B.; Barbul, A. General principles of wound healing. Surg. Clin. N. Am. 1997, 77, 509–528. [Google Scholar] [CrossRef]
  104. Robson, M.C.; Steed, D.L.; Franz, M.G. Wound healing: Biologic features and approaches to maximize healing trajectories. Curr. Probl. Surg. 2001, 38, A1–A140. [Google Scholar] [CrossRef] [PubMed]
  105. Broughton, G., 2nd; Janis, J.E.; Attinger, C.E. The basic science of wound healing. Plast. Reconstr. Surg. 2006, 117, 12S–34S. [Google Scholar] [CrossRef] [PubMed]
  106. Enoch, S.; Leaper, D.J. Basic science of wound healing. Surgery (Oxford) 2008, 26, 31–37. [Google Scholar] [CrossRef]
  107. Pohlman, T.H.; Stanness, K.A.; Beatty, P.G.; Ochs, H.D.; Harlan, J.M. An endothelial cell surface factor (s) induced in vitro by lipopolysaccharide, interleukin 1, and tumor necrosis factor-alpha increases neutrophil adherence by a CDw18-dependent mechanism. J. Immunol. 1986, 136, 4548–4553. [Google Scholar] [PubMed]
  108. Branski, R.C.; Verdolini, K.; Rosen, C.A.; Hebda, P.A. Acute vocal fold wound healing in a rabbit model. Ann. Otol. Rhinol. Laryngol. 2005, 114, 19–24. [Google Scholar] [CrossRef] [PubMed]
  109. Li-Jessen, N.Y.; Powell, M.; Choi, A.J.; Lee, B.J.; Thibeault, S.L. Cellular source and proinflammatory roles of high-mobility group box 1 in surgically injured rat vocal folds. Laryngoscope 2016, 127, E193–E200. [Google Scholar] [CrossRef] [Green Version]
  110. King, S.N.; Chen, F.; Jetté, M.E.; Thibeault, S.L. Vocal fold fibroblasts immunoregulate activated macrophage phenotype. Cytokine 2013, 61, 228–236. [Google Scholar] [CrossRef]
  111. Hanson, S.E.; Kim, J.; Johnson, B.H.Q.; Bradley, B.; Breunig, M.J.; Hematti, P.; Thibeault, S.L. Characterization of mesenchymal stem cells from human vocal fold fibroblasts. Laryngoscope 2010, 120, 546–551. [Google Scholar] [CrossRef] [Green Version]
  112. Jiang, P.; Wu, H.; Wang, W.; Ma, W.; Sun, X.; Lu, Z. MiPred: Classification of real and pseudo microRNA precursors using random forest prediction model with combined features. Nucleic Acids Res. 2007, 35, W339–W344. [Google Scholar] [CrossRef]
  113. Virag, J.I.; Murry, C.E. Myofibroblast and endothelial cell proliferation during murine myocardial infarct repair. Am. J. Pathol. 2003, 163, 2433–2440. [Google Scholar] [CrossRef]
  114. Butterfield, T.A.; Best, T.M.; Merrick, M.A. The dual roles of neutrophils and macrophages in inflammation: A critical balance between tissue damage and repair. J. Athl. Train. 2006, 41, 457. [Google Scholar]
Figure 1. Workflow of VF-ABM Calibration. A total of 36 top parameters were ranked for each cell type at each time point. Twelve of them were overlapped across conditions. As a result, a total of 24 unique parameters were used for calibration.
Figure 1. Workflow of VF-ABM Calibration. A total of 36 top parameters were ranked for each cell type at each time point. Twelve of them were overlapped across conditions. As a result, a total of 24 unique parameters were used for calibration.
Applsci 09 02974 g001
Figure 2. Percentage of each cell type over time in Panel A and B with gating method. Behaviour of four cell types in two different panels: (A) Panel A and (B) Panel B. (C) Combination of two panels listed above in one curve.
Figure 2. Percentage of each cell type over time in Panel A and B with gating method. Behaviour of four cell types in two different panels: (A) Panel A and (B) Panel B. (C) Combination of two panels listed above in one curve.
Applsci 09 02974 g002
Figure 3. Top 25 most influential parameters for Day 1. (A) Neutrophils; (B) Macrophages; (C) Fibroblasts. The biological function of each parameter is coded by the shade of the bar.
Figure 3. Top 25 most influential parameters for Day 1. (A) Neutrophils; (B) Macrophages; (C) Fibroblasts. The biological function of each parameter is coded by the shade of the bar.
Applsci 09 02974 g003
Figure 4. Empirical and ABM-simulated trajectories of cell populations in injured vocal folds. (A) Number of Neutrophils. (B) Number of Macrophages. (C) Number of Fibroblasts.
Figure 4. Empirical and ABM-simulated trajectories of cell populations in injured vocal folds. (A) Number of Neutrophils. (B) Number of Macrophages. (C) Number of Fibroblasts.
Applsci 09 02974 g004
Table 1. Parameter estimation methods in agent-based modeling (ABM).
Table 1. Parameter estimation methods in agent-based modeling (ABM).
StudyCalibration MethodNumber of Unknown ParametersField of ABM Application
Folcik et al., 2007 [41]Parameter Sweeping87Basic Immune Simulator
Grimm et al., 2005 [42]Pattern-Oriented ApproachunreportedEcology
Gallaher et al., 2017 [43]Genetic algorithm with random weighted sampling16Glioblastoma multiforme model
Hussain et al., 2015 [44]Bayesian Approach24Dynamics of acute inflammation
Li et al., 2017 [45]Particle swarm optimization (PSO)50Immune System
Moedomo et al., 2010 [46]Genetic Algorithm 6Avian Influenza (H5N1) viruses mutation
Tong et al., 2015 [47]Greedy algorithm and Regression4Immune System
Wise et al., 2008 [48]Nonlinear multigrid/finite difference method20Three-dimensional multispecies nonlinear tumor growth
Table 2. Rat cell surface marker profile for flow cytometry in Panel A. “+” denotes positive expression and “-” denotes negative expression.
Table 2. Rat cell surface marker profile for flow cytometry in Panel A. “+” denotes positive expression and “-” denotes negative expression.
MarkerFluorochromeNeutrophil MacrophageEndothelial CellFibroblastReferences
CD11b/cFITC++--[68,69,70,71]
CD29PE-Cy7++++[68,71,72,73]
CD44HAPC-Cy7-+++[68,71,73]
CD45PerCP-Cy5.5++--[68,71,72,73,74,75]
CD68PE-Texas Red++-+[68,71,74,76]
CD105PE-+++[68,71,73,75]
CD106Brilliant Violet 421--+-[68,73]
His48APC+---[77,78]
Cell Viability DyeAmCyan++++[79]
Table 3. Rat cell surface marker profile for flow cytometry in Panel B. “+” shows positive expression and “-” shows negative expression.
Table 3. Rat cell surface marker profile for flow cytometry in Panel B. “+” shows positive expression and “-” shows negative expression.
MarkerFluorochromeNeutrophilMacrophageEndothelial CellFibroblastReferences
CD31APC+++-[68,71,73]
CD45 PerCP-Cy5.5++--[68,71,72,73,74,75]
CD90FITC--++[68,71,73,75]
CD163PE-Cy7-+--[68,71,75]
His48PE+---[77,78]
Cell Viability DyeAmCyan++++[79]
Table 4. Initial configurations of VF-ABM.
Table 4. Initial configurations of VF-ABM.
ItemUnitSizeReference
Vocal Fold Widthmm1[35]
patches142
Vocal Fold Heightmm1[35]
patches142
Vocal Fold Thicknessmm0.2[35]
patches28
Vocal Fold Epithelium Thicknessmm0.01[31]
patches1
Capillary Diameterµm7[28]
patches1
Capillary Gapµm12.89[28]
patches1
Patch sizeµm7 × 7 × 7
Total number of patchespatches564,592-
Total number of non-epi patchespatches544,428-
Total number of capillary patchespatches138,450-
Total number of tissue patchespatches405,978-
Simulated time-stepMinutes30-
Neutrophilsµm7[29]
Cells517Flow data
Macrophagesµm6[30]
Cells316Flow data
Fibroblastsµm6[29]
Cells3594Flow data
Table 5. Parameters used for model calibration.
Table 5. Parameters used for model calibration.
Time PointCellParameterBiological Significance
Day 1Neutrophils156--WhSproutAmount5Determine the effect of damage on the rate of neutrophil extravasation
164--WhSproutAmount13Determine the baseline rate of fibroblast sprouting
39--MacCytSynth14Determine the macrophage’s baseline synthesis rate of TNF-alpha
Macrophages148--WhSproutFreq3Frequency macrophages enter vocal fold capillaries when damage is present
117--FibCytSynth10Determine a fibroblast’s baseline synthesis rate of FGF
133--FibCytSynth26Determine a fibroblast’s baseline synthesis rate of IL-8
Fibroblasts200--FibProlif0Frequency of inactivated fibroblast proliferation in hours
150--WhSproutFreq5Frequency of fibroblasts sprouting in tissue in the presence of damage
166--FibActivat1Probability (in percent) of fibroblast activation in the absence of TGF-beta
Day 2Neutrophils156--WhSproutAmount5Determine the effect of damage on the rate of neutrophil extravasation
181--FibECMsynth0Determine the fibroblast’s baseline synthesis rate of collagen
103--MacCytSynth78Determine the macrophage’s baseline synthesis rate of IL-10
Macrophages159--WhSproutAmount8Determine the baseline rate of macrophage extravasation in the presence of damage
158--WhSproutAmount7Determine the effect of damage on the rate of macrophage extravasation
177--NeuActivat2Probability (in percent) of neutrophil activation in the presence of high TNF-alpha concentration
Fibroblasts14--NeuCytSynth10Effect of local TGF-beta concentration on neutrophil MMP-8 synthesis
72--MacCytSynth47Effect of local IL-1beta concentration on macrophage IL-6 synthesis
11--NeuCytSynth7Determine the neutrophil’s baseline synthesis rate of MMP-8
Day 3Neutrophils156--WhSproutAmount5Determine the effect of damage on the rate of neutrophil extravasation
200--FibProlif0Frequency of inactivated fibroblast proliferation in hours
148--WhSproutFreq3Frequency macrophages enter vocal fold capillaries when damage is present
Macrophages149--WhSproutFreq4Frequency macrophages extravasate in the presence of damage
160--WhSproutAmount9 Determine the effect of damage on the rate of macrophage extravasation
103--MacCytSynth78Determine the macrophage’s baseline synthesis rate of IL-10
Fibroblasts110--FibCytSynth3Determine the fibroblast’s baseline synthesis rate of TNF-alpha
51--MacCytSynth26Effect of local TNF-alpha concentration on macrophage IL1-beta synthesis
117--FibCytSynth10Determine a fibroblast’s baseline synthesis rate of FGF
Day 5Neutrophils156--WhSproutAmount5Determine the effect of damage on the rate of neutrophil extravasation
155--WhSproutAmount4Determine the baseline rate of neutrophil extravasation in the presence of damage
153--WhSproutAmount2Determine the rate of neutrophil extravasation in the absence of damage
Macrophages154--WhSproutAmount3Determine the rate of macrophage extravasation in the absence of damage
159--WhSproutAmount8Determine the baseline rate of macrophage extravasation in the presence of damage
149--WhSproutFreq4Frequency macrophages extravasate in the presence of damage
Fibroblasts200--FibProlif0Frequency of inactivated fibroblast proliferation in hours
51--MacCytSynth26Effect of local TNF-alpha concentration on macrophage IL1-beta synthesis
110--FibCytSynth3Determine the fibroblast’s baseline synthesis rate of TNF-alpha
Table 6. Cell quantities from empirial and simulation experiments from Day 7 to Week 4. * denotes the empirical data within the 95% confidence interval (CI) of VF-ABM outputs.
Table 6. Cell quantities from empirial and simulation experiments from Day 7 to Week 4. * denotes the empirical data within the 95% confidence interval (CI) of VF-ABM outputs.
Time-PointCell TypesEmpirical Data of Cell Quantities95% Confidence Interval (and Mean)
of ABM-Simulated Cell Quantities
Day 7Neutrophils214198.5–228 * (212.29)
Macrophages359310–406 * (359.15)
Fibroblasts57514920.5–7541.5 * (6083.45)
Empirical Data within 95% CI100%
Week 2Neutrophils180188.5–222.5 (204.87)
Macrophages186103.5–155.5 (126.8)
Fibroblasts38541759–2640 (2142)
Empirical Data within 95% CI0%
Week 4Neutrophils214192.5–227 * (210.16)
Macrophages25210–27 (17.07)
Fibroblasts4347254–378 (304.45)
Empirical Data within 95% CI33%

Share and Cite

MDPI and ACS Style

Garg, A.; Yuen, S.; Seekhao, N.; Yu, G.; Karwowski, J.A.C.; Powell, M.; Sakata, J.T.; Mongeau, L.; JaJa, J.; Li-Jessen, N.Y.K. Towards a Physiological Scale of Vocal Fold Agent-Based Models of Surgical Injury and Repair: Sensitivity Analysis, Calibration and Verification. Appl. Sci. 2019, 9, 2974. https://doi.org/10.3390/app9152974

AMA Style

Garg A, Yuen S, Seekhao N, Yu G, Karwowski JAC, Powell M, Sakata JT, Mongeau L, JaJa J, Li-Jessen NYK. Towards a Physiological Scale of Vocal Fold Agent-Based Models of Surgical Injury and Repair: Sensitivity Analysis, Calibration and Verification. Applied Sciences. 2019; 9(15):2974. https://doi.org/10.3390/app9152974

Chicago/Turabian Style

Garg, Aman, Samson Yuen, Nuttiiya Seekhao, Grace Yu, Jeannie A. C. Karwowski, Michael Powell, Jon T. Sakata, Luc Mongeau, Joseph JaJa, and Nicole Y. K. Li-Jessen. 2019. "Towards a Physiological Scale of Vocal Fold Agent-Based Models of Surgical Injury and Repair: Sensitivity Analysis, Calibration and Verification" Applied Sciences 9, no. 15: 2974. https://doi.org/10.3390/app9152974

APA Style

Garg, A., Yuen, S., Seekhao, N., Yu, G., Karwowski, J. A. C., Powell, M., Sakata, J. T., Mongeau, L., JaJa, J., & Li-Jessen, N. Y. K. (2019). Towards a Physiological Scale of Vocal Fold Agent-Based Models of Surgical Injury and Repair: Sensitivity Analysis, Calibration and Verification. Applied Sciences, 9(15), 2974. https://doi.org/10.3390/app9152974

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop