An agent-based modeling approach for lung fibrosis in response to COVID-19

The severity of the COVID-19 pandemic has created an emerging need to investigate the long-term effects of infection on patients. Many individuals are at risk of suffering pulmonary fibrosis due to the pathogenesis of lung injury and impairment in the healing mechanism. Fibroblasts are the central mediators of extracellular matrix (ECM) deposition during tissue regeneration, regulated by anti-inflammatory cytokines including transforming growth factor beta (TGF-β). The TGF-β-dependent accumulation of fibroblasts at the damaged site and excess fibrillar collagen deposition lead to fibrosis. We developed an open-source, multiscale tissue simulator to investigate the role of TGF-β sources in the progression of lung fibrosis after SARS-CoV-2 exposure, intracellular viral replication, infection of epithelial cells, and host immune response. Using the model, we predicted the dynamics of fibroblasts, TGF-β, and collagen deposition for 15 days post-infection in virtual lung tissue. Our results showed variation in collagen area fractions between 2% and 40% depending on the spatial behavior of the sources (stationary or mobile), the rate of activation of TGF-β, and the duration of TGF-β sources. We identified M2 macrophages as primary contributors to higher collagen area fraction. Our simulation results also predicted fibrotic outcomes even with lower collagen area fraction when spatially-localized latent TGF-β sources were active for longer times. We validated our model by comparing simulated dynamics for TGF-β, collagen area fraction, and macrophage cell population with independent experimental data from mouse models. Our results showed that partial removal of TGF-β sources changed the fibrotic patterns; in the presence of persistent TGF-β sources, partial removal of TGF-β from the ECM significantly increased collagen area fraction due to maintenance of chemotactic gradients driving fibroblast movement. The computational findings are consistent with independent experimental and clinical observations of collagen area fractions and cell population dynamics not used in developing the model. These critical insights into the activity of TGF-β sources may find applications in the current clinical trials targeting TGF-β for the resolution of lung fibrosis.


TGF-β-dependent functions for fibroblasts
The function for TGF-β-dependent recruitment of fibroblasts from [2] was fit to experimental data of [3], represented by solid red circles, and is shown in   [6].We extracted data between day 4 and day 9 of the fibroblast population and the ECM concentration because we observed fibroblast activation and recruitment in our model during that period.Also, in Hao et al. [6], the population dynamics and ECM deposition remained linear during that period.We assumed fibroblast cell weight was 2 × 10 −9 g to calculate the number of fibroblasts.We used Eq S1 to estimate fibroblast collagen production rate (k F C ):

1 . 6 .Neutrophils 1 ./15 4 . 1 . 1 .Fibroblasts 1 . 4 . 5 .
Live epithelial cells undergo apoptosis after sufficient cumulative contact time with adhered CD8+ T cells 2. Dead epithelial cells produce debris 3. Virus adheres to unbound external ACE2 receptor to become external (virus)-bound ACE2 receptor 4. Bound external ACE2 receptor is internalized (endocytosed) to become internal bound ACE2 receptor 5. Internalized bound ACE2 receptor releases its virion and becomes unbound internalized receptor.The released virus is available for use by the viral lifecycle model Internalized unbound ACE2 receptor is returned to the cell surface to become external unbound receptor 7. Each receptor can bind to at most one virus particle 8. Internalized virus is uncoated 9. Uncoated virus (viral contents) lead to release of functioning RNA 10.RNA creates viral protein at a constant rate unless it degrades 11.Viral RNA is replicated at a rate that saturates with the amount of viral RNA 12. Viral RNA undergoes constitutive (first order) degradation December 14, 2023 1/15 13.Viral protein is transformed to an assembled virus state 14.Assembled virus is released by the cell (exocytosis) 15.After infection, cells secrete chemokine 16.As a proxy for viral disruption of the cell, the probability of cell death increases with the total number of assembled virions 17.Apoptosed cells lyse and release some or all their contents 18.Once viral RNA exceeds a particular threshold, the cell enters the pyroptosis cascade 19.Once pyropotosis begins, the intracellular cascade is modeled by a system of ODEs monitoring cytokine production and cell volume swelling 20.Cell secretion rate for pro-inflammatory increases to include secretion rate of IL-18 21.Cell secretes IL-1β which causes a bystander effect initiating pyroptosis in neighboring cells 22. Cell lyses (dying and releasing its contents) once its volume has exceeded 1.5× the homeostatic volume 23.Infected epithelial cells secrete pro-inflammatory cytokine 24.Antigen presentation in infected cells is a function of intracellular viral proteinMacrophages 1. Resident (unactivated) and newly recruited macrophages move along debris gradients 2. Macrophages phagocytose dead cells.Time taken for material phagocytosis is proportional to the size of the debris 3. Macrophages break down phagocytosed materials 4.After phagocytosing dead cells, macrophages activate and secrete pro-inflammatory cytokines 5. Activated macrophages can decrease migration speed 6. Activated macrophages have a higher apoptosis rate 7. Activated macrophages migrate along chemokine and debris gradients 8. Macrophages are recruited into tissue by pro-inflammatory cytokines 9. Macrophages can die and become dead cells only if they are in an exhausted state 10.Macrophages become exhausted (stop phagocytosing) if internalized debris is above a threshold 11.CD4+ T cell contact induces activated macrophage phagocytosis of live infected cells Neutrophils are recruited into the tissue by pro-inflammatory cytokines 2. Neutrophils die naturally and become dead cells 3. Neutrophils migrate locally in the tissue along chemokine and debris gradients December 14, 2023 2Neutrophils phagocytose dead cells and activate 5. Neutrophils break down phagocytosed materials 6. Activated neutrophils reduce migration speed 7. Neutrophils uptake virus 8. Neutrophils secrete ROS upon phagocytosis Dendritic cells (DCs) 1. Resident DCs exist in the tissue 2. DCs are activated by infected cells and/or virus 3. Portion of activated DCs leave the tissue to travel to the lymph node 4. DCs chemotaxis up chemokine gradient 5. Activated DCs present antigen to CD8+ T cells increasing their proliferation rate and killing efficacy (doubled proliferation rate and attachment rate) 6. Activated DCs also regulate the CD8+ T cell levels in within a threshold by enhancing CD8+ T cell clearance CD8+ T cells 1. CD8+ T cells are recruited into the tissue by pro-inflammatory cytokines 2. CD8+ T cells apoptose naturally and become dead cells 3. CD8+ T cells move locally in the tissue along chemokine gradients 4. CD8+ T cells adhere to infected cells.Cumulated contact time with adhered CD8+ T cells can induce apoptosis 5. Activated DCs present antigen to CD8+ T cells, which increases the CD8+ T cell proliferation rate 6.Activated DCs also regulate the CD8+ T cell levels in within a threshold by enhancing CD8+ T cell clearance 7. CD8+ T cells have a max generation counter and will not proliferate after the set generation CD4+ T cells 1. CD4+ T cells are recruited into the tissue by the lymph node 2. CD4+ T cells apoptose naturally and become dead cells 3. CD4+ T cells move locally in the tissue along chemokine gradients 4. CD4+ T cells are activated in the lymph node by three signals: antigenic presentation by the DCs, direct activation by cytokines secreted by DCs, and direct activation by cytokines secreted by CD4+ T cells 5. CD4+ T cells are suppressed directly by cytokines secreted by CD4+ T cells 6. CD4+ T cells have a max generation counter and will not proliferate after the set generation Tissue microenvironment Virus diffuses in the microenvironment 2. Virus adhesion to a cell stops its diffusion (acts as an uptake term) 3. Pro-inflammatory cytokine diffuses in the microenvironment December 14, 2023 3/15 4. Pro-inflammatory cytokine is taken up by recruited immune cells 5. Pro-inflammatory cytokine is eliminated or cleared 6. Chemokine diffuses in the microenvironment 7. Chemokine is taken up by immune cells during chemotaxis 8. Chemokine is eliminated or cleared 9. Debris diffuses in the microenvironment 10.Debris is taken up by macrophages and neutrophils during chemotaxis 11.Debris is eliminated or cleared 12. Immunoglobulin (Ig) diffuses in the microenvironment 13.Virions and Ig react in tissue removing both components 14.Reactive oxidative species (ROS) diffuses in the microenvironment Fibrosis model We extended the overall model to include the mechanisms of fibrosis.The details of the cell types and biological hypotheses are described in detail in the manuscript.Here, we summarize the cell types and biological hypotheses (agent-based model rules) for the fibrosis model.M2 Macrophages 1. CD8+ T cell contact stops activated macrophage secretion of pro-inflammatory cytokine and switches to M2 phase 2. M2 macrophages secrete TGF-β Secreting agent Secreting agents are created at the site that a CD8+ T cell kills an epithelial cell to mimic latent TGF-β activation embedded in the tissue 2. TGF-β cytokine is eliminated or cleared 3. Secreting agents secrete TGF-β for a set amount of time Resident inactive fibroblasts exist in the tissue 2. Resident inactive fibroblasts do not apoptose to maintain homeostasis 3. TGF-β activates the resident inactive fibroblasts and absence of TGF-β inactivates fibroblasts Active fibroblasts are recruited into the tissue by TGF-β Active fibroblasts and inactive fibroblasts from the active state apoptose naturally and become dead cells 6. Fibroblasts move locally in the tissue along up gradients of TGF-β 7. Fibroblast cells deposit collagen continuously Tissue microenvironment 1. TGF-β diffuses in the microenvironment 2. TGF-β is degraded or removed Text B. Methods Fig A.A. This fitted curve is used as F g (T β ) for recruitment of new fibroblasts in Eq 5 of the fibrosis model.We assumed constant response (green line) when TGF-β concentration exceeds the threshold value (10 ng/mL, marked by the vertical black line in Fig A.A) because of the polynomial nature of the recruitment signal.Fig A.B shows the function used for TGF-β-dependent collagen deposition rate from fibroblasts, F c (T β ).The experimental data from [4, 5], represented by solid red circles, were normalized and fitted with a Michaelis-Menten kinetic function represented by a blue curve to estimate V T β and k T β for F c (T β ) in Eq 7 of the fibrosis model.

Fig
Fig A. TGF-β concentration dependency for functions used in fibrosis model to describe (A) recruitment of fibroblasts (Eq 5) and (B) collagen deposition from fibroblasts (Eq 7).In A) the green curve is used as the F g (T β ) constant value for T β > 10.

Fig D.
Fig D. Case DM: equal TGF-β activation rate from stationary damaged sites and secretion rate from mobile macrophages.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.

Fig E.
Fig E. Case D: TGF-β activation from stationary damaged sites only, and secretion from mobile macrophages is turned off.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.

Fig F.
Fig F. Case M: TGF-β secretion from mobile macrophages only, and activation rate from stationary damaged sites is turned off.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.

Fig G.
Fig G. Effects of sources on collagen deposition.Dynamics of collagen concentration spatially averaged over the domain varying with the sources of TGF-β from (A) case DM, (B) case D, and (C) case M. 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.

Fig H.
Fig H. Case DHMH: high TGF-β activation rate from stationary damaged sites and high TGF-β secretion rate from mobile macrophages.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.

Fig I.
Fig I. Case DLML: low TGF-β activation rate from stationary damaged sites and low TGF-β secretion rate from mobile macrophages.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.

Fig J.
Fig J. Case DHML: high TGF-β activation rate from stationary damaged sites and low TGF-β secretion rate from mobile macrophages.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.

Fig K.
Fig K. Case DLMH: low TGF-β activation rate from stationary damaged site and high TGF-β secretion rate from mobile macrophages.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.

Fig M.
Fig M. Case DA with Gaussian distribution for initial placement of virions.(A) Cell populations; (B) TGF-β concentration spatially averaged over the domain; (C) collagen concentration spatially averaged over the domain; and (D) heat map showing collagen area fraction (AF ) above the threshold value of 1 × 10 −7 µ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.

Fig N.
Fig N. Case DA with Gaussian distribution for initial placement of virions.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.

Fig O.
Fig O. Case MA: extension of case M with longer duration of TGF-β secretion from mobile macrophages.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.

Fig P.
Fig P. Dynamics of populations of macrophages in different states.The overall model has five states: inactivated state (MI), M1 phenotype (M1), M2 phenotype (M2), hyperactive state (MH), and exhausted state (ME).Both M1 and M2 phenotypes can also exist in MH and ME states.Details are available in Getz et al. [1].(A) Case D, (B) case M, (C) case DA, and (D) case MA.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.

Fig R.
Fig R. Case MAU: extension of case MA with the uptake of TGF-β by fibroblasts.Simulated dynamics of cell population with (A-G) collagen and (H-N) TGF-β concentration fields, shown behind the cells (circles) and corresponding to the color bars.Representative model results (one replicate) at days (A, H) 0, (B, I) 2, (C, J) 5, (D, K) 6, (E, L) 8, (F, M) 9, and (G, N) 15.Epithelial cells are blue, macrophages are green, CD8+ T cells are red, secreting agents are black, and fibroblasts are purple.The color codes for other immune cells are the same as in Getz et al. [1].The cases are defined in Fig 2.