Overall model

In our earlier work, we developed a SARS-CoV-2 lung tissue simulator [46], which we refer to as the “overall model”. Here, we briefly recap the most salient features of this overall model before describing our new contribution: the “fibrosis model”.

The overall model structure is developed in an open-source multiscale ABM framework PhysiCell [47]. The model simulates a tissue section of 800 μm × 800 μm × 20 μm, representing a monolayer of stationary epithelial cells on an alveolar surface of lung tissue. Initially, 2793 epithelial cells, 50 resident macrophages, 28 dendritic cells (DCs), and 57 fibroblasts are present in the simulated tissue. These initializations are directly from the code used in the overall model manuscript [46]; here, we maintain the same cellular initializations. Cells are recruited to tissue through voxels that represent vasculature points; we randomly assign 8.8% voxels as vasculature points for arrivals of recruited cells because approximately this percentage of the tissue consists of the vasculature [48]. Cellular movement and diffusion of substrates occur along a 2D plane through or above the tissue, represented by a single layer of 3D voxels; this is consistent with the transport through thin alveolar tissue. We infect the simulated tissue by placing SARS-CoV-2 viral particles in the extracellular space using a uniform random distribution. The initial number of virions is selected based on the multiplicity of infection (MOI), defined as the ratio of initial virions to the number of epithelial cells. We used a default MOI = 0.1. Viral particles diffuse through tissue, bind with unoccupied ACE2 receptors on the epithelial cell surfaces, internalize via endocytosis, replicate through intracellular viral replication kinetics, and export back to the extracellular domain by exocytosis.

PhysiCell is coupled to BioFVM [49], which solves partial differential equations for chemicals considered to be continuous (as opposed to the discrete agents) with terms for diffusion, decay, uptake, and release of multiple substances simultaneously. The extracellular viral transport, uptake by adhering to cell surface ACE2, and virion export from the cell are modeled using BioFVM formulation with Neumann boundary conditions: (1) where ρ_v is the concentration or population density of the viruses, D_v is the diffusion coefficient, U_i is the virion uptake rate by the cell, V_i is the cell volume, and E_i is the virion export rate from the cell. We used D_v = 2.5 μm²min^-1, which is within the range of influenza viral diffusion in the airway calculated in [50] (from 0.02 μm²min^-1 to 190.8 μm²min^-1). Moses et al. [51] also used a similar diffusion coefficient (1.88 μm²min^-1) for SARS-CoV-2 lung infection.

Virions bind with unoccupied external ACE2 (R_EU) using a continuum to discrete transition rule based on the binding flux: (2) where r_b is the binding rate between virion and R_EU, V_c(t) is the volume of the cell, v_e(t) is the density of the virions near a cell, and R_EU(t) is the number of R_EU at time t. The continuum to discrete transition rule allowed all nearby whole virions uptake when dv_t > 1. For fractional virion or for dv_t < 1, which resulted from the continuum rule in Eq 1, virion uptake depends on the probabilistic rule. The fractional virion uptake occurs when a random number drawn from a uniform distribution is greater than the fractional binding flux (U(0, 1) > dv_t, where 0 < dv_t < 1).

The external virus-bound ACE2 receptors endocytose, release the virions, and recycle back to the cell surface. The intracellular released virions uncoated and synthesized viral RNA and proteins and assembled for exocytosis. The receptor kinetics and intracellular viral replication kinetics are modeled as a set of ordinary differential equations (ODEs), and corresponding parameters are selected from the experimental characterization of SARS-CoV entry into host cells [52]. The assembled virions also initiate cell death response, which is modeled as Hill function to relate the apoptosis rate of a cell to assembled virions. The infected cells secrete chemokine until the cell dies through lysis or CD8+ T cell-induced apoptosis. Details of the intracellular virus model for replication kinetics, viral response, and receptor trafficking are described in much greater detail in the overall model manuscript [53].

The resident and newly recruited macrophages move along the gradients of chemokine and debris, phagocytose dead cells, and secrete pro-inflammatory cytokines. We consider five states of macrophages: inactivated, M1 phenotype, M2 phenotype, hyperactive, and exhausted. Both M1 and M2 phenotypes can also exist in hyperactive and exhausted states. Neutrophils are recruited into the tissue by pro-inflammatory cytokines, move along the gradients of chemokine and debris, phagocytose dead cells, and uptake virions. The resident DCs chemotaxis along the gradients of chemokines, activated by infected cells and viruses, and a portion of the activated DCs egress out of the tissue to the lymph node to induce activation of virus-specific CD4+ and CD8+ T cells. The activated T cells are recruited in the tissue from the lymph node. CD8+ cells attempt to attach with infected cells based on PhysiCell’s mechanical interaction and cumulative attachment above a threshold time (25 min) causes apoptosis of the infected cells. The contact between CD4+ T cells and macrophages induces a hyperactive state in macrophages, enabling them to phagocytose live infected cells. The out-of-bound cells from the simulated domain are pushed back into the tissue section by changing the direction away from the edge.

Extracellular densities of chemokines and pro-inflammatory cytokines are modeled using standard BioFVM formulation with Neumann boundary conditions: (3) where D is the diffusion coefficient, λ is the decay rate, S_i is the secretion rate, U_i is the uptake rate, and is the saturation density. For chemokines, pro-inflammatory cytokines, and debris, diffusion coefficients are assumed to be the same as the diffusion coefficient for monoclonal antibodies [54] and set at 555.56 μm²min^-1; decay rates are assumed to be the same as the decay of IL-6 [55] and set at 1.02 × 10⁻² min^-1. Further details of parameters for immune cells, equations for chemotaxis and recruitment, rules for phagocytosis and CD8+ T cell-induced apoptosis, and proliferation, activation, and clearance of the T cells in the lymph nodes using ODEs are available elsewhere [46, 53, 56]. The equations for chemotaxis and cell recruitment are also defined in the following section as they are directly used in the fibrosis model. The overall model has 5 cell types and 76 biological hypotheses, listed in the Text A in S1 Appendix, as adapted from the overall model manuscript [46]. The simulation update time for the microenvironment using BioFVM formulation is 0.01 min, cell mechanics is 0.1 min, cell processes is 6 min, and cell recruitment is 10 min. The total simulation time in the default case is 21,600 min (15 days), and the data output for saving every 60 min.

Fibrosis model

We extended the overall model to include the mechanisms for fibrosis (Fig 1). We used the ABM platform PhysiCell to simulate the rules for phenotypic and chemical changes, including cell cycle, volume changes, death processes, cellular movement, and recruitment of the agents in Fig 1, in addition to the other processes included in the overall model. The rules for the following cell types in the fibrosis model (Fig 1)—epithelial cells, infected cells, CD8+ T cells, and M1 macrophages—are all used directly from the overall model [46].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Cellular and chemical mediators in lung fibrosis during SARS-CoV-2 infection included in the fibrosis model.

Secreting agents (S) mimic the activation of latent stores of TGF-β by acting as stationary sources of TGF-β for a certain period at the site and time of death of infected epithelial cells (I). Infected cells recruit M1 macrophages initially as an innate immune response and CD8+ T cells in the later phase of infection as an adaptive immune response. M1 macrophages stop secreting pro-inflammatory cytokines when located within close proximity to CD8+ T cells and instantaneously convert to M2 macrophages. M2 macrophages secrete and activate TGF-β and are considered mobile sources. Fibroblasts become activated and recruited by TGF-β, chemotax towards the gradient of TGF-β, and deposit collagen continuously. Cells with rectangles drawn around them represent stationary cells in tissue sections, while all other cells are considered to be mobile cells. Solid arrows denote the transformation of one cell species to another, dotted lines denote interactions that influence processes without the source being produced or consumed, and the X over an arrow indicates the process does not occur. Triangles represent secreted chemical factors that can diffuse and interact in the system. Cells with different color codes represent the corresponding cells in the simulation.

https://doi.org/10.1371/journal.pcbi.1011741.g001

Several additional processes are included in the fibrosis model. These involve the secreting agents, M2 macrophages, inactive and activated fibroblasts, TGF-β, and collagen. The death of an infected epithelial cell (I in Fig 1) activates the latent TGF-β embedded in the tissue ECM. The activation of latent stores of TGF-β is represented by creating a stationary secreting agent (S in Fig 1) at the site and time of epithelial cell death. The secreting agents secrete TGF-β for a certain period as the amount of TGF-β stored in the ECM is limited; a death rate of secreting agents (AT) is used to represent terminating the activation of latent stores of TGF-β. Here, death rate AT is analogous to the extinction rate of the ECM-bound TGF-β sources and is inversely proportional to the mean duration of those sources.

Macrophages are the mobile sources of TGF-β in this fibrosis model. M1 macrophages are recruited in the initial phase, and CD8+ T cells are recruited in the later stage of infection. Usually, regulatory T cells (Tregs), a subset of CD4+ T cells, are involved in the phenotypic transformation of M1 to M2 macrophages [15]. Although we do not have Tregs in the current model, CD8+ T cells exhibit similar behavior by inhibiting the secretion of pro-inflammatory cytokines from M1 macrophages in the later phase of infection [46]. Since the fibroproliferative phase of DAD occurs in the later phase of the disease, a conditional rule is added in the fibrosis model to follow this pathology: if a CD8+ T cell and an M1 macrophage are within close proximity, then the M1 macrophage stops secreting pro-inflammatory cytokines, transforms to an M2 macrophage, and starts secreting TGF-β. The interaction distance for this transition is calculated using a multiple (ϵ) of the radii of the two interacting immune cells [46]. The transformation of the M1 to M2 phenotype is considered to occur instantaneously when the distance between an M1 macrophage and a CD8+ T cell is ≤ ϵ multiplied by the sum of the radii of the M1 macrophage and the CD8+ T cell. The cell radii are computed in PhysiCell from cell volumes over time as the cell volume changes with the progression of the cell cycle. The initial volume of macrophages is 4849 μm³ and that for CD8+ T cell is 478 μm³ with nuclear volumes of 10% of the cell volumes; these are the same values specified in the overall model [53].

We simulate extracellular concentrations of TGF-β (T_β) using the standard BioFVM formulation with Neumann boundary conditions: (4) where t is time, D_Tβ is the diffusion coefficient of TGF-β, μ_Tβ is the net decay and bulk removal rate, δ(x) is the discrete Dirac delta function, x is the center of the voxel, x_i is the position of the center of discrete cell i, and E_Tβ,i is the net export rate of TGF-β from cell i. Net export rate E_Tβ,i denotes activation rate of latent stores of TGF-β at the damaged site, which we term damaged-site secretion (DS), or TGF-β secretion and activation rate from macrophages, which we term macrophage secretion (MS).

Initial fibroblasts remain in an inactive state. The experiments of mice model of MI showed a homeostatic value of fibroblasts without external perturbation or injury [57]. To maintain the homeostatic population of fibroblasts, we assume that the initial inactive fibroblasts move randomly throughout the domain and do not undergo apoptosis (i.e., μ_F = 0). TGF-β activates the resident inactive fibroblasts. We assume that inactive fibroblasts become active fibroblasts in the presence of TGF-β (T_β > 0) and switch back to the inactive state in the absence of TGF-β (T_β = 0). The active fibroblasts and inactive fibroblasts from the active state undergo apoptosis naturally at a rate of μ_F and become dead cells. Both inactive and active fibroblasts are subject to the cycle cell model in PhysiCell where their cell volumes are over updated time as the cell cycle progresses. The initial volume of a fibroblast (F_V) and the corresponding nucleus (F_VN) are specified.

The activated fibroblasts chemotax up the gradient of TGF-β. The chemotaxis velocity of a fibroblast () depends on the TGF-β concentration in its neighboring regions, fibroblast speed (s_mot), migration bias (b), and time for persistence in a specific trajectory (Δt_mot) [46]: (5) where s_mot is the speed of chemotaxis, 0 ≤ b ≤ 1 is the level of migration bias (i.e., b = 0 represents no influence of chemotaxis and only random cell migration), is a random unit vector direction, and is the migration bias direction. Fibroblasts also persist on their given trajectory for Δt_mot before a new trajectory is computed. For fibroblasts, is estimated from the gradient of TGF-β concentrations in the neighboring regions [46]: (6) Chemotaxis of other immune cells along the gradients of pro-inflammatory cytokines, chemokines, and debris in the overall model follows this formulation [46].

New fibroblasts are recruited into the tissue by TGF-β. The number of fibroblasts recruited to the tissue (N_r) is determined by integrating the recruitment signal over the whole domain: (7) where Ω is the computational domain, r_r is the recruitment rate (per volume), s_min is the minimum recruitment signal, s_sat is the saturating or maximum signal, and s_cytokine is the recruitment signal that depends on the cytokine concentration. The volume of each grid voxel is dV = 20 × 20 × 20 μm³, and the time interval for recruitment of fibroblasts is Δt_r = 10 min, the same as in the overall model [46]. Eq 7 gives the number of fibroblasts recruited between t and t + Δt_r. Here, s_cytokine = F_g(T_β) for the recruitment signal of fibroblasts via the cytokine TGF-β. Recruitment of any immune cell by a specific cytokine in the overall model follows this formulation; proinflammatory cytokine-dependent recruitment of neutrophils and macrophages is already included in the overall model [46].

The function for TGF-β-dependent recruitment of fibroblasts (F_g(T_β))—particularly its specific parameters—comes from previously published works [40, 43]. The parameters are estimated based on experimental data and illustrated in Fig A.A in S1 Appendix. We selected a threshold value to set the upper range of TGF-β concentration beyond which F_g(T_β) is assumed to be constant because of the polynomial nature of the recruitment signal. (8)

The activated fibroblasts deposit collagen continuously. The newly deposited collagen is assumed to not diffuse or degrade. TGF-β-dependent collagen deposition from fibroblasts is described by (9) and (10) where C is the concentration of collagen, E_C,i is the net export rate of collagen from fibroblast cell i, and k_FC is the collagen production rate by fibroblasts. F_c(T_β) defines the TGF-β dependence of fibroblast-mediated collagen deposition (Fig A.B in S1 Appendix), which is assumed to follow a Michaelis-Menten form, where V_Tβ and k_Tβ are the Michaelis-Menten limiting rate and half saturation parameters, respectively.

Parameter estimation

Experimental data from the literature used in this section were extracted from graphs in the original sources (referenced in the text) using the open-source Java program Plot Digitizer [58]. The initial values of the species in the uninfected condition are listed in Table 1. The estimated and calculated parameters are listed in Table 2. Several parameters were taken from the literature without new parameter estimation: AT, ϵ, D_Tβ, μ_Tβ, μ_F, F_V, F_VN, V_F, s_mot, b, Δt_mot, and r_r. Their values and sources are also listed in Table 2. AT was also varied in one of the in silico experiments.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Initial conditions for model variables.

https://doi.org/10.1371/journal.pcbi.1011741.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Parameters for the fibrosis model.

https://doi.org/10.1371/journal.pcbi.1011741.t002

The initial number of fibroblasts (F₀) was selected based on in vivo experiments of mice models of MI [57]. We extracted the fibroblast population per surface area from the experimental control condition and multiplied it by the surface area of the computational domain (800 × 800 μm²) to evaluate the initial number of inactive resident fibroblasts in the simulated tissue. No fibroblasts are activated at t = 0. We considered the concentrations of TGF-β (T_β) and collagen deposited (C) to be changed from homeostatic values and set them to 0 as initial values (Table 1).

We evaluated the signal of TGF-β to recruit new fibroblasts (s_cytokine) in Eq 7 using TGF-β-dependent recruitment of fibroblasts rate from Eq 8 (Fig A.A in S1 Appendix), i.e., s_cytokine = F_g(T_β). The minimum and maximum (i.e., saturating) recruitment signals (s_min and s_sat, respectively) were calculated using Eq 8 for the lower (0 ng/mL) and upper (10 ng/mL) ranges of TGF-β, respectively.

The fibroblast collagen production rate (k_FC) was estimated from Hao et al. [41] (Text B in S1 Appendix). F_c(T_β) is the TGF-β dependent collagen deposition rate. We normalized published experimental data [62, 63] and fitted the data to a Michaelis-Menten relationship to estimate the parameters V_Tβ and k_Tβ for the F_c(T_β) formulation in Eq 10 (Fig A.B in S1 Appendix).

The TGF-β-dependent functions (Fig A in S1 Appendix) for recruitment rate of fibroblasts (F_g(T_β)) and collagen deposition rate from fibroblasts (F_c(T_β)) were defined for TGF-β ranges of 0–10 ng/mL from experiments [62–64] and mathematical models [40, 43]. We calibrated the fibrosis model by adjusting the free parameters DS and MS to keep TGF-β concentration within the defined range (Fig B in S1 Appendix). We iteratively changed DS while MS = 0 (turning off secretion from macrophages) and changed MS while DS = 0 (turning off secretion/activation from damaged sites) to estimate the concentration range of TGF-β (Fig B in S1 Appendix); our goal was to keep TGF-β concentration within the range of experimental observations (0–10 ng/mL) where the TGF-β-dependent functions for fibroblasts were defined. We selected ranges of DS from 5 × 10⁻¹⁰ to 2 × 10⁻⁹ ng/min and MS from 5 × 10⁻¹⁰ to 5 × 10⁻⁹ ng/min (Table 2).

Statistical analysis

The statistical significance is assessed by the Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction using the statannot package in Python [65]. p values are represented by: not significant (ns) p > 0.05; *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001.

Each in silico experiment is replicated 15 times to capture the variability. The variability arises from the initial random placement of virions, random assignment of vasculature points for immune cell entry, uptake of virion by epithelial cells, distribution of cytokines, chemokines, and debris in the microenvironment and number of recruited immune cells based on that, and migration bias in cell velocity. We used mean, 5th percentile, and 95th percentile to show the variations in replications.

Computational implementation

All the simulations are performed in a Dell Precision 3640 tower workstation: Intel Core Processor i9–10900K (10 core, 20 MB cache, base 3.7 GHz, up to 5.3GHz, and 32GB RAM) using hyperthreading, for 6 total execution threads. For a single run with 21,600 minutes of simulation, the total wall clock time was around 18 minutes. The repository for version 5.0 of the overall model is available at https://github.com/pc4covid19/pc4covid19/releases/tag/v5.0 [66]. We have provided the code used to produce the simulation results and analysis for the fibrosis model in a repository at https://github.com/ashleefv/covid19fibrosis [67].

Spatial behaviors of TGF-β sources affect collagen dynamics

For the spatial behaviors, we considered the cases of DM, D, and M (Fig 2). Figs D–F in S1 Appendix show the population spatial distributions over time for each agent type for these three cases. The stationary sources caused localized collagen deposition, whereas the distribution of collagen was more dispersed for mobile sources (Fig 3). We also observed two distinct peaks in TGF-β dynamics depending on the sources. Initially, the death of infected epithelial cells created secreting agents, which were stationary sources that activated latent TGF-β stored in the ECM at the damage sites. The number of these secreting agents fixed at damaged sites had a maximum around 5 days post-infection. The population reduced after that as the removal of the virions from the system occurred around 6 days (Fig C.A in S1 Appendix). As a result, we observed the first peak in TGF-β concentration around 5 days. The constant rate of activation and diffusion of TGF-β from stationary sources created localized fields with higher TGF-β concentrations. Active fibroblasts crowded at the sites of higher concentration of TGF-β and deposited collagen to cause localized distribution. In this model, M2 macrophages appear after the arrival of CD8+ T cells in the later phase of infection as we have encoded the rule that the interactions between M1 macrophages and CD8+ T cells transform M1 macrophages into the M2 phenotype, which acts as mobile sources of TGF-β. We observed the maximum number of macrophages around 9 days, which caused the second peak in the TGF-β profile. In the later phase of infection, damage was already present in the tissue when M2 macrophages arrived. M2 macrophages moved along the edge of damaged sites, phagocytosed dead and infected cells, cleared viral particles, and simultaneously secreted TGF-β. The movement of M2 macrophages along the edge of damaged sites caused higher concentrations of TGF-β at those sites. As a result, fibroblasts chemotaxed to those sites and deposited collagen. The deposited collagen was more dispersed with a higher collagen area fraction in the cases (DM and M) with mobile sources of TGF-β secretion compared to those (case D) with only stationary sources of TGF-β activation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Effects of stationary and mobile sources of TGF-β on the dynamics of collagen deposition.

(A–D) Case DM: equal TGF-β activation rate from stationary damaged sites and secretion rate from mobile macrophages, (E–H) Case D: TGF-β activation from stationary damaged sites only, and secretion from mobile macrophages is turned off, (I–L) Case M: TGF-β secretion from mobile macrophages only, and activation rate from stationary damaged sites is turned off. (A, E, I) Dynamics of populations of CD8+ T cells, macrophages, latent TGF-β activation sites (secreting agents), and fibroblasts; (B, F, J) dynamics of TGF-β concentration spatially averaged over the domain; (C, G, K) collagen deposition at damaged sites of tissue at day 15; (D, H, L) heat map showing collagen area fraction (AF) above the threshold value of 1 × 10⁻⁷ μg μm^-3. The solid curves represent the mean of predictions, and shaded areas represent the predictions between the 5th and 95th percentile of 15 replications of the agent-based model. The cases are defined in Fig 2.

https://doi.org/10.1371/journal.pcbi.1011741.g003

Comparison between simulated and clinical studies of collagen area fraction.

In summary for cases DM, D, and M, we observed changes in fibrotic patterns depending on the partial creation and deletion of sources, mobility of the sources, and days post infection. The combination of stationary and mobile sources of TGF-β in different phases of infection created higher TGF-β concentration zones at more sites than a single source type, resulting in higher collagen area fractions (Fig 3D, 3H and 3L). The distribution of collagen area fractions for cases DM, D, and M are shown in the violin plots in Fig 4. The collagen area fraction was lower in case D, intermediate in case M, and higher in case DM; Fig G in S1 Appendix shows the dynamics of the spatially averaged collagen concentrations for cases DM, D, and M. We observed exponential increases of spatially averaged collagen concentration in cases DM and M with a steeper slope in case DM; case D reached a plateau.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Quantification of collagen area fraction for different cases.

Violin plots represent distribution, red solid circles represent data for 15 replications of each of the cases, and box plots show medians and quartiles. Case DM: equal TGF-β production rate from both damaged sites and macrophages; case D: TGF-β activation from damaged sites only and secretion from macrophages are turned off; case M: TGF-β secretion from macrophages only and activation rate from damaged sites are turned off; case DHMH: TGF-β activation rate from damaged sites is high and TGF-β secretion rate from macrophages is high; case DLMH: TGF-β activation rate from damaged sites is low and TGF-β secretion rate from macrophages is high; case DHML: TGF-β activation rate from damaged sites is high and TGF-β secretion rate from macrophages is low; case DLML: TGF-β activation rate from damaged sites is low and TGF-β secretion rate from macrophages is low; case DA: longer activation length of latent TGF-β in case D; case MA: longer activation length of M2 macrophages in case M. The horizontal dashed lines and shading indicate the ranges of collagen extension reported in Ball et al. [34]. The cases listed along the x-axis are defined in Fig 2. p values are represented by: ns p > 0.05; *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001.

https://doi.org/10.1371/journal.pcbi.1011741.g004

Literature values for measurements performed on postmortem transbronchial cryobiopsy samples of COVID-19 pneumonia diseased tissue were reported with 2%–38% collagen extension considering all the lobes of the lung [34]. Analogously to our simulation domain, the percentages of collagen in the samples were determined by restricting the samples to the alveolar tissue while excluding airway walls, vessels, and airspaces. Our simulated collagen area fractions (Fig 4) are in similar ranges (2%–40%) compared to the experimentally observed percentage of collagen extensions at the endpoint outcomes. Our simulated collagen area fractions are also in a similar range to the reported increase of collagen volume fraction in human lung with different stages of IPF [68] and predicted collagen area fractions over time in the multiscale model and experiments of MI [45]. Although the mechanisms for fibrosis may differ in these analyses, our findings suggest that the activity of TGF-β sources is an important aspect, which alone can generate fibrotic areas similar to experimental observations.

Higher TGF-β secretion rate from M2 macrophages leads to higher collagen area fraction

We explored the effects of varying the TGF-β activation rate from stationary damaged sites and the secretion rate from mobile macrophages at different levels in combination: cases DHMH, DLML, DHML, and DLMH (Fig 2). Fig 5 shows the impacts of these variations, particularly on fibroblasts, TGF-β dynamics spatially averaged over the domain, and the TGF-β impact on the distribution of collagen deposition. Figs H–K in S1 Appendix show the population spatial distributions over time for each agent type for these four cases.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Effects of varying the TGF-β activation rate from stationary damaged sites DS and the secretion rate from mobile macrophages MS at different levels in combination.

H and L denote higher and lower values, respectively, of DS when they come after D and of MS when they come after M in the case names. (A–D) Case DHMH, (E–H) case DLML, (I–L) case DHML, and (M-P) case DLMH. (A, E, I, M) Dynamics of populations of CD8+ T cells, macrophages, latent TGF-β activation sites (secreting agents), and fibroblasts; (B, F, J, N) dynamics of TGF-β concentration spatially averaged over the domain; (C, G, K, O) dynamics of collagen concentration spatially averaged over the domain; and (D, H, L, P) heat map showing collagen area fraction (AF) above the threshold value of 1 × 10⁻⁷μg μm^-3. The solid curves represent the mean of predictions, and shaded areas represent the predictions between the 5th and 95th percentile of 15 replications of the agent-based model. The cases are defined in Fig 2.

https://doi.org/10.1371/journal.pcbi.1011741.g005

Longer duration of stationary TGF-β sources causes localized fibrotic pattern

We explored the impact of longer duration for both stationary and mobile sources on fibrotic outcomes: cases DA and MA (Fig 2).

Case DA: Longer duration of TGF-β activation from stationary sources.

Latent stores of TGF-β from damaged sites may continue to be activated for a longer duration for severe lung tissue damage. The effects of the longer duration of stationary sources (case DA) are shown in Fig 6A–6D. In this case, we turned off the secretion from mobile sources (MS = 0) to observe the response from damaged sites alone. We increased the released TGF-β activation duration 100-fold by decreasing the secreting agent death rate (AT) by a factor of 100. Fig 6A shows a greater accumulation of secreting agents over time. The concentration profile of TGF-β is proportional to the number of secreting agents. The TGF-β concentration has an approximately 2-fold larger maximum value compared to all previous cases with a peak value at an earlier phase of infection (Fig 6B). The spatial profile showed localized collagen fields with high concentrations of deposited collagen (Fig 6D). However, the violin plot of the collagen area fraction in Fig 4 showed results higher than case D and intermediate collagen area fraction compared to other cases. Therefore, even with a lower collagen fraction, localized damage with high collagen deposition can lead to long-term fibrosis by the dysregulation in the duration of activation from spatially-localized latent TGF-β sources.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Effects of longer duration of stationary and mobile sources of TGF-β.

(A–D) Case DA: longer duration of stationary damaged site sources and (E–H) case MA: longer duration of mobile macrophage sources. (A, E) Dynamics of populations of CD8+ T cells, macrophages, latent TGF-β activation sites (secreting agents), and fibroblasts; (B, F) dynamics of TGF-β concentration spatially averaged over the domain; (C, G) dynamics of collagen concentration spatially averaged over the domain; and (D, H) heat map showing collagen area fraction (AF) above the threshold value of 1 × 10⁻⁷ μg μm^-3. The solid curves represent the mean of predictions, and shaded areas represent the predictions between the 5th and 95th percentile of 15 replications of the agent-based model. The cases are defined in Fig 2.

https://doi.org/10.1371/journal.pcbi.1011741.g006

We observed localization of fibroblasts around 2 days with slow creation of stationary sources (Fig L.B and L.I in S1 Appendix), which is disrupted around 5 days by the rapid creation of new stationary sources (Fig L.C and L.J in S1 Appendix). The longer duration of stationary sources resulted in persistent localization of fibroblasts at the damaged sites from 5 days to 15 days (Fig L.C–L.G and L.I–L.N in S1 Appendix). The localization pattern remained the same during that period due to a slow decrease in stationary sources. The number of accumulated fibroblasts also increased in localized positions.

As a further in silico experiment, we changed the initial placement of virions in case DA. All cases thus far used uniform random distribution of virions. In this scenario, we considered a Gaussian distribution at the tissue center for case DA. Figs M and N in S1 Appendix showed higher fibroblast accumulation and denser collagen deposition near the center compared with the default random distribution. While only a single major fibrotic region developed (Fig M.D in S1 Appendix), the collagen concentration was the highest among all cases considered here. Thus, fibrotic patterns vary depending on the location of initial infection sites.

Our simulated results are similar to the features of focal fibrosis. Pulmonary focal fibrosis is a non-specific tissue response to injury, including infection. The feature of focal fibrosis is the sharply demarcated nodular ground-glass opacity with a maximal diameter of less than 2 cm on a thin-section computed tomography scan. The pathogenesis of focal fibrosis still needs to be explored [70]. Some studies speculated that focal fibrosis was associated with microscopic arterio-venous fistula or overexpression of platelet-derived growth factor-B [70, 71]. A multicenter observational study using 12 specimens from 11 adult COVID-19 survivors reported focal fibrosis in 50% of the patients [72]. Since we observed similar characteristic features of focal fibrosis, we hypothesized that focal fibrosis can be produced by longer activation of latent TGF-β sources.

Model validation for the dynamics of TGF-β, collagen area fraction, and macrophage cell population

Experimental data from the literature used in this section were extracted from graphs in the original sources (referenced in the text) using the web-based program WebPlotDigitizer [75]. Strobel et al. [76] developed a mouse model based on Adeno-associated-virus (AAV)-mediated expression of TGF-β1 to study progressive fibrosing interstitial lung diseases. The AAV-TGF-β1 model targeted TGF-β1 expression to the bronchial epithelium and alveolar type II cells, leading to persistent expression of TGF-β1 and histological features of lung fibrosis. The model characterized the dynamics of TGF-β1 protein levels in the bronchoalveolar lavage (BAL) at days 3, 7, 14, 21, and 28 after the start of the experiment. A higher level of TGF-β1 was observed on day 3 with a peak value on day 7, and the concentration of TGF-β1 remained high until 28 days (1.95–5.94 ng mL^-1). However, they only reported the representative tissue section showing fibrotic areas at day 21. Contrary, Sefik et al. [73] reported the percentage of affected tissue in mice models before infection and at days 2, 4, 7, 14, 28, and 35 after infection with SARS-CoV-2. However, they did not characterize the TGF-β dynamics in their experiments. Similar to the dynamics of macrophage in the experimental studies of Sefik et al. [73], Strobel et al. [76] also reported an increased number of monocytes at day 3 with a peak value at day 14, and the number of monocytes remained high until 28 days.

From our in silico experiments, we observed a higher concentration of TGF-β in the earlier phase of infection from the activation of latent TGF-β sources (case D) and persistent expression of TGF-β for longer activation of macrophages (case MA). Here, we considered a new case DHMA (DS = 2 × 10⁻⁹ ng min^-1, MS = 2 × 10⁻⁹ ng min^-1, and no apoptosis of M2 macrophages) and simulated it for 28 days to validate our simulated results with experimental observations.

Our simulated result showed similar ranges and dynamics of TGF-β concentration compared to the BAL TGF-β1 protein levels in Strobel et al. [76] (Fig 7A). In addition, we observed a peak value in the TGF-β concentration around day 5. Although TGF-β concentration at day 5 was not reported, the mathematical modeling of fibrosis during MI reported peak latent TGF-β concentration at a similar time [44], and the corresponding TGF-β mRNA level in the infracted mice model reported a peak value at day 7 [77]. The dynamics of simulated collagen area fractions are also in similar ranges compared to the percentage of affected lung tissue in Sefik et al. [73] except on day 2 and day 14 (Fig 7B). In addition, we quantified collagen area fractions using the open-source platform Fiji [78] from the day 21 representative images of NaCl control and AAV-TGF-β1 treated lung tissue section in Strobel et al. [76]. We considered the area fraction from NaCl control as an uninfected condition. Contrary to the experimental studies of Strobel et al. [76], we assumed no collagen in our simulation before infection. However, we observed a similar collagen area fraction at day 21 compared to the area fraction quantified from the AAV-TGF-β1 treated tissue section in Strobel et al. [76]. To compare the simulated M2 macrophages dynamics with the experimental data, we calculated the mean at each time point t_j for n number of samples as and scaled the cell count by dividing all the data with the maximum value of means (max(x_mean(t_j)). Our simulated dynamics of M2 macrophages were also similar to the population dynamics of alveolar macrophages reported in Sefik et al. [73] and monocyte dynamics reported in Strobel et al. [76] (Fig 7C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Model validation with experimental data of TGF-β, collagen, and macrophage cell population.

A) Comparison between experimental data of TGF-β in Strobel et al. [76] and simulated TGF-β concentration spatially averaged over the domain. B) Comparison of simulated collagen area fractions with experimental data from Sefik et al. [73] and Strobel et al. [76]. C) Comparison of the dynamics of simulated M2 macrophage count with the population dynamics of alveolar macrophages reported in Sefik et al. [73] and monocytes reported in Strobel et al. [76]. The solid curves represent the mean of predictions, and shaded areas represent the predictions between the 5th and 95th percentile of 8 replications of the agent-based model. In the box plots from Sefik et al. [73], the whiskers go down to the smallest value (minimum) and up to the largest value (maximum), the box extends from the 25th to 75th percentiles, and the median is shown as a line inside the box.

https://doi.org/10.1371/journal.pcbi.1011741.g007

Partial removal of TGF-β changes fibrotic pattern

We considered the impact of TGF-β uptake by fibroblasts as they chemotax through the domain by adding an uptake term in Eq 4: (11) where U_F,i is the uptake rate of TGF-β by fibroblast cell i and F_V,i is the volume of fibroblast cell i. The uptake rate of TGF-β by fibroblasts was assumed equivalent to the uptake rate of debris in the overall model. Here we considered cases DAU and MAU (Fig 2). Figs Q and R in S1 Appendix show the population spatial distributions over time for each agent type for these two cases. Fig S in S1 Appendix compares the collagen area fractions for these final cases to other cases.

Case DAU: Uptake of TGF-β by fibroblasts in case DA.

We considered the uptake of TGF-β in case DA and termed this scenario as case DAU. Although the population dynamics and mean collagen concentrations (Fig 8A and 8C) remained similar to those in case DA (Fig 6A and 6C), we observed variations in the spatially averaged TGF-β dynamics with a reduced peak value (Fig 8B). The uptake of TGF-β also changed the concentration gradient of TGF-β in the ECM (Eq 6), which caused delocalization and redistribution of fibroblasts in other tissue areas (Fig Q in S1 Appendix). This effectively maintained a chemotactic gradient to drive fibroblast movement. As a result, we observed a significantly higher collagen area fraction in case DAU compared to cases D and DA (Fig 8D and Fig S in S1 Appendix).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Effects of the uptake of TGF-β by fibroblasts.

(A–D) Case DAU: uptake of TGF-β by fibroblasts in case DA and (E–H) case MAU: uptake of TGF-β by fibroblasts in case MA. (A, E) Dynamics of populations of CD8+ T cells, macrophages, latent TGF-β activation sites (secreting agents), and fibroblasts; (B, F) dynamics of TGF-β concentration spatially averaged over the domain; (C, G) dynamics of collagen concentration spatially averaged over the domain; and (D, H) heat map showing collagen area fraction (AF) above the threshold value of 1 × 10⁻⁷ μg μm^-3. The solid curves represent the mean of predictions, and shaded areas represent the predictions between the 5th and 95th percentile of 15 replications of the agent-based model. The cases are defined in Fig 2.

https://doi.org/10.1371/journal.pcbi.1011741.g008