A Mathematical Model of Idiopathic Pulmonary Fibrosis

Idiopathic pulmonary fibrosis (IPF) is a disease of unknown etiology, and life expectancy of 3-5 years after diagnosis. The incidence rate in the United States is estimated as high as 15 per 100,000 persons per year. The disease is characterized by repeated injury to the alveolar epithelium, resulting in inflammation and deregulated repair, leading to scarring of the lung tissue, resulting in progressive dyspnea and hypoxemia. The disease has no cure, although new drugs are in clinical trials and two agents have been approved for use by the FDA. In the present paper we develop a mathematical model based on the interactions among cells and proteins that are involved in the progression of the disease. The model simulations are shown to be in agreement with available lung tissue data of human patients. The model can be used to explore the efficacy of potential drugs.


Introduction
Idiopathic pulmonary fibrosis (IPF) is a disease in which scar tissue in the lung is deposited; the deposition of the scar tissue is called fibrosis. As the disease progresses, alveolar-capillary units are impacted, oxygen and carbon dioxide exchange is impaired, ultimately leading to respiratory failure. IPF usually affects older people [1], but its etiology is unknown. IPF has no cure yet, and life expectancy is 3-5 years after diagnosis [2]. IPF is characterized by repeated injury to alveolar epithelium. The injury results in loss of alveolar epithelial cells (AECs) due to increased apoptosis, epithelial to mesenchymal transition (EMT), and abnormal tissue repair [3]. Oxidative stress is associated with the disregulation of the AECs [4,5], and inflammation is initiated by damaged AECs [6]. Fibrocytes, bone marrow mesenchymal progenitor cells circulating in the blood, play a role in wound repair and are increased in lungs of patients with IPF. However, fibrocyte numbers do not correlate with disease severity [7,8].
TNF-α is produced by the proinflammatory macrophages as well as by activated AEC, and it induces polarization of M2 into M1 [21] which helps to resolve the fibrosis. This polarization by TNF-α is resisted by IL-13 [22,23,24] which is produced by M2 macrophages and TH2 lymphocytes [25]. On the other hand, MMP28 [26] and other extracellular matrix (ECM) molecules (e.g. monomeric collagen type 1 interacting with CD204 on M1 [17]) activate polarization of M1 into M2 macrophages. TGF-β, along with AEC-derived basic fibroblast growth factor (bFGF) increase the proliferation of interstitial fibroblasts [6,20]. PDGF and TGF-β transform fibroblasts into myofibroblasts [27,28,29,30], which together with fibroblasts produce ECM. Imbalance between MMP and its inhibitor TIMP facilitates the accumulation of ECM and the formation of fibrosis [31].
Fibrosis is a disease in which scar tissue develops in an organ resulting in loss of functionality of the organ. Although this process evolves in nearly identical way in all organs, there may be some aspects which are organ specific. Recently Hao et. al. [32] developed a mathematical model of renal interstitial fibrosis and demonstrated that the model can be used to monitor the effect of treatment by anti-fibrotic drugs that are currently being used, or undergoing clinical trials, in non-renal fibrosis. The present paper is based on the model developed in [32] but in addition in includes two features that are unique to pulmonary fibrosis. The first one is the fact that in lung fibrosis we need to deal with two phenotypes of macrophages: monocyte-derived inflammatory macrophages (M1) and anti-inflammatory alveolar macrophages (M2). The network shown in  [32], but in the present figure the macrophages are divided into M1 and M2 phenotypes, and they play different roles in the fibrotic process.
The second unique feature in lung fibrosis is the geometry of the lung which includes a very large number of alveoli. This complex geometry is represented, in a simplified form, in The present paper develops for the first time a mathematical model of IPF. The model is based on the experimental and clinical information referenced above, schematically summarized in the network shown in Fig 1. The model is represented by a system of partial differential equations. The model is validated by comparing the simulation results with patients data and may be used to test the efficacy of potential drugs in stopping the patient's growth of fibrosis. Table 1 lists all the variables of the model in units of g/cm 3 . For the purpose of mathematical modeling we use a simple representation of the lung geometry, whose 2-dimensional projection is shown in Fig 2. The tissue under consideration is a cube R with edge-size 1 cm. The cube is partitioned by periodically arranged small cubes T ε with edge-size ε, and in each ε-cube there is a concentric cube A ε of edge-size (1 − θ)ε; the A ε represent the alveoli air space, and the domains T ε /A ε represent the alveolar tissue. An alveolar diameter is approximately 140 μm [33] and the thickness of the arterial wall which contains the capillaries, epithelial cells and fibroblasts is 10 μm. We correspondingly take 1Ày y ¼ 6, i.e., θ = 1/7. The dimensions of a lung are 12 × 31 × 41 cm 3 , and there are approximately 350 million alveoli in a lung. Hence ε is extremely small. We first write down all the differential equation in T ε /A ε , and then take ε ! 0 to obtain the homogenized system in the cube R. The variables that will be used in the model are given in Table 1.

Mathematical model
Equation for macrophage density. The equation for macrophage density in T ε /A ε (coming from the blood) is given by Macrophages are terminally differentiated cells; they do not proliferate. They differentiate from monocytes that are circulating in the blood and are attracted by MCP-1 into the lung tissue. Hence they satisfy the boundary condition wherebðPÞ depends on MCP-1 concentration, P. Here M 0 denotes the density of monoctyes in the blood, i.e., the source of M 1 macrophages from the vascular system. We note that the above Robin boundary condition arises from boundary homogenization of the vascular system, as done, for example, in [34]. The term l MT T a K T a þT a M 2 accounts for transformation from M2 to M1 induced by TNF-α [21]. The term −rÁ(M 1 χ P rP) is the chemotactic effect of MCP-1 on M1 macrophages; χ P is the chemotactic coefficient. As noted in the Introduction, macrophages from blood monocytes evolve into AM [12,16] and, in IPF, there is a shift from AM to profibrotic M2 macrophages. There is also a polarization from M1 to M2 induced by MMP28 [26], and by collagen type I via CD204 receptor on M1 [17]. The term λ M1 M1 represents polarization from M1 to M2 by the above processes and possibly other processes (e.g. [35]).
We want replace the boundary condition of M 1 , by a spatial distribution f. If D M r 2 u = f in T ε /A ε , @u @n ¼ 0 on @A ε , D M @u @n ¼ g on @T ε , then, by integration where f and g are the mean values of f and g. Since where γ = 127/343, and ε is small so that f % f and g % g, we can replace the boundary condition of M1 by the spatial distribution 6εbðPÞ=g ¼ bðPÞ. Hence, the equation for M1 density in T ε /A ε is given by with zero boundary flux. We take bðPÞ ¼ b P K P þP , where β is a constant The M2 macrophage density satisfies the equation where the first and last terms on the right-hand side are complimentary to the corresponding terms in Eq (1).
Equation for AEC density (E 0 and E). The equation of the inactivated AEC density is given by In normal healthy, the production of E 0 is represented by the term A E 0 and the death rate is represented by d E0 E 0 . The equation for the activated AEC is In homeostasis, I D = ;, δ = 0 and activated TGF-β concentration is very small. The injury to the epithelium is expressed in two ways: (i) by activation of AEC, which is represented by term where D is the damaged region and I D = 1 on D and I D = 0 elsewhere, and (ii) by increased apoptosis caused by oxidative stress [4,5] (the term δ) and by TGF-β [20,3]. In IPF, the damaged epithelium is partially repaired by fibrocytes, and this is expressed by the term [7]. The second term of the right-hand side in Eq (4) accounts for EMT due to injury [3].
Equations for fibroblast density (f) and myofibroblast density (m). The fibroblasts and myofibroblasts equations are given by: The first term on the right-hand side of Eq (5) is a source from E 0 -derived bFGF, which for simplicity we take to be in the form λ Ef E 0 . As in [32], TGF-β and PDGF transform fibroblasts into myofibroblasts [27,28,29,30]. Furthermore, TGF-β and IL-13 [22,23,24], along with Ederived bFGF, increase proliferation of fibroblasts [6,32,27]. For simplicity, we do not include bFGF specifically in the model, but instead represent it by E. The production of fibroblasts in healthy normal tissue depends on the density of AECs in homeostasis, and is represented by the term λ Ef E 0 [6,20]. Equation for ECM density (ρ) and scar (S). The ECM, produced by fibroblasts and myofibroblasts [27,28,29,30], is degraded by MMP [36], and TGF-β enhances the production of ECM by myofibroblasts [27,28,29,30]. The equation for the density of ECM is then given (as in [32]) by: where 1 À r Excessive accumulation of ECM components (particularly collagen) associated with tissue injury and inflammation, results in permanent scar formation [37]. Within each type of scar, there is considerable heterogeneity: an imbalance between MMP and TIMP activity has been implicated in the development of scar [31]. Thus a scar depends on production and deposition of ECM and disruption of normal, healthy protein cross-linking. We define the scar simply by the equation where ρ Ã is the ECM density in homeostasis and λ S is a constant, but this definition is a simplified characterization of a scar since it does not account for disruption in protein cross-linking. Equation for MCP-1 (P). The MCP-1 equation is given by where λ PE represents the growth rate by activated AEC following damage to the endothelium [32,7,14,15,1]. The last term accounts for the internalization of MCP-1 by macrophage, which may be limited due to the limited rate of receptor recycling.
As in [32], TGF-β is produced and activated by M2 macrophages while enhanced by IL-13 [22,23,24]; in addition, TGF-β is produced and activated by fibroblasts and AEC [12,20]: TNF-α is produced by M1 macrophages [21], and is also produced by AEC [12,13]: IL-13 is produced by M2 macrophages [22,23], and follows the equation Actually, IL-13 is also produced by TH2 cells [25]; for simplicity we do not include TH2 cells in our model but accounts for their production of IL-13 by the term λ I 13 .
The homogenized equations. On the boundary of T ε /A ε all the variables are assumed to have zero flux. Hence, each of the Eqs (1)- (14), if written in the form takes, after homogenization [38] (Sec. 3.1 and p.31), the following form: where γ is the volume fraction of the tissue in each ε-cube, g ¼ 127

343
. Herer 2 ¼ P a ij @ 2 @x i @x j , where the coefficient a ij are computed by where χ i satisfies the equation , n i is the i-th component of the outward normal n, and χ i is periodic in the directions of the three axes x j (j = 1,2,3). Computing a ij by finite element discretization, we find (similarly to [32]) that a ii = 0.11 (i = 1,2,3) and a ij = 0 if i 6 ¼ j.

Boundary conditions
All variables are assumed to satisfy the zero flux boundary condition on @R, the boundary of the cube R.

Initial conditions
We assume initial homeostasis, that is, λ E 0 E 0 I D = 0, but with a small amount of inflammation, represented by the term λ PE E in Eq (9). We take this term to be 10 −10 and compute the initial values by solving the steady state equations.

Numerical scheme
We briefly describe the technique used in the simulations, and for simplicity take R to be the unit cube, i.e., R = [0, 1] × [0, 1] × [0, 1]. Consider the following general diffusion equation in R with zero flux on @R. Given three positive integers K 1 , K 2 , K 3 , let Then we denote c i, j, k (t) the numerical approximation of C(x i , y j , z k , t), and get the following ODE system by semi-discretization: dc i;j;k ðtÞ @t ¼ D C ðK 2 1 ½c iþ1;j;k ðtÞ þ c iÀ1;j;k ðtÞ À 2c i;j;k ðtÞ þ K 2 2 ½c i;jþ1;k ðtÞ þ c i;jÀ1;k ðtÞ À 2c i;j;k ðtÞ þ K 2 3 ½c i;j;kÀ1 ðtÞ þ c i;j;kþ1 ðtÞ À 2c i;j;k ðtÞÞ þ F C ðc i;j;k ðtÞÞ: The Runge-kutta method is employed to solve this ODE system. The above method is used to solve the coupled system of equations of the complete model.

Treatment studies
We can use the model to explore potential drugs. Such drugs could be, for instance, anti-TGFβ, anti-PDGF, anti-IL-13 or anti-TNF-α. Anti TNF-α. To implement the effect of anti-TNF-α (TNF-α receptor that inactivates TNF-α and thus blocks TNF-α activity [44]), we need to modify the model replacing λ MT in Eqs (1) (2) by λ MT /(1 + B 1 ) to represent the inhibition of the activity of TNF-α. We assume that the drug is administered starting at day 100 from the beginning of the disease. The red curve in Figs 7 and 8 show the effect of the drug on the ECM average concentration (with B 1 = 1) over a period of 300 days. The corresponding scar has a similar curve and hence it is not given here. We see that the drug has no effect on reducing the ECM. This is in agreement with clinical phase 2 trials with Etanercept reported in [44].
The effect of the drug is introduced gradually over a period of 20 days, that is, we actually take θ(t)B 1 instead of B 1 , where θ(t) increases linearly from 0 to 1 over a period of 20 days. The same procedure is used in treatment of the subsequent drugs.
Anti-PDGF. We next consider anti-PDGF treatment, by Imatinib, an inhibitor of PDGFR and thus a blocker of PDGF activity [45]. In our model this corresponds to replacing, in Eqs (5) and (6), λ mfG by λ mfG /(1 + B 2 ). The green curve in Figs 7 and 8 show the effect of the drug on ECM for B 2 = 1. We see that the drug does not confer significant benefit, which is in agreement with phase 2 study with Imatinib.
Anti-IL-13. We next consider anti-IL-13, monoclonal antibody, a drug currently in early phase clinical trials. Tralokinamab and lebrikizumab are two drugs delivering antibody that blocks the action of IL-13. To implement their effect in our model we need to replace λ T β I 13 in Eq (13) by λ T β I 13 / (1 + B 3 ). With the choice of B 3 = 1, the blue curve in Figs 7 and 8 show no significant benefits; this seems to suggest that a moderate level of dosing will not be effective.
Anti-TGF-β. We finally consider an anti-TGFβ drug, such as Pirfenidone [46] which was recently approved in the United States. In our model we need to replace λ T β M and λ T β f by λ T β M / (1 + A) and λ T β f /(1 + A), and T β by T β /(1 + B) in all terms where T β acts to promote fibrosis. In the previous examples we showed that the drug has no benfits even at the level B = 1. For the present anti-TGF-β drug we demonstrate a clear benefit already with small A and B. Indeed, the cyan curve in Figs 7 and 8 show the effect of the drug on ECM for A = B = 0.1. We see that in terms of ECM, the drug could be effective in stopping, or even slowly decreasing fibrosis.

Discussion
IPF is a disease which exhibits, as in cutaneous wounds, both pro-inflammatory features when the alveolar epithelium is damaged and AECs begin to secrete pro-inflammatory mediators, and anti-inflammatory features associated with unsuccessful repair processes.
In this paper we developed for the first time a mathematical model for IPF. The model includes many of the principal players of cells and cytokines associated with the disease. The complex geometry of the lung alveoli is simplified by using the averaging method of homogenization, which provides a way to calculate the effective interactions among the cells and cytokines. The simulations of the model agree with lung tissue data that are available from human patients. The model can be used to explore the effect of drug treatment. Indeed, we used the model to explore the treatment of IPF by anti-TNF-α, anti-PDGF, anti-IL-13 and anti-TGF-β. We found that the first three drugs did not confer any benefits, while the last drug, pirfenidone, could be effective in stopping, or even slowly decreasing fibrosis.   We can use the model to explore novel therapeutic approaches to the treatment of IPF. For example, what will be the effect of combining two anti-fibrotic drugs? From Figs 5 and 7 we see that anti-TGF-β is the most effective drug to slow the IPF progression (with A = B = 0.1) and anti-IL-13 has only very mild benefits (with B 3 = 1). However if we combine these two drugs (at the same respective levels) we obtain significant improvement of over anti-TGF-β alone, especially in the case of severe case of IPF, as seen in the bottom curves in Figs 5 and 7. We propose this result as an hypothesis that could be checked in clinical trials.