Multiscale Coupling of an Agent-Based Model of Tissue Fibrosis and a Logic-Based Model of Intracellular Signaling

Wound healing and fibrosis following myocardial infarction (MI) is a dynamic process involving many cell types, extracellular matrix (ECM), and inflammatory cues. As both incidence and survival rates for MI increase, management of post-MI recovery and associated complications are an increasingly important focus. Complexity of the wound healing process and the need for improved therapeutics necessitate a better understanding of the biochemical cues that drive fibrosis. To study the progression of cardiac fibrosis across spatial and temporal scales, we developed a novel hybrid multiscale model that couples a logic-based differential equation (LDE) model of the fibroblast intracellular signaling network with an agent-based model (ABM) of multi-cellular tissue remodeling. The ABM computes information about cytokine and growth factor levels in the environment including TGFβ, TNFα, IL-1β, and IL-6, which are passed as inputs to the LDE model. The LDE model then computes the network signaling state of individual cardiac fibroblasts within the ABM. Based on the current network state, fibroblasts make decisions regarding cytokine secretion and deposition and degradation of collagen. Simulated fibroblasts respond dynamically to rapidly changing extracellular environments and contribute to spatial heterogeneity in model predicted fibrosis, which is governed by many parameters including cell density, cell migration speeds, and cytokine levels. Verification tests confirmed that predictions of the coupled model and network model alone were consistent in response to constant cytokine inputs and furthermore, a subset of coupled model predictions were validated with in vitro experiments with human cardiac fibroblasts. This multiscale framework for cardiac fibrosis will allow for systematic screening of the effects of molecular perturbations in fibroblast signaling on tissue-scale extracellular matrix composition and organization.

Introduction

Approximately 605,000 Americans experience their first myocardial infarction (MI) each year, and another 200,000 experience a recurrent MI (Benjamin et al., 2018). Approximately 82% of males and 77% of females survive at least 1 year following their MI (Benjamin et al., 2018), making management of post-MI recovery an increasingly important topic.

Wound healing and scar remodeling following MI is a dynamic process involving many cell types, extracellular matrix, and inflammatory cues. Myocyte death due to prolonged ischemia initiates an inflammatory response led by cytokines such as IL-1β and TNFα (Frantz et al., 2009). Neutrophils and macrophages are recruited to the wound site within 24 h and begin to phagocytose debris and propagate the inflammatory response. Inflammatory cells peak within the first week of wound healing and then begin to subside as the proliferative phase begins (Czubryt, 2012). Inflammatory macrophages secrete TGFβ, which stimulates fibroblast recruitment and proliferation (Lambert et al., 2008). The release of TGFβ may also contribute to the conversion of macrophages to a more anti-inflammatory phenotype (Lambert et al., 2008). The proliferative phase may last for days to weeks and is marked by the proliferation of fibroblasts and transition to a myofibroblast phenotype, along with synthesis of many ECM components including collagen (Czubryt, 2012). ECM deposition produces a scar in the infarct region that contributes to its structural stability during wound healing. This proliferative phase is followed by weeks to months of scar remodeling and significant ECM turnover.

Post-MI cardiac wound healing is a complex and dynamic process with many overlapping phases. The cardiac fibroblast is the key effector cell throughout the phases of wound healing that creates and remodels scar tissue (Spinale et al., 2016; Mouton et al., 2019). However, fibroblasts are a highly dynamic and plastic cell type that can transition from a pro-inflammatory phenotype in the early phases of wound healing to an anti-inflammatory and pro-fibrotic phenotype later in the wound healing cascade (Chen and Frangogiannis, 2013; Mouton et al., 2019). Fibroblast response to single cytokine inputs are well-documented (Fredj et al., 2005; Fix et al., 2011; Turner, 2014), but fibroblast activation and cytokine secretion in response to multiple cytokines and other stimuli in vivo that shift over the time course of MI wound healing are not well-described (Ma et al., 2017). This lack of understanding of activation shifts over the time course of healing is at the core of the failure of many attempts to improve post-MI wound healing by modulating scar formation (Clarke et al., 2016). Inhibition of inflammation too early in the wound healing cascade can lead to thinning of the LV wall and scar rupture (Brown et al., 1983; Hammerman et al., 1983a,b). Aberrant fibrosis can lead to LV dilation and heart failure. This inherent complexity of the biological phenomenon necessitates the development of computational models to design and test therapeutic interventions that potentially have opposite effects at different phases throughout the wound healing cascade. Previous computational models have extensively characterized cardiac fibroblast signaling pathways and expression profiles to provide information about fibroblast activation and kinetics (Nim et al., 2015; Zeigler et al., 2016a,b), but fibroblast activation has generally been studied in response to single stimuli in vitro. Other researchers in the field have noted the need to understand fibroblast activation in response to mixed stimuli, and have called for the development of computational models that can integrate the effects of spatial and temporal shifts in fibroblast activation, with the cell-cell interactions and cell-matrix interactions that coordinate the short and long-term remodeling of scar tissue (Ma et al., 2017). A multiscale model that can translate cardiac fibroblast gene and protein expression to tissue level functional remodeling with spatial and temporal precision could provide an invaluable platform for identifying, testing, and validating new therapeutic interventions for inducing functional regeneration and mitigating fibrosis.

Our group has recently developed computational models to study distinct scales of cardiac wound healing, including a logic-based differential equation (LDE) model of intracellular signaling in individual cardiac fibroblasts and an agent-based model (ABM) of collagen remodeling by multiple cells in the infarct (Rouillard and Holmes, 2012; Zeigler et al., 2016a). Each model represents a different spatial and temporal scale of the wound healing process. The LDE model provides detailed information about the network state of 91 different signaling nodes in an individual fibroblast, while the ABM predicts fibroblast number, collagen area fraction, and collagen alignment at the tissue level. In the work presented here, we couple these LDE and ABM models in order to capture the dynamic interplay between fibroblast intracellular signaling and spatially heterogeneous extracellular cues such as cytokines and ECM composition, which themselves are modulated by individual fibroblast behaviors. Verification tests confirmed that the coupled model and network model alone exhibit consistent behavior in response to constant cytokine and growth factor inputs, allowing for the establishment of a framework that can readily incorporate updates from either the network model or ABM without affecting the integrity of the individual model predictions. Furthermore, a subset of coupled model predictions was validated by comparison to measurements of pro-collagen 1, αSMA, and F-actin expression in human cardiac fibroblasts treated with combinations of cytokines and growth factors in vitro. We believe this work demonstrates the first coupling of a large-scale network model to predict tissue-level changes in ECM composition in the setting of fibrosis with feedback from environmental cues (e.g., diffusible cytokines) to regulate the signaling of individual cells. Predictions about cytokine and growth factor production from fibroblasts are computed in physical units, which were not previously possible with a logic-based network model alone. This coupled model provides a platform for systematically testing molecular interventions with the ability to measure their effects on single cell signaling and ECM composition with detailed spatial resolution.

Materials and Methods

Description of Individual Models

Agent-Based Model

An agent-based model (ABM) is comprised of value layers and agents (An et al., 2009). The value layers in this two-dimensional ABM represent features of the extracellular space, including collagen, latent TGFβ, active TGFβ, IL-1β, IL-6, and TNFα. All cytokines are stored as concentrations in pg/mL, and collagen is quantified as an area fraction. The value layers are divided into a 10 × 10 grid, where each individual grid space measures 10 × 10 μm. A volume for each grid space is approximated based on cell culture conditions in a 96 well-plate, which is the primary source of experimental data used to inform this model. For soluble cytokines (active TGFβ, IL-1β, IL-6, and TNFα), it is assumed that these cytokines are uniformly distributed in the media above each cell, resulting in a compartment of 10 × 10 × 3125 μm (3.125e-7 mL). Latent TGFβ binds to the extracellular matrix (Horiguchi et al., 2012), and is thus assumed to occupy the space immediately surrounding the cell, or 10 × 10 × 10 μm (1e-9 mL). The individual grid space approximates the footprint of a single fibroblast, allowing the model to simulate a maximum of 100 fibroblasts simultaneously. The total number of fibroblasts is kept relatively low to allow for calculation of the entire network state of each fibroblast while minimizing computational time for the purposes of method development. Value layers store a unique quantity in each grid space that can be modulated by parameters including degradation rates, activation rates, and the agents that move over them. The agents in this model represent cardiac fibroblasts that migrate and modulate their extracellular space by depositing and degrading collagen, and secreting cytokines.

Logic-Based Network Model

The logic-based differential equation (LDE) network model is a previously published (Zeigler et al., 2016a) model of cardiac fibroblast signaling that integrates 10 signaling pathways with 11 biochemical or mechanical stimuli that are important for myofibroblast activation and ECM remodeling. These stimuli include IL-1 (interleukin 1), IL-6 (interleukin 6), TNFα (tissue necrosis factor α), NE (norepinephrine), NP (natriuretic peptide), β-integrins, TGFβ (tissue growth factor β), angiotensin II, PDGF (platelet derived growth factor), ET1 (endothelin 1), mechanical stimulation, and forskolin. The network includes 91 nodes connected by 142 reactions, which are supported by in vitro data collected from cardiac fibroblasts. The network was constructed using a logic-based ordinary differential equation modeling approach, where the activity of each node is modeled using a normalized Hill ODE with default parameters and logic gating. Default reaction parameters include weight (0.9), Hill coefficient (1.4), and EC₅₀ (0.6), and species parameters include y_init(0), y_max(1), and τ. The τ parameter (time constant) was scaled according to the type of reaction: 6 min for signaling reactions, 1 h for transcription reactions, and 10 h for translation reactions. The baseline level of input is defined as 25% activity for all input nodes. The system of ODEs is generated using the Netflux software available at: https://github.com/saucermanlab/Netflux, and implemented in MATLAB.

Coupled Model

Interactions That Drive the Coupled Model

Figure 1 provides an overview of the components and interactions between the LDE network model and ABM. The ABM contains the value layers that represent the extracellular space and the cardiac fibroblasts that migrate over and interact with these value layers. The time step for this coupled model is 1 h, representing the approximate timescale for a change in input to the cell signaling network to affect production of cytokines and ECM proteins that will be deposited in the ABM (Enríquez-de-Salamanca et al., 2008; Azghani et al., 2014). Agents execute a series of methods at each time step: receive input from value layers, update network state, secrete latent TGFβ and IL-6, deposit collagen, migrate. Migration occurs randomly for all simulations, and cell proliferation and death are not simulated. One agent is allowed to occupy an individual grid space, and agent migration is confined to the borders of the simulation space. This series of methods is repeated for 1,008 time steps (6 weeks).

FIGURE 1

Figure 1. Components of individual ABM and network models. The ABM is comprised of agents that store information about attributes and perform methods. Value layers can be modified independently by defined parameters or by the activity of agents. Individual agents store a network state, which is updated by the fibroblast network model.

Interactions between the network model and ABM are described by Equations (1–10). These equations are used to define the behavior at the interface of the two models and are distinct from the equations that define the network model alone. The network model operates using normalized values between 0 and 1, whereas the ABM stores values in terms of physical concentrations. This set of equations act as a translator between these two systems. Figure 2 describes how these equations interact with components of each model and the order in which these methods are executed.

FIGURE 2

Figure 2. Coupled model process diagram. A detailed process diagram illustrates the methods and order in which they occur at each time step (1 h), and components of the ABM and network model that interact. Boxed numbers refer to the equation number which describes that process.

Network Model Inputs

Equations (1–4) are used to translate the cytokine levels stored as concentrations in the value layers of the ABM into inputs for the network model. Input weights for the network model range from 0 to 1, representing receptor activation between 0 and 100%. These weights are determined using the quantitative dissociation constants (Table 1) for the inputs of interest (IL-6, IL-1β, TNFα, and TGFβ). The dissociation constant is the concentration of ligand at which approximately half of the free ligand is bound to receptor at equilibrium. Receptor activity is described by a Hill equation, where a concentration of ligand equal to the K_d is considered 50% activation of the input node. Values for each of these dissociation constants are listed in Table 1.

$\begin{array}{ll}{w}_{IL\text{-}6}=\frac{\left[IL\text{-}{6}_{ABM}\right]}{[IL\text{-}{6}_{ABM}]+K_{d,IL\text{-}6}}& (1)\end{array}$

$\begin{array}{ll}{w}_{IL\text{-}1\beta }=\frac{[IL\text{-}1{\beta }_{ABM}]}{[IL\text{-}1{\beta }_{ABM}]+K_{d,IL\text{-}1\beta }}& (2)\end{array}$

$\begin{array}{ll}{w}_{TNF\alpha }=\frac{\left[TNF{\alpha }_{ABM}\right]}{[TNF{\alpha }_{ABM}]+K_{d,TNF\alpha }}& (3)\end{array}$

$\begin{array}{ll}{w}_{TGF\beta }=\frac{[TGF{\beta }_{ABM}]}{[TGF{\beta }_{ABM}]+K_{d,TGF\beta }}& (4)\end{array}$

TABLE 1

Table 1. Dissociation constants.

Active TGFβ and Latent TGFβ

Equation (5) describes the production of latent TGFβ from sources other than fibroblasts (k_gen), secretion of latent TGFβ from fibroblasts (k_sec) based on the network activity of latent TGFβ (latentTGFβ_net), degradation of latent TGFβ (k_deg), and activation of latent TGFβ (k_act) based on the concentration of latent TGFβ in the ABM (latentTGFβ_ABM). The generation rate (k_gen) describes the production of latent TGFβ from sources that are not currently represented in this model (e.g., macrophages, neutrophils, etc.) and is used to maintain the gradient setup as described below under Initial Conditions. The secretion rate (k_sec) describes the maximum physiological secretion of latent TGFβ from fibroblasts under stimulated conditions and this rate is scaled based on the network activity level (0–1) of latent TGFβ for each fibroblast. The degradation rate is a first-order rate based on the stability of latent TGFβ in vitro, and the activation rate describes the proportion of latent TGFβ that is converted to active TGFβ. Based on literature review, we chose a value for k_{act,latentTGFb} that maintains active TGFβ at 4–5% of total TGFβ, which is consistent with values measured in both in vitro and in vivo studies (Maeda et al., 2002; Hawinkels et al., 2007). In Equation (6), we use a rapid equilibrium assumption for active TGFβ concentration because the degradation rate of active TGFβ is on the order of minutes, much faster than our model time step of 1 h. Thus, we assume that the kinetics of active TGFβ are rate limited by the kinetics of latent TGFβ and come to a rapid quasi-equilibrium based on current latent TGFβ concentrations. Parameter values for Equations (5) and (6) can be found in Table 2.

$\begin{array}{l}\frac{{\partial latentTGF\beta }_{ABM}}{\partial t}={k_{gen,latentTGF\beta }+k}_{sec,latentTGF\beta }\\ *latentTGF{\beta }_{net}\\ -k_{deg,latentTGF\beta }*latentTGF{\beta }_{ABM}\\ -k_{act,latentTGF\beta }*latentTGF{\beta }_{ABM}& (5)\end{array}$

$\begin{array}{llll}TGF{\beta }_{ABM}(t)& =& k_{act,latentTGF\beta }*latentTGF{\beta }_{ABM}(t)& (6)\end{array}$

TABLE 2

Table 2. Parameters for active and latent TGFβ kinetics.

Inflammatory Cytokines

Equations (7–9) describe the production and degradation of IL-1β, IL-6, and TNFα. IL-1β and TNFα are not secreted by the current fibroblast network model, so these equations simply consist of a generation rate and first-order degradation rate. The generation rates are selected to maintain prescribed cytokine gradients as described below under Initial Conditions. The equation for IL-6 kinetics simply has the addition of a secretion rate that represents the maximum physiological secretion of IL-6 from fibroblasts under stimulated conditions and is scaled based on the network activity level (0–1) of IL-6 for each fibroblast. Parameter values for Equations (7–9) can be found in Table 3.

$\begin{array}{llll}\frac{{\partial IL-1\beta }_{ABM}}{\partial t}& =& k_{gen,IL\text{-}1\beta }-k_{deg,IL\text{-}1\beta }*{IL-1\beta }_{ABM}& (7)\end{array}$

$\begin{array}{l}\frac{{\partial IL-6}_{ABM}}{\partial t}=k_{gen,IL\text{-}6}+k_{sec,IL6}*{IL-6}_{net}-\\ k_{deg,IL6}*{IL-6}_{ABM}& (8)\end{array}$

$\begin{array}{llll}\frac{{\partial TNF\alpha }_{ABM}}{\partial t}& =& k_{gen,TNF\alpha }-k_{deg,TNF\alpha }*TNF{\alpha }_{ABM}& (9)\end{array}$

TABLE 3

Table 3. Parameters for inflammatory cytokine kinetics.

Collagen

Equation (10) describes the deposition and degradation of collagen in the ABM based on the collagen I and III mRNA nodes in the network model. Deposition of collagen occurs only where a fibroblast is present and is based on the value of the collagen I and III mRNA nodes in the network model for each fibroblast. Degradation is modeled as a first-order process based on the current collagen concentration in the ABM and thus occurs at every grid location, regardless of the presence of a fibroblast. This assumes evenly distributed MMP activity since we are not explicitly representing MMP production in this model.

These two parameters (Table 4) were fit based on previously published data in a rat model of MI (Fomovsky and Holmes, 2010). Baseline collagen area fraction was considered to be 4%, based on typical measurements from a healthy rat prior to an infarction (MacKenna et al., 1994). Collagen I and III mRNA activity levels from the network model run at a baseline condition for 6 weeks (0.25 for all input nodes) were used to fit the baseline experimental data. This was done by analytically solving equation 10 for the ratio of k_dep/k_deg that produces a steady state collagen area fraction of 4%. Then, the network model was run for 6 weeks to simulate a stimulated condition (0.5 input for TGFβ, IL-1β, IL-6, and TNFα network nodes) and used to fit the infarct experimental data. This was accomplished by performing a parameter sweep of values for k_dep while constraining k_deg to satisfy the ratio determined previously and minimizing the sum of squared error (SSE) between the model fit and infarct experimental data.

$\begin{array}{l}\frac{{\partial Collagen}_{ABM}}{\partial t}=k_{dep,Collagen}(ColIRN{A}_{net}+ColIIIRN{A}_{net})\\ -k_{deg,Collagen}*Collage{n}_{ABM}& (10)\end{array}$

TABLE 4

Table 4. Parameters for collagen deposition and degradation.

Initial Conditions

To evaluate the effects of spatial gradients in fibrotic and inflammatory cues, the value layers are initialized with a gradient of TGFβ increasing from bottom to top, and a gradient of inflammatory cytokines (IL-1β, IL-6, and TNFα) increasing from left to right. Thus, each individual grid space contains a unique combination of fibrotic and inflammatory cues. Cytokine gradients were specified by scaling the generation rate (k_gen) along the x or y axis to result in concentrations ranging up to twice the dissociation constant for that particular cytokine or growth factor, corresponding to receptor activation rates between ~16 and 67% (Figures 3A–C). The purpose is 2-fold: to explore a dynamic range of inputs to the network model and to create an environment where fibroblasts migrate through spatially varying environmental cues.

FIGURE 3

Figure 3. Agent-based model is initialized with cytokine gradients. (A) Four phenotypic regions are created by a combination of fibrotic and inflammatory cues. (B) TGFβ is initialized with an increasing gradient from bottom to top. (C) IL-1β, IL-6, and TNFα are initialized with an increasing gradient from left to right.

Parameter Estimation and Fitting

A total of 17 parameters are defined in this set of 10 equations. Of these parameters, 11 are based on literature and 6 are estimated, or fit to experimental data, as noted in Tables 1–4. The dissociation constants, degradation rates, and production rates from fibroblasts are based on literature. The generation rates are estimated based on a mass balance constraint in order to create the specified initial gradient, which are ultimately based on the dissociation constants found in literature. Coefficients for collagen deposition and degradation are fit based on experimental data.

Sensitivity Analysis

A sensitivity analysis of all 17 model parameters was conducted. Each parameter was decreased individually by an order of magnitude (0.1×) and compared to the results of a model run with all parameters at baseline values. A sensitivity coefficient was calculated using Equation (11), where y_o and y_i are the measured state variable when parameters are at baseline or perturbed, respectively, and p_o and p_i are the values of the baseline parameter and perturbed parameter.

$\begin{array}{ll}S=\frac{y_i-y_o}{p_i-p_o}*\frac{p_o}{y_o}& (11)\end{array}$

The state variables measured in the sensitivity analysis are the total collagen content, global semivariance (r_xy), and semivariance in either the x (r_x) or y (r_y) dimension. Collagen content is measured by summing collagen area fraction across all individual grid spaces. Global semivariance is (r_xy) defined as:

$\begin{array}{ll}{r}_{xy}=\frac{1}{2s_0}{\displaystyle \sum _i}{\displaystyle \sum _j}{W}_{ij}{(x_i-x_j)}^2& (12)\\ s_0={\displaystyle \sum _i}{\displaystyle \sum _j}{W}_{ij}\end{array}$

where x_i is an observed data point, x_j is an adjacent observation, W_ij is a matrix of spatial weights, and s₀ is the sum of all W_ij. If two data points are immediate neighbors, W_ij is assigned as 1, otherwise W_ij is set to 0 (Lee et al., 2019). Semivariance in the x dimension is calculated by assigning 1 to W_ij for adjacent observations in the x direction and 0 otherwise, and vice versa for the y dimension.

Model Implementation

This model was implemented using Repast Simphony 2.6 with a java engine to connect to MATLAB R2018b, which was used to run the network model and perform all data analysis. All simulations here were performed on a single CPU (Intel^® Xeon^® E5-2640 v4@2.4GHz). The approximate runtime to simulate 100 fibroblasts for a period of 6 weeks is 1 h and 17 min.

Cardiac Fibroblast in vitro Experiments

Primary human ventricular cardiac fibroblasts were purchased from PromoCell (PromoCell C-12375; PromoCell GmbH, Germany). Cells were cultured in DMEM containing 10% FBS and 1% Pen/Strep, and were kept in an incubator maintained at 5% CO₂. Cells were plated in a 96-well plate at 5,000 cells/well and then grown in 10% FBS for 24 h, serum starved for 24 h, and then treated with the following conditions for 96 h: 0%FBS control media, 0%FBS media with 20 ng/mL TGFβ1 (Cell Signaling Technology, 8915LC), and 0% FBS media with 10 ng/mL human IL1β (Cell Signaling Technology, 8900SC). Cells were then fixed in 4% PFA in PBS for 30 min, permeabilized and blocked for 1 h in a solution containing 3% BSA and 0.2% Triton, and then stained overnight at 4°C with a 1:500 primary Anti-Collagen I antibody (Abcam, ab34710). After an overnight incubation, cells were washed 3 × in PBS and stained with 1:5,000 Dapi, 1:1,000 Phalloidin CruzFluor 647 stain (Santa Cruz Biotechnology, sc-363797), 1:250 α-Smooth Muscle Actin preconjugated antibody (Sigma-Aldrich, C6198), and 1:1000 Goat-anti-Rabbit (secondary for Anti-Collagen I) (ThermoFisher Scientific, A-11034).

Microscopy and Image Analysis

96-well-plates we imaged using the Operetta CLS High-Content Analysis System with confocal view, and Dapi, Alexa 488, TRITC, and Alexa 647 imaging channels (Perkin Elmer). Three wells for each condition were imaged and quantified. To identify individual cells, an automated image analysis pipeline was employed in CellProfiler (Carpenter et al., 2006). Fibroblast nuclei were identified by DAPI signal, and fibroblast boundaries corresponding to each nuclei were segmented based on collagen and phalloidin (actin) signals using the “propagate” algorithm. αSMA, pro-collagen I, and phalloidin signals were integrated within each cell's boundary to determine fluorescence per cell. To reduce error from edge effects, only cells in the center tile of each well were measured. The median fluorescence for all cells in a given well was reported (n = 3 replicate wells per treatment group, 250–450 cells per well were imaged). Significance between groups was determined by one-way ANOVA with Tukey HSD post-hoc test, p ≤ 0.05 considered significant. Data are provided in Supplemental Table 1.

Results

Coupled Model Can Reproduce Predictions Made by Network Model Alone

Verification tests were performed to evaluate whether coupling of LDE and ABM models affected results obtained from each model individually. This was accomplished by seeding one fibroblast in each grid space with no migration and simulating either an unstimulated condition (0.25 input for all nodes) or stimulated condition (0.5 input for TGFβ, IL-1β, IL-6, and TNFα). We compared the activity level of all 91 nodes of the fibroblast network state at steady state for each condition and calculated the sum of squared error (SSE) between the network-only and coupled models. The SSE for the unstimulated condition is 3.865e-7 (Figure 4A) and 1.168e-6 for the stimulated condition (Figure 4B). We next examined the SSE between the network states of the coupled model and network model alone for all 100 combinations of cytokine inputs from Figure 3, resulting in SSEs in the range of [2.07e-8–8.88e-5]. Thus, the network-only and coupled models produce equivalent network states for non-migrating fibroblasts when cytokine inputs are maintained at a constant level for individual cells.

FIGURE 4

Figure 4. Verification tests confirm that coupled model and network model produce equivalent fibroblast network states. (A) In an unstimulated condition (0.25 input for all nodes), and a (B) stimulated condition (0.5 input for TGFβ, TNFα, IL-1β, and IL-6 nodes), the network state of a fibroblast using the coupled model or network model alone are comparable with a SSE of 3.865e-7 (A) and 1.168e-6 (B).

Coupled Model Predicts That Inflammatory Cytokines Antagonize TGFβ-Induced Collagen Accumulation

Model parameters for collagen deposition and degradation were fit to match experimental data obtained previously from a rat model of MI (Fomovsky and Holmes, 2010). For these simulations, a single fibroblast was placed in every grid space and not allowed to migrate. Model simulations were initialized to match the measured rise in collagen area fraction in healing infarcts when fibrotic and inflammatory inputs to the coupled model (TGFβ, IL-1β, IL-6, and TNFα) were maintained at an elevated level of 0.5, and to match the normal myocardial collagen area fraction when the same inputs were maintained at their baseline values of 0.25 (Figure 5A). Expanding to a broader range of cytokine combinations (100 combinations of fibrotic vs. inflammatory cytokines), the coupled model predicted biologically plausible variations in steady-state collagen content. As shown in Figure 5B, collagen content is highest in areas with high TGFβ input and low inflammatory input, and is reduced as inflammatory input increases for the same magnitude of TGFβ input.

FIGURE 5

Figure 5. Coupled model predicts collagen profile over a range of physiological conditions. (A) Collagen area fraction for an unstimulated condition (0.25 input for all nodes) is compared to baseline collagen area fraction (4%) in a healthy rat. Model predictions for a stimulated condition (0.5 input for TGFβ, TNFα, IL-1β, and IL-6 nodes) are compared to results from a rat model of myocardial infarction up to 6 weeks post-MI. Error bars = SEM. (B) Collagen area fraction predictions at 6 weeks from a model simulation with gradient initial conditions and a fibroblast in each grid space (n = 100).

Crosstalk Between TGFβ and Inflammatory Extracellular Cues Produces Complex Signaling Behaviors

The fibroblast network model exhibits complex behaviors due to its integration of 10 interdependent signaling pathways. Figure 6 illustrates representative nodes in the fibroblast network and their activity level under conditions of constant cytokine inputs for many combinations of inflammatory and fibrotic inputs, as described previously. Network receptors display a range of activation patterns based on their extracellular cues (Figures 6A–C). IL-1β receptor activation closely follows the gradient created by the initial conditions. In contrast, TGFβR1 and endothelin-1 (ET-1) are influenced by autocrine feedback loops. Inflammatory cytokines cause inhibition of TGFβR1 that increases along the x-axis as the concentration of inflammatory cytokines increase. TGFβR1 activation is also influenced by the rate of latent TGFβ activation, which occurs in the ABM value layers. Fibroblasts secrete latent TGFβ, which is then activated to active TGFβ. But as noted in Figure 6G there is minimal latent TGFβ produced in environments of low TGFβ and inflammatory input. As a result, there is decreased TGFβR1 activation in this quadrant. In contrast to the gradual applied input gradients, ET-1 receptor displays switch-like activation, due to autocrine feedback of activator protein 1 (AP1) downstream of both TGFβ and IL-1β.

FIGURE 6

Figure 6. Signaling network exhibits a range of activation patterns in response to extracellular cues. Node activity level of individual network nodes for each fibroblast at steady state (6 weeks). Model simulation with gradient initial conditions and a fibroblast in each grid space (n = 100). Heat maps show network states for input receptors (A–C), intermediate network nodes (D–F), and network outputs (G–I).

Some nodes downstream of each of these inputs display similar activation patterns (Figures 6D–F), while others integrate multiple unique upstream inputs. Smad7, which is immediately downstream of STAT and IL-1β receptor, displays an activity pattern similar to the IL-1β receptor. Smad3 is downstream of TGFβR1, and it regulates many network outputs including collagen mRNA, fibronectin, periostin, and αSMA. NF-κB activity is regulated by many inputs, including IL-1βR, ERK, p38, and AKT. But ERK and p38 (which are immediately downstream of the ET-1 receptor) dominate the response of NF-κB, so it displays an activation pattern most similar to ET-1 receptor. NF-κB contributes to the expression of MMPs, fibronectin, and provides feedback to IL-6 input.

Network outputs represent the integration of many upstream inputs (Figures 6G–I). Latent TGFβ expression is primarily influenced by AP1 transcriptional activity, which itself is regulated by ERK and JNK. The model predicts that IL-1β antagonizes TGFβ-induced collagen I mRNA and αSMA mRNA, which is validated by experimental studies in lung and dermal fibroblasts (Mia et al., 2014). Expression of collagen I mRNA is predicted to be a product of input from Smad3, SRF (serum response factor), and CBP (CREB binding protein). proMMP 1 expression is a prime example of integration of multiple upstream inputs that each exhibit distinct activation patterns including AP1, Smad3, and NF-κB. Visualization of how this combination of transcription factors regulates proMMP 1 expression is shown in Supplementary Figure 1. In summary, the coupled model provides a platform to investigate how combinations of dynamic inputs affect downstream intermediate network nodes and network outputs.

Key Parameters Affect Spatial Gradient of Collagen Deposition

A sensitivity analysis was conducted to determine the relative influence of decreasing the values of parameters associated with ABM-network coupling on overall collagen content (area fraction) and collagen heterogeneity (semivariance, either globally, or in the x or y dimension). Parameters were individually decreased by an order of magnitude (0.1×), and normalized sensitivity coefficients were computed, in which positive coefficients indicate positive correlation of the parameter with the output measured (see Equation 11). Parameters were ranked by their positive influence on collagen area fraction (Figure 7A). Parameters related to TGFβ production, activation, and degradation are the most influential in determining collagen content, because TGFβ input is important in altering downstream collagen I and III mRNA activity in the fibroblast signaling network. As expected, the coupled model is also highly sensitive to the two parameters in Equation (10) that govern deposition and degradation in the collagen layer. In terms of inflammatory inputs, this analysis reveals that IL-1β and IL-6 are more influential in determining the collagen profile than TNFα input. This is likely because IL-6 has a downstream effect on Smad3, which promotes collagen mRNA activity, and IL-1β upregulates NF-κB, which has a positive feedback on IL-6. TNFα has a smaller effect on NF-κB signaling and no direct connection to Smad3 signaling. It is also noted that the two parameters related to secretion of latent TGFβ and IL-6 (k_{sec,latentTGFβ} and k_sec,IL−6) from the fibroblast have little effect on the overall collagen profile, yet TGFβ and IL-6 inputs themselves seem to be very influential. This is likely because the production rates from fibroblasts are not high enough to significantly impact the gradients that are created in the initial conditions.

FIGURE 7

Figure 7. Key parameters affect spatial gradient of collagen deposition. (A) Sensitivity coefficients calculated based on Equation (11) with individual parameter perturbations of 0.1×. State variable outputs include total collagen area fraction, global semivariance, semivariance in the x direction, and semivariance in the y direction. Collagen area fraction heat maps at 6 weeks with (B) all parameters at baseline, (C) k_gen,IL−1β parameter multiplied by 0.1, and (D) K_d,TGFβ parameter multiplied by 0.1.

Interestingly, some parameters may have different effects on overall collagen content and collagen spatial heterogeneity. Decreasing parameters such as the K_d of TGFβ for its receptor, increase both collagen content and collagen heterogeneity (Figures 7A,D), compared to a collagen profile at 6 weeks when all parameters are at baseline values (Figure 7B). Decreasing parameters associated with synthesis of latent TGFβ or collagen, or degradation of IL-6, cause a decrease in both collagen content and heterogeneity in both dimensions. In contrast, decreasing the IL-1 generation term (k_gen,IL−1) has little effect on the total collagen content, but has opposite effects on collagen heterogeneity in the x and y dimension, as measured by semivariance in either direction (Figure 7C).

Single Cell Dynamics in Response to a Changing Extracellular Environment

To test the role of fibroblast migration on collagen remodeling, fibroblasts were seeded sparsely in the coupled model and allowed to migrate stochastically at a rate of one grid space per hour. Fibroblasts were seeded at moderate density (20 fibroblasts) within the fibrotic vs. inflammatory cytokine grid and responses simulated for 6 weeks. Fibroblasts experience changes in their extracellular environment as they migrate, which causes their intracellular signaling network state and rate of collagen deposition to change accordingly. Figure 8 illustrates single cell migration trajectories (panels A and C), local cytokine inputs, and gene expression (panels B and D) for two representative fibroblasts migrating within the cytokine gradient environment. The fibroblast shown in panels A and B remains in areas with high to moderate TGFβ levels and with increasing levels of inflammatory cytokines. Correspondingly, this fibroblast exhibited relatively high levels of collagen mRNA expression that mirrored the level of TGFβ input. As this fibroblast migrated to regions of increasing IL-6, there was a delayed but then rapid increase in MMP mRNA expression, consistent with switch-like responses seen in Figure 6I.

FIGURE 8

Figure 8. Individual fibroblasts respond dynamically to extracellular environment. (A,C) Fibroblast migration path for a single fibroblast over a period of 50 h. Fibroblast starting location indicated by filled black triangle and end location indicated by open white triangle. (B,D) Corresponding TGFβ and IL-6 inputs for the fibroblasts tracked in (A,C), and their respective collagen and MMP mRNA activity over the time course of 50 h.

Figures 8C,D track a separate fibroblast that remains in areas of low to moderate TGFβ levels, but migrates from a region of high to low inflammatory inputs. This simulation shows similarly that collagen mRNA expression closely follows TGFβ inputs. Interestingly, it also demonstrates that exposure to high IL-6 levels triggers a rapid and sustained increase in MMP mRNA expression, which persists well after the cell migrates to a region with lower levels of inflammatory cytokines. Thus, some network nodes respond with close coordination to particular cytokine inputs, whereas other nodes may be activated in a switch-like manner consistent with the activation patterns seen in Figure 6. Supplementary Videos 1, 2 offer a visual of how the entire network state changes over time for the individual fibroblasts presented in Figure 8.

Fibroblast Migration Speed and Density Affect Collagen Spatial Heterogeneity

Under normal conditions, fibroblasts migrate at a speed of ~10 μm/h, but this can vary significantly in the presence of growth factors and cytokines (Ware et al., 1998; Pérez-Rodríguez et al., 2018). As noted above, individual grid spaces are 10 × 10 μm, so the baseline migration speed was set at 1 grid space per hour. To test the impact of altered migration speed, migration speeds were set to default values (1 grid/h), decreased (1 grid space/10 h) or increased (10 grid spaces/h). Slower migration speed resulted in greater heterogeneity in the collagen profile, while faster migration speed resulted in a more homogenous collagen profile (Figures 9A–C). In contrast, overall collagen content was linearly dependent on fibroblast density, but not migration speed or initial cell location (n = 10 simulations per condition) (Figure 9D). As with the sensitivity analyses, heterogeneity in collagen was quantified by global, x-direction, or y-direction semivariance (see Equation 12). Consistent with qualitative observations from Figures 9A–C, decreasing the migration speed enhanced both the average magnitude and run-to-run variance in collagen heterogeneity (Figure 9E). Faster migration decreased collagen heterogeneity globally and in the y dimension, but not in the x dimension. The coupled model predicts that fibroblast migration can have a substantial impact on the spatial heterogeneity and stochasticity of collagen deposition.

FIGURE 9

Figure 9. Fibroblast migration speed and density affect spatial heterogeneity of collagen. Collagen area fraction heat map at 6 weeks for simulations with 20 randomly migrating fibroblasts with a migration speed of (A) 1 grid/10 h, (B) 1 grid/h, and (C) 10 grids/h. (D) Average collagen area fraction at 6 weeks for each migration speed and simulations with 20, 40, 60, 80, or 100 fibroblasts. Mean reported for 5 runs. Error bars = standard deviation. (E) Semivariance calculated globally, in the x direction, and y direction for each migration speed. Mean reported for 10 runs. Error bars = standard deviation. ^*p < 0.05.

A Subset of Coupled Model Predictions Was Validated by Comparison to in vitro Experiments

In order to perform experimental validation of our coupled model predictions, we ran a series of simulations wherein input cytokine (IL-1β) and growth factor (TGFβ1) levels were varied from baseline to simulate those tested with in vitro experiments using primary human cardiac fibroblasts (HCFs). HCFs were treated with either TGFβ1 (20 ng/mL), IL-1β (10 ng/mL), or TGFβ1 (20 ng/mL) + IL-1β (10 ng/mL), and compared to a control condition in media without FBS (since this is also a source of TGFβ1). Pro-collagen I, αSMA, and F-actin expression were quantified using immunocytochemistry and image processing to quantify the median fluorescence in each of these experimental conditions (Figure 10A). Treatment with TGFβ1 significantly increased expression of pro-collagen I, αSMA, and F-actin compared to the control condition, while IL-1β treatment alone had no significant effect on pro-collagen I, αSMA, or F-actin expression when compared to the control condition. The combination of TGFβ1 and IL-1β treatment decreased expression of pro-collagen I, αSMA, and F-actin when compared to TGFβ1-only treatment, and was statistically significant in the cases of pro-collagen I and F-actin. Representative images of pro-collagen 1 (green), αSMA (orange), and F-actin (blue) expression in HCFs for each of these treatment conditions indicate the trends described above (Figure 10C). These experimental measurements were compared to in silico predictions that simulated the addition of these factors at the same concentrations tested experimentally: TGFβ1 (20 ng/mL), IL-1β (10 ng/mL), or TGFβ1 (20 ng/mL) + IL-1β (10 ng/mL). As with the experimental results, predictions were compared to a control simulation in which all parameters were set to baseline levels (Figure 10B). Most model predictions qualitatively agreed with the trends observed with in vitro experiments. For example, network expression of collagen I mRNA, αSMA, and F-actin were increased relative to the control simulation in response to TGFβ1 stimulation. Similar to observed experimental results, IL-1β treatment alone had no effect on collagen I mRNA, αSMA, or F-actin network expression compared to the control simulation, but the combination of TGFβ1 and IL-1β treatment was predicted to decrease the expression of collagen I mRNA and αSMA when compared to TGFβ1 only treatment. Unlike experimental results, however, simulating this combined treatment predicted no change in F-actin network expression when compared to simulating TGFβ1 only treatment. These model predictions were further validated by published experimental studies wherein IL-1β attenuated TGFβ1-induced collagen I synthesis and αSMA expression of lung and dermal fibroblasts (Mia et al., 2014).

FIGURE 10

Figure 10. Coupled model predictions were compared to independent in vitro experiments using human cardiac fibroblasts treated with TGFβ1 and/or IL-1β. (A) Pro-collagen 1, αSMA, and F-actin expression from in vitro experiments with human cardiac fibroblasts were quantified by image analysis to measure the median fluorescence for all individual cells in each well (n = 3). Treatment conditions included control, TGFβ1 (20 ng/mL), IL-1β (10 ng/mL), and TGFβ1 (20 ng/mL) + IL-1β (10 ng/mL). Error bars = standard deviation. ^*p < 0.05 with reference to control condition. ^∧p < 0.05 with reference to TGFβ1 condition. (B) Coupled model predicts network expression of collagen I mRNA, αSMA, and F-actin when simulating the addition of TGFβ1 (20 ng/mL), IL-1β (10 ng/mL), and TGFβ1 (20 ng/mL) + IL-1β (10 ng/mL), compared to a simulation with all parameters at baseline. (C) Representative images of human cardiac fibroblast expression of pro-collagen 1 (green), αSMA (orange), and F-actin (purple) when treated with TGFβ1 (20 ng/mL), IL-1β (10 ng/mL), and TGFβ1 (20 ng/mL) +IL-1β (10 ng/mL), compared to control. Nuclei are stained with DAPI (blue). Scale bar = 500 microns.

Discussion

A Novel Hybrid Multiscale Model of Tissue Fibrosis

Here we present a novel hybrid multiscale model of tissue fibrosis that couples a LDE model of cardiac fibroblast intracellular signaling with an ABM of multi-cellular tissue remodeling. Prescribed gradients of inflammatory cues (IL-1β, IL-6, and TNFα) and fibrotic cues (TGFβ) stimulate migrating fibroblasts to respond dynamically to their locally varying extracellular environment. Under conditions with no fibroblast migration and constant cytokine input levels, the coupled model was verified to exhibit consistent network states predicted by the network model alone (Figure 4). In contrast, the addition of fibroblast migration across a gradient of cytokine inputs demonstrated that fibroblasts respond dynamically to both their local cytokine environment and their previous history of cytokine exposure. Spatial heterogeneity of collagen was dependent on the speed of fibroblast migration and key parameters (e.g., Il-1β synthesis rate) identified in a sensitivity analysis as having distinct effects on semivariance in the x or y dimension. Additionally, several parameters were identified to be influential in contributing to the overall amount of collagen deposition, including cell density and model parameters related to TGFβ production and activation.

The effects of pro-fibrotic stimuli (TGFβ) on increasing collagen expression and other myofibroblast markers is well-established (Reed et al., 1994; Chen et al., 2000; Thannickal et al., 2003). However, the crosstalk of IL-1β with other fibrotic signaling pathways is not as well-described. This coupled model provides a framework for investigating the effects of combined inflammatory and pro-fibrotic cues on spatial fibroblast activity and ECM composition. A subset of coupled model predictions was validated by comparison to experiments with human cardiac fibroblasts treated with combinations of TGFβ1 and IL-1β in vitro. The model accurately predicted TGFβ1-enhanced expression pro-collagen I, αSMA, and F-actin, as well as negative crosstalk on pro-collagen I and αSMA by IL-1β. However, the model did not predict the experimentally-observed attenuation of TGFβ1-enhanced F-actin expression by IL-1β, suggesting that additional cross-talk mechanisms may need to be explored in future experiments and model revisions.

Previous mathematical models of fibrosis have used deterministic, continuum methods to holistically represent the complex processes of fibrosis (Hao et al., 2015; Friedman and Hao, 2017). For example, Hao et al., describe a model of liver fibrosis using a system of 24 partial differential equations (PDEs) that represent many different cell types, cytokines and growth factors, and interactions between cells, and then use this model to interrogate different treatment options (Friedman and Hao, 2017). Other continuum-based models have focused more on specific mechanisms that contribute to the progression of fibrosis, such as macrophage activation and polarization (Jin et al., 2011; Wang et al., 2012). A model developed by Wang et al., for example, explores how the timing of monocyte recruitment and macrophage differentiation affects left ventricular remodeling following MI (Wang et al., 2012). However, an interesting study by Figueredo et al. suggested that stochastic differential equation approaches that assume continuous space and time could not capture the individual variability and spatial heterogeneity predicted by an agent-based modeling approach applied to the same biological case study, and that emergent behavior of the ABM contributed additional insight about the system (Figueredo et al., 2014).

An increasing number of hybrid models couple continuum with discrete approaches. These hybrid models typically couple ABMs, which use a discrete representation of 2D space or 3D volumes, with continuum based approaches that represent cytokine gradients and/or receptor-ligand kinetics (Warsinske et al., 2016, 2017; Virgilio et al., 2018). For example, Warsinske et al. simulated granuloma-associated fibrosis by incorporating a system of ODEs and PDEs that describe molecular level diffusion of chemokines (TGFβ and prostaglandin) and receptor ligand signaling, coupled with discrete cellular agents whose behaviors were defined by a set of rules that related receptor activation levels to cell proliferation, differentiation, chemotaxis, and secretion of ECM proteins. In these hybrid continuum-ABM models, outcomes at the tissue scale are the emergent product of actions of the individual agents governed by rules that are informed by molecular scale interactions simulated using continuum assumptions (Warsinske et al., 2017).

The hybrid multiscale model presented here represents the coupling of a large-scale intracellular network model, comprising 10 cytokine/neurohormonal inputs and 134 reactions, with an ABM that maps physiologically relevant in vitro concentrations of cytokines and ECM components to normalized network activity levels and vice versa. We believe that this represents the first coupling of a large-scale network model to make predictions about tissue-level changes in extracellular matrix composition in the setting of fibrosis. This coupled model and its use of concentration scaling between the logic-based model and physical units enables the quantitative prediction of fibroblast production of cytokines and growth factors and spatial gradients of cytokine concentrations, which was not previously possible with the network model alone. This coupled model framework will ultimately enable quantitative comparisons of model predictions to in vivo experimental data such as measurements of multiple cytokine concentrations over time, spatial profiles and gradients of ECM components, cell densities, and single cell mRNA expression.

Impact of Spatially Varying Environmental Cues on Fibroblast Signaling

The response of a complex signaling network to multiple simultaneous cues is rarely intuitive, and we have demonstrated that individual nodes of the signaling network respond with distinct patterns of activation (Figure 6). Some receptors respond in sync with their input, such as IL-1βR, whose activity level mimics the gradient initial conditions of IL-1β input. Meanwhile other receptors, such as TGFβR1, display a more complex pattern of activation reflecting not only the gradient initial conditions but also feedback from latent TGFβ activation and inhibition by IL-1 activity. Intermediate nodes often display a similar pattern of activation to their immediate upstream receptors (e.g., Smad7 and IL-1βR, Smad3 and TGFβR1, NF-κB and Endothelin1-R), while network outputs integrate the effects of many upstream network nodes that represent a combination of stimulatory and inhibitory inputs. MMP1, for example, is upregulated by NF-κB and AP1 (activator protein 1) activity, and inhibited by Smad3 activity. Tracking the response of individual fibroblasts moving through varying levels of inflammatory and fibrotic inputs revealed a complex kinetic relationship between the locally sensed extracellular environment and network state of a migrating fibroblast (Figure 8 and Supplementary Videos 1, 2). For example, a fibroblast that experiences a high inflammatory context will upregulate its MMP activity, which remains elevated even if the fibroblast moves to an environment with low inflammatory and fibrotic inputs. The fibroblast's network state is highly dependent on the current extracellular environment in some cases (e.g., collagen mRNA expression in response to TGFβ input) but displays history-dependence of previous environments in other cases.

Processes That Contribute to Spatial Heterogeneity of Collagen Deposition

The artificially prescribed cytokine gradient environment employed in these simulations (Figure 3) was not intended to represent a particular in vivo situation, but was used to evaluate the ability of the coupled signaling and multicellular model to predict the progression of fibrosis across a wide range of signaling contexts. Thus, changes in heterogeneity discussed here reflect the range of responses a population of fibroblasts would be expected to generate across those varied signaling contexts (Figure 5B). For example, we found fibroblast migration speed to be an important determinant of collagen heterogeneity in our simulations (Figure 9). Slower migration speed leads to pockets of high collagen deposition and overall higher heterogeneity. Faster migration produces a more uniform collagen distribution. In healing wounds where cytokine concentrations vary in both space and time, we expect that high migration speeds could similarly blur the effects of variable cytokine levels while slow migration speeds could accentuate them. In contrast, migration speed did not substantially affect average collagen accumulation across the entire simulated range of cytokine combinations. Rather, overall collagen accumulation was strongly dependent on fibroblast density. In addition to fibroblast density, model parameters related to TGFβ production, activation, and degradation are among the most important model parameters in determining total collagen content as well as the gradient of collagen deposition in either dimension (Figure 7), which agrees with the findings from similar models of fibrosis (Warsinske et al., 2016). Other parameters such as the degradation or synthesis of IL-1β had opposite effects on collagen heterogeneity in two dimensions. One advantage of coupling an ABM is that it produces stochastic predictions as a result of individual-based rule sets and a spatial context (Supplementary Figure 2). Repeated runs of the coupled model may help to capture individual variability of spatial fibrosis seen in animal models.

Computational Requirements for Scaling Up

Simulations were performed on a single CPU (Intel^® Xeon^® E5-2640 v4 @2.4GHz). The runtime for 100 fibroblasts for a period of 6 weeks with access to 10 cores is 1 h and 17 min. This model is structured such that parallel computing can be implemented across multiple CPUs to reduce simulation runtime. The limit of computational efficiency is reached when the number of cores is equal to the number of fibroblasts in a simulation, allowing the network state of each fibroblast to be updated simultaneously at each time step. We anticipate that with access to 100 cores, this model with 100 fibroblasts for 6 weeks of simulated time would be computed in ~8 min. Simulation on a high-performance computing system with thousands of cores is expected to enable simulation of up to 100,000 fibroblasts (comparable to a myocardial infarct).

Limitations and Sources of Error

The prescribed input cytokine gradients employed here were used to explore the dynamic range of the network model and create an environment where fibroblasts migrate through a rapidly changing extracellular environment. However, this environment is not representative of a specific physiological environment. Cytokine diffusion was not enabled in the current simulations, in order to maintain the prescribed cytokine input gradients over small spatial area (100 × 100 μm). Furthermore, the current model did not include proliferation and apoptosis, which have been simulated in previous work (Warsinske et al., 2016, 2017). Future applications will incorporate cell migration and proliferation rates that are driven by the dynamic network state of individual fibroblasts (Bailey et al., 2007, 2009). Additionally, this model focused specifically on the contributions of fibroblasts in the progression of fibrosis, through the coupling of a fibroblast signaling network, but future work will incorporate inflammatory cells that serve as local sources of many of the inflammatory cytokines that affect fibroblast signaling (Virgilio et al., 2018).

State of the Multiscale Modeling Field and Contributions of This Multiscale Model

A primary focus of the computational modeling community is to develop methods for integrating biological data across spatial, temporal, and functional scales (Walpole et al., 2013). Continued advancements in the capabilities and availability of high-performance computing has allowed models to tackle more complex problems with greater resolution. Perturbations to fine-grained parameters, such as protein or gene expression data, can predict observable changes to coarse-grained parameters (e.g., cell distributions, tissue patterning) (Stern et al., 2012; Walpole et al., 2017; Warsinske et al., 2017). Furthermore, the use of multiscale models allows for in silico predictions for a wide range of parameter values in a high-throughput manner that would otherwise not be feasible with experimental assays, either due to time or cost constraints, or lack of the appropriate technology. For example, the multiscale model presented here allows for real time tracking of individual fibroblasts and continuous measurements of their network states, which would not be feasible in vivo. While models cannot fully replace experimental studies, they can offer insight into unexpected predictions that can then be experimentally tested or lead to new hypotheses entirely, as demonstrated by Martin et al., who predicted a new therapeutic approach as a result of their in silico experiments of muscle regeneration following injury (Martin et al., 2017).

In summary, we have contributed a hybrid multiscale model of tissue fibrosis by coupling models across spatial and temporal scales. This represents the coupling of a large-scale network model with an ABM to make predictions about fibroblast production of cytokines and growth factors and tissue-level changes in ECM composition. This coupled model makes predictions about fibroblast production of cytokines and growth factors in physical units, which was not possible previously with the logic-based model alone. Verification tests confirmed that the model coupling did not disrupt the behavior of the individual models, allowing for future model revisions or software implementations of individual modules. Application of this coupled model in the context of post-MI wound healing will allow for further investigation and validation of cytokine concentrations, collagen content and heterogeneity, and cell behaviors with both fine spatial and temporal resolution. Experimental studies suggest that collagen density alone may have effects on fibroblast behavior, including adhesion, migration, and gene expression (Loftis et al., 2003), and that furthermore, collagen density and fibroblast density play an important role in the mechanical properties of the myocardium (Fomovsky et al., 2012; van Spreeuwel et al., 2017). This multiscale model framework allows for further investigation and understanding of emergent phenomena that result from the dynamic interplay between molecular signaling, cell behavior, ECM composition, and tissue mechanics. For example, previous computational models have demonstrated that simultaneous targeting of multiple cells types rather than fibroblasts alone can enhance the efficacy of therapies for pulmonary fibrosis (Warsinske et al., 2016). Inflammatory cells, including macrophages and neutrophils, will be incorporated into a model of post-MI wound healing as the primary source and modulators of inflammatory cytokines and TGFβ input, as has been demonstrated previously in simulations of skeletal muscle and lung fibrosis (Martin et al., 2016; Warsinske et al., 2016; Virgilio et al., 2018). This will add another layer of complexity to the spatial heterogeneity of the coupled model by representing cytokine production from individual cells, diffusion of soluble cytokines and growth factors, and migration that is driven by chemokine gradients. Our goal is to develop a hybrid multiscale model that can systematically screen the effect of therapeutic interventions on the progression of cardiac fibrosis, from the level cell signaling to a tissue level of ECM remodeling, with both spatial and temporal resolution.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Author Contributions

SR and TA developed the code for the coupled model. SR, AN, J-JL, SP, JH, and JS contributed to the development of model parameters and rule sets. SC and AN performed in vitro experiments and image analysis. SR wrote the initial draft of the manuscript. JH, SP, and JS contributed to manuscript revision. All authors read and approved the submitted version of the manuscript.

Funding

This work was supported by the National Institutes of Health (HL116449, HL137755, HL137100, and AR069393), the National Science Foundation (1252854 and 1842490), and the Allen Discovery Center at Stanford University.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank Angela Zeigler for the contribution with the fibroblast network model.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2019.01481/full#supplementary-material

References

An, G., Mi, Q., Dutta-Moscato, J., and Vodovotz, Y. (2009). Agent-based models in translational systems biology. Wiley Interdiscip. Rev. Syst. Biol. Med. 1, 159–171. doi: 10.1002/wsbm.45

PubMed Abstract | CrossRef Full Text | Google Scholar

Ancey, C., Corbi, P., Froger, J., Delwail, A., Wijdenes, J., Gascan, H., et al. (2002). Secretion of IL-6, IL-11 and LIF by human cardiomyocytes in primary culture. Cytokine 18, 199–205. doi: 10.1006/cyto.2002.1033

PubMed Abstract | CrossRef Full Text | Google Scholar

Azghani, A. O., Neal, K., Idell, S., Amaro, R., Baker, J. W., Omri, A., et al. (2014). Mechanism of fibroblast inflammatory responses to Pseudomonas aeruginosa elastase. Microbiology 160, 547–555. doi: 10.1099/mic.0.075325-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Bailey, A. M., Lawrence, M. B., Shang, H., Katz, A. J., and Peirce, S. M. (2009). Agent-based model of therapeutic adipose-derived stromal cell trafficking during ischemia predicts ability to roll on P-selectin. PLoS Comput. Biol. 5:e1000294. doi: 10.1371/journal.pcbi.1000294

PubMed Abstract | CrossRef Full Text | Google Scholar

Bailey, A. M., Thorne, B. C., and Peirce, S. M. (2007). Multi-cell agent-based simulation of the microvasculature to study the dynamics of circulating inflammatory cell trafficking. Ann. Biomed. Eng. 35, 916–936. doi: 10.1007/s10439-007-9266-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Baran, P., Hansen, S., Waetzig, G. H., Akbarzadeh, M., Lamertz, L., Huber, H. J., et al. (2018). The balance of interleukin (IL)-6, IL-6.soluble IL-6 receptor (sIL-6R), and IL-6.sIL-6R.sgp130 complexes allows simultaneous classic and trans-signaling. J. Biol. Chem. 293, 6762–6775. doi: 10.1074/jbc.RA117.001163

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamin, E. J., Virani, S. S., Callaway, C. W., Chamberlain, A. M., Chang, A. R., Cheng, S., et al. (2018). Heart disease and stroke statistics-2018 update: a report from the American Heart Association. Circulation 137, e67–e492. doi: 10.1161/CIR.0000000000000558

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolívar, S., Santana, R., Ayala, P., Landaeta, R., Boza, P., Humeres, C., et al. (2017). Lipopolysaccharide activates toll-like receptor 4 and prevents cardiac fibroblast-to-myofibroblast differentiation. Cardiovasc. Toxicol. 17, 458–470. doi: 10.1007/s12012-017-9404-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, E. J., Kloner, R. A., Schoen, F. J., Hammerman, H., Hale, S., and Braunwald, E. (1983). Scar thinning due to ibuprofen administration after experimental myocardial infarction. Am. J. Cardiol. 51, 877–883. doi: 10.1016/S0002-9149(83)80148-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Campaner, A. B., Ferreira, L. M., Gragnani, A., Bruder, J. M., Cusick, J. L., and Morgan, J. R. (2006). Upregulation of TGF-beta1 expression may be necessary but is not sufficient for excessive scarring. J. Invest. Dermatol. 126, 1168–1176. doi: 10.1038/sj.jid.5700200

CrossRef Full Text | Google Scholar

Campbell, S. E., and Katwa, L. C. (1997). Angiotensin II stimulated expression of transforming growth factor-beta1 in cardiac fibroblasts and myofibroblasts. J. Mol. Cell. Cardiol. 29, 1947–1958. doi: 10.1006/jmcc.1997.0435

PubMed Abstract | CrossRef Full Text | Google Scholar

Carpenter, A. E., Jones, T. R., Lamprecht, M. R., Clarke, C., Kang, I. H., Friman, O., et al. (2006). CellProfiler: Image analysis software for identifying and quantifying cell phenotypes. Genome Biol. 7:R100. doi: 10.1186/gb-2006-7-10-r100

PubMed Abstract | CrossRef Full Text | Google Scholar

Cartledge, J. E., Kane, C., Dias, P., Tesfom, M., Clarke, L., Mckee, B., et al. (2015). Functional crosstalk between cardiac fibroblasts and adult cardiomyocytes by soluble mediators. Cardiovasc. Res. 105, 260–270. doi: 10.1093/cvr/cvu264

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M. M., Lam, A., Abraham, J. A., Schreiner, G. F., and Joly, A. H. (2000). CTGF expression is induced by TGF- beta in cardiac fibroblasts and cardiac myocytes: a potential role in heart fibrosis. J. Mol. Cell. Cardiol. 32, 1805–1819. doi: 10.1006/jmcc.2000.1215

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W., and Frangogiannis, N. G. (2013). Fibroblasts in post-infarction inflammation and cardiac repair. Biochim. Biophys. Acta 1833, 945–953. doi: 10.1016/j.bbamcr.2012.08.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, S. A., Richardson, W. J., and Holmes, J. W. (2016). Modifying the mechanics of healing infarcts: Is better the enemy of good? J. Mol. Cell. Cardiol. 93, 115–124. doi: 10.1016/j.yjmcc.2015.11.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Czubryt, M. P. (2012). Common threads in cardiac fibrosis, infarct scar formation, and wound healing. Fibrogenesis Tissue Repair 5:19. doi: 10.1186/1755-1536-5-19

PubMed Abstract | CrossRef Full Text | Google Scholar

Dower, S. K., Kronheim, S. R., March, C. J., Conlon, P. J., Hopp, T. P., Gillis, S., et al. (1985). Detection and characterization of high affinity plasma membrane receptors for human interleukin 1. J. Exp. Med. 162, 501–515. doi: 10.1084/jem.162.2.501

PubMed Abstract | CrossRef Full Text | Google Scholar

Enríquez-de-Salamanca, A., Calder, V., Gao, J., Galatowicz, G., García-Vázquez, C., Fernández, I., et al. (2008). Cytokine responses by conjunctival epithelial cells: an in vitro model of ocular inflammation. Cytokine 44, 160–167. doi: 10.1016/j.cyto.2008.07.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Fallahi-Sichani, M., Schaller, M. A., Kirschner, D. E., Kunkel, S. L., and Linderman, J. J. (2010). Identification of key processes that control tumor necrosis factor availability in a tuberculosis granuloma. PLoS Comput. Biol. 6:e1000778. doi: 10.1371/journal.pcbi.1000778

PubMed Abstract | CrossRef Full Text | Google Scholar

Figueredo, G. P., Siebers, P. O., Owen, M. R., Reps, J., and Aickelin, U. (2014). Comparing stochastic differential equations and agent-based modelling and simulation for early-stage cancer. PLoS ONE 9:e95150. doi: 10.1371/journal.pone.0095150

PubMed Abstract | CrossRef Full Text | Google Scholar

Fix, C., Bingham, K., and Carver, W. (2011). Effects of interleukin-18 on cardiac fibroblast function and gene expression. Cytokine 53, 19–28. doi: 10.1016/j.cyto.2010.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Fomovsky, G. M., and Holmes, J. W. (2010). Evolution of scar structure, mechanics, and ventricular function after myocardial infarction in the rat. Am. J. Physiol. Heart Circ. Physiol. 298, H221–228. doi: 10.1152/ajpheart.00495.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Fomovsky, G. M., Rouillard, A. D., and Holmes, J. W. (2012). Regional mechanics determine collagen fiber structure in healing myocardial infarcts. J. Mol. Cell. Cardiol. 52, 1083–1090. doi: 10.1016/j.yjmcc.2012.02.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Frantz, S., Bauersachs, J., and Ertl, G. (2009). Post-infarct remodelling: contribution of wound healing and inflammation. Cardiovasc. Res. 81, 474–481. doi: 10.1093/cvr/cvn292

PubMed Abstract | CrossRef Full Text | Google Scholar

Fredj, S., Bescond, J., Louault, C., Delwail, A., Lecron, J. C., and Potreau, D. (2005). Role of interleukin-6 in cardiomyocyte/cardiac fibroblast interactions during myocyte hypertrophy and fibroblast proliferation. J. Cell. Physiol. 204, 428–436. doi: 10.1002/jcp.20307

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, A., and Hao, W. (2017). Mathematical modeling of liver fibrosis. Math. Biosci. Eng. 14, 143–164. doi: 10.3934/mbe.2017010

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, A., and Siewe, N. (2018). Chronic hepatitis B virus and liver fibrosis: a mathematical model. PLoS ONE 13:e0195037. doi: 10.1371/journal.pone.0195037

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerhartz, C., Dittrich, E., Stoyan, T., Rose-John, S., Yasukawa, K., Heinrich, P. C., et al. (1994). Biosynthesis and half-life of the interleukin-6 receptor and its signal transducer gp130. Eur. J. Biochem. 223, 265–274. doi: 10.1111/j.1432-1033.1994.tb18991.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Grell, M., Wajant, H., Zimmermann, G., and Scheurich, P. (1998). The type 1 receptor (CD120a) is the high-affinity receptor for soluble tumor necrosis factor. Proc. Natl. Acad. Sci. USA. 95, 570–575. doi: 10.1073/pnas.95.2.570

PubMed Abstract | CrossRef Full Text | Google Scholar

Hammerman, H., Kloner, R. A., Hale, S., Schoen, F. J., and Braunwald, E. (1983a). Dose-dependent effects of short-term methylprednisolone on myocardial infarct extent, scar formation, and ventricular function. Circulation 68, 446–452. doi: 10.1161/01.CIR.68.2.446

PubMed Abstract | CrossRef Full Text | Google Scholar

Hammerman, H., Kloner, R. A., Schoen, F. J., Brown, E. J., Hale, S., and Braunwald, E. (1983b). Indomethacin-induced scar thinning after experimental myocardial infarction. Circulation 67, 1290–1295. doi: 10.1161/01.CIR.67.6.1290

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, W., Marsh, C., and Friedman, A. (2015). A mathematical model of idiopathic pulmonary fibrosis. PLoS ONE 10:e0135097. doi: 10.1371/journal.pone.0135097

PubMed Abstract | CrossRef Full Text | Google Scholar

Hawinkels, L. J., Verspaget, H. W., van Duijn, W., van der Zon, J. M., Zuidwijk, K., Kubben, F. J., et al. (2007). Tissue level, activation and cellular localisation of TGF-beta1 and association with survival in gastric cancer patients. Br. J. Cancer 97, 398–404. doi: 10.1038/sj.bjc.6603877

PubMed Abstract | CrossRef Full Text | Google Scholar

Hazuda, D. J., Lee, J. C., and Young, P. R. (1988). The kinetics of interleukin 1 secretion from activated monocytes. Differences between interleukin 1 alpha and interleukin 1 beta. J. Biol. Chem. 263, 8473–8479.

PubMed Abstract | Google Scholar

Horiguchi, M., Ota, M., and Rifkin, D. B. (2012). Matrix control of transforming growth factor-beta function. J. Biochem. 152, 321–329. doi: 10.1093/jb/mvs089

PubMed Abstract | CrossRef Full Text | Google Scholar

Issafras, H., Corbin, J. A., Goldfine, I. D., and Roell, M. K. (2014). Detailed mechanistic analysis of gevokizumab, an allosteric anti-IL-1beta antibody with differential receptor-modulating properties. J. Pharmacol. Exp. Ther. 348, 202–215. doi: 10.1124/jpet.113.205443

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Y. F., Han, H. C., Berger, J., Dai, Q., and Lindsey, M. L. (2011). Combining experimental and mathematical modeling to reveal mechanisms of macrophage-dependent left ventricular remodeling. BMC Syst. Biol. 5:60. doi: 10.1186/1752-0509-5-60

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambert, J. M., Lopez, E. F., and Lindsey, M. L. (2008). Macrophage roles following myocardial infarction. Int. J. Cardiol. 130, 147–158. doi: 10.1016/j.ijcard.2008.04.059

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. J., Talman, L., Peirce, S. M., and Holmes, J. W. (2019). Spatial scaling in multiscale models: methods for coupling agent-based and finite-element models of wound healing. Biomech. Model Mechanobiol. 18, 1297–1309. doi: 10.1007/s10237-019-01145-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Loftis, M. J., Sexton, D., and Carver, W. (2003). Effects of collagen density on cardiac fibroblast behavior and gene expression. J. Cell. Physiol. 196, 504–511. doi: 10.1002/jcp.10330

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Y., Iyer, R. P., Jung, M., Czubryt, M. P., and Lindsey, M. L. (2017). Cardiac fibroblast activation post-myocardial infarction: current knowledge gaps. Trends Pharmacol. Sci. 38, 448–458. doi: 10.1016/j.tips.2017.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

MacKenna, D. A., Omens, J. H., McCulloch, A. D., and Covell, J. W. (1994). Contribution of collagen matrix to passive left ventricular mechanics in isolated rat hearts. Am. J. Physiol. 266, H1007–1018. doi: 10.1152/ajpheart.1994.266.3.H1007

PubMed Abstract | CrossRef Full Text | Google Scholar

Maeda, S., Dean, D. D., Gomez, R., Schwartz, Z., and Boyan, B. D. (2002). The first stage of transforming growth factor beta1 activation is release of the large latent complex from the extracellular matrix of growth plate chondrocytes by matrix vesicle stromelysin-1 (MMP-3). Calcif. Tissue Int. 70, 54–65. doi: 10.1007/s002230010032

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, K. S., Kegelman, C. D., Virgilio, K. M., Passipieri, J. A., Christ, G. J., Blemker, S. S., et al. (2017). In silico and in vivo experiments reveal M-CSF injections accelerate regeneration following muscle laceration. Ann. Biomed. Eng. 45, 747–760. doi: 10.1007/s10439-016-1707-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, K. S., Virgilio, K. M., Peirce, S. M., and Blemker, S. S. (2016). Computational modeling of muscle regeneration and adaptation to advance muscle tissue regeneration strategies. Cells Tissues Organs. 202, 250–266. doi: 10.1159/000443635

PubMed Abstract | CrossRef Full Text | Google Scholar

Mia, M. M., Boersema, M., and Bank, R. A. (2014). Interleukin-1beta attenuates myofibroblast formation and extracellular matrix production in dermal and lung fibroblasts exposed to transforming growth factor-beta1. PLoS ONE 9:e91559. doi: 10.1371/journal.pone.0091559

PubMed Abstract | CrossRef Full Text | Google Scholar

Mouton, A. J., Ma, Y., Rivera Gonzalez, O. J., Daseke, M. J., Flynn, E. R., Freeman, T. C., et al. (2019). Fibroblast polarization over the myocardial infarction time continuum shifts roles from inflammation to angiogenesis. Basic Res. Cardiol. 114:6. doi: 10.1007/s00395-019-0715-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Nim, H. T., Furtado, M. B., Costa, M. W., Kitano, H., Rosenthal, N. A., and Boyd, S. E. (2015). CARFMAP: a curated pathway map of cardiac fibroblasts. PLoS ONE 10:e0143274. doi: 10.1371/journal.pone.0143274

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez-Rodríguez, S., Tomás-González, E., and García-Aznar, J. M. (2018). 3D cell migration studies for chemotaxis on microfluidic-based chips: a comparison between cardiac and dermal fibroblasts. Bioengineering 5:45. doi: 10.3390/bioengineering5020045

PubMed Abstract | CrossRef Full Text | Google Scholar

Reed, M. J., Vernon, R. B., Abrass, I. B., and Sage, E. H. (1994). TGF-beta 1 induces the expression of type I collagen and SPARC, and enhances contraction of collagen gels, by fibroblasts from young and aged donors. J. Cell. Physiol. 158, 169–179. doi: 10.1002/jcp.1041580121

PubMed Abstract | CrossRef Full Text | Google Scholar

Rollins, B. J., O'Connell, T. M., Bennett, G., Burton, L. E., Stiles, C. D., and Rheinwald, J. G. (1989). Environment-dependent growth inhibition of human epidermal keratinocytes by recombinant human transforming growth factor-beta. J. Cell. Physiol. 139, 455–462. doi: 10.1002/jcp.1041390302

PubMed Abstract | CrossRef Full Text | Google Scholar

Rouillard, A. D., and Holmes, J. W. (2012). Mechanical regulation of fibroblast migration and collagen remodelling in healing myocardial infarcts. J. Physiol. 590, 4585–4602. doi: 10.1113/jphysiol.2012.229484

PubMed Abstract | CrossRef Full Text | Google Scholar

Spinale, F. G., Frangogiannis, N. G., Hinz, B., Holmes, J. W., Kassiri, Z., and Lindsey, M. L. (2016). Crossing into the next frontier of cardiac extracellular matrix research. Circ. Res. 119, 1040–1045. doi: 10.1161/CIRCRESAHA.116.309916

PubMed Abstract | CrossRef Full Text | Google Scholar

Stern, J. R., Christley, S., Zaborina, O., Alverdy, J. C., and An, G. (2012). Integration of TGF-beta- and EGFR-based signaling pathways using an agent-based model of epithelial restitution. Wound Repair Regen. 20, 862–871. doi: 10.1111/j.1524-475X.2012.00852.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Thannickal, V. J., Lee, D. Y., White, E. S., Cui, Z., Larios, J. M., Chacon, R., et al. (2003). Myofibroblast differentiation by transforming growth factor-beta1 is dependent on cell adhesion and integrin signaling via focal adhesion kinase. J. Biol. Chem. 278, 12384–12389. doi: 10.1074/jbc.M208544200

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, N. A. (2014). Effects of interleukin-1 on cardiac fibroblast function: relevance to post-myocardial infarction remodelling. Vascul. Pharmacol. 60, 1–7. doi: 10.1016/j.vph.2013.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, N. A., Das, A., Warburton, P., O'Regan, D. J., Ball, S. G., and Porter, K. E. (2009). Interleukin-1alpha stimulates proinflammatory cytokine expression in human cardiac myofibroblasts. Am. J. Physiol. Heart Circ. Physiol. 297, H1117–H1127. doi: 10.1152/ajpheart.00372.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, N. A., Mughal, R. S., Warburton, P., O'Regan, D. J., Ball, S. G., and Porter, K. E. (2007). Mechanism of TNFalpha-induced IL-1alpha, IL-1beta and IL-6 expression in human cardiac fibroblasts: effects of statins and thiazolidinediones. Cardiovasc. Res. 76, 81–90. doi: 10.1016/j.cardiores.2007.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

van Spreeuwel, A. C. C., Bax, N. A. M., van Nierop, B. J., Aartsma-Rus, A., Goumans, M. T. H., and Bouten, C. V. C. (2017). Mimicking cardiac fibrosis in a dish: fibroblast density rather than collagen density weakens cardiomyocyte function. J. Cardiovasc. Transl. Res. 10, 116–127. doi: 10.1007/s12265-017-9737-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Virgilio, K. M., Martin, K. S., Peirce, S. M., and Blemker, S. S. (2018). Agent-based model illustrates the role of the microenvironment in regeneration in healthy and mdx skeletal muscle. J. Appl. Physiol. (1985) 125, 1424–1439. doi: 10.1152/japplphysiol.00379.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Wakefield, L. M., Smith, D. M., Masui, T., Harris, C. C., and Sporn, M. B. (1987). Distribution and modulation of the cellular receptor for transforming growth factor-beta. J. Cell Biol. 105, 965–975. doi: 10.1083/jcb.105.2.965

PubMed Abstract | CrossRef Full Text | Google Scholar

Walpole, J., Mac Gabhann, F., Peirce, S. M., and Chappell, J. C. (2017). Agent-based computational model of retinal angiogenesis simulates microvascular network morphology as a function of pericyte coverage. Microcirculation 24:e12393. doi: 10.1111/micc.12393

PubMed Abstract | CrossRef Full Text | Google Scholar

Walpole, J., Papin, J. A., and Peirce, S. M. (2013). Multiscale computational models of complex biological systems. Annu. Rev. Biomed. Eng. 15, 137–154. doi: 10.1146/annurev-bioeng-071811-150104

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Yang, T., Ma, Y., Halade, G. V., Zhang, J., Lindsey, M. L., et al. (2012). Mathematical modeling and stability analysis of macrophage activation in left ventricular remodeling post-myocardial infarction. BMC Genomics 13(Suppl. 6):S21. doi: 10.1186/1471-2164-13-S6-S21

PubMed Abstract | CrossRef Full Text | Google Scholar

Ware, M. F., Wells, A., and Lauffenburger, D. A. (1998). Epidermal growth factor alters fibroblast migration speed and directional persistence reciprocally and in a matrix-dependent manner. J. Cell. Sci. 111(Pt 16), 2423–2432.

PubMed Abstract | Google Scholar

Warsinske, H. C., DiFazio, R. M., Linderman, J. J., Flynn, J. L., and Kirschner, D. E. (2017). Identifying mechanisms driving formation of granuloma-associated fibrosis during Mycobacterium tuberculosis infection. J. Theor. Biol. 429, 1–17. doi: 10.1016/j.jtbi.2017.06.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Warsinske, H. C., Wheaton, A. K., Kim, K. K., Linderman, J. J., Moore, B. B., and Kirschner, D. E. (2016). Computational modeling predicts simultaneous targeting of fibroblasts and epithelial cells is necessary for treatment of pulmonary fibrosis. Front. Pharmacol. 7:183. doi: 10.3389/fphar.2016.00183

PubMed Abstract | CrossRef Full Text | Google Scholar

Zahn, G., and Greischel, A. (1989). Pharmacokinetics of tumor necrosis factor alpha after intravenous administration in rats. Dose dependence and influence of tumor necrosis factor beta. Arzneimittelforschung. 39, 1180–1182.

PubMed Abstract | Google Scholar

Zeigler, A. C., Richardson, W. J., Holmes, J. W., and Saucerman, J. J. (2016a). A computational model of cardiac fibroblast signaling predicts context-dependent drivers of myofibroblast differentiation. J. Mol. Cell. Cardiol. 94, 72–81. doi: 10.1016/j.yjmcc.2016.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeigler, A. C., Richardson, W. J., Holmes, J. W., and Saucerman, J. J. (2016b). Computational modeling of cardiac fibroblasts and fibrosis. J. Mol. Cell. Cardiol. 93, 73–83. doi: 10.1016/j.yjmcc.2015.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar