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].

Inflammation and injury activate AECs [9, 10, 11], and activated AECs secrete a number of pro-inflammatory mediators including tumor necrotic factor alpha (TNF-α) [12, 13] and chemoattractant monocyte chemotactic protein-1 (MCP-1) [7, 14, 15]. MCP-1 recruits circulating monocytes from the blood into damaged lung tissue, where they differentiate into classically activated macrophages M1.In normal lung tissue (homeostasis), macrophages from blood monocytes develop into alveolar macrophages (AM) [12, 16]. Alveolar macrophages are often referred to as alternatively activated macrophages, or M2 macrophages. However M2 macrophages are heterogeneous, and in IPF there appears to be a shift from monocyte-derived M1 macrophages to pro-fibrotic M2 macrophages [17, 18]. These M2 macrophages are responsible for the progression from inflammation to interstitial fibrosis [2, 18] by secreting platelet-derived growth factor (PDGF) [19, 17], transforming growth factor-beta (TGF-β) [17], matrix metalloproteinase (MMP) and tissue inhibitor of metalloproteinase (TIMP) [17], all of which are involved in the regulation of tissue fibrosis. TGF-β is produced also by fibroblasts activated by AEC [12, 20]. Both TGF-β and reactive oxygen species increase AEC apoptosis [20].

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 Fig 1 is similar to the network in Fig 1 of [32], but in the present figure the macrophages are divided into M1 and M2 phenotypes, and they play different roles in the fibrotic process.

Fig 1

Schematic network of cells and proteins in IPF.

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 Fig 2. Our mathematical model of IPF is based on Fig 1 combined with ‘homogenization’ method associated with Fig 2.

Fig 2

Lung geometry consists of a periodically arranged cubes with smaller cubes representing the air space of alveoli.

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.

Materials and Methods

Mathematical model

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 , 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.

Table 1

The variables of the model in units of g/cm 3.

M 1:	density of M1 macrophages	M 2:	density of M2 macrophages
E 0:	density of AEC	E	density of activated AEC
f:	density of fibroblast	m:	density of myofibroblast
ρ:	density of ECM	P:	concentration of MCP-1
G:	concentration of PDGF	T β:	concentration of activated TGF-β
Q:	concentration of MMPs	Q r:	concentration of TIMP
T α	concentration of TNF-α	I 13	concentration of IL-13
S	scar density

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.

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 where 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 accounts for transformation from M2 to M1 induced by TNF-α [21]. The term −∇⋅(M 1 χ ∇P) is the chemotactic effect of MCP-1 on M1 macrophages; χ 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 pro-fibrotic 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 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 ∇2 u = f in T /A , on ∂A , on ∂T , then, by integration ∫ fdV = ∫∂ gdS. Hence where and are the mean values of f and g. Since where γ = 127/343, and ɛ is small so that , we can replace the boundary condition of M1 by the spatial distribution . Hence, the equation for M1 density in T /A is given by with zero boundary flux. We take , 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 and the death rate is represented by d E 0.

The equation for the activated AEC is In homeostasis, I = ∅, δ = 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 λ E 0 I , where D is the damaged region and I = 1 on D and I = 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 λ 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 E-derived 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 λ 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 if ρ < ρ 0, if ρ ≥ ρ 0.

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 λ 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 λ 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.

Equations for concentrations of PDGF (G), MMP (Q), TIMP (Q ), TGF-β (T ), TNF-α (T ) and IL-13 (I 13)

As in [32], the following sets of diffusion equations hold for G, Q and Q : Note that in Eq (11), MMP is lost by binding with TIMP (second term).

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 λ .

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, . Here , where the coefficient a are computed by where χ satisfies the equation here , , n is the i-th component of the outward normal n, and χ is periodic in the directions of the three axes x (j = 1,2,3). Computing a by finite element discretization, we find (similarly to [32]) that a = 0.11 (i = 1,2,3) and a = 0 if i ≠ 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 I = 0, but with a small amount of inflammation, represented by the term λ E in Eq (9). We take this term to be 10−10 and compute the initial values by solving the steady state equations.

In particular we find the initial values of T = 2.5 × 10−8, T = 2.51 × 10−12 and I 13 = 3.2 × 10−8 in units of gm/ml. Taking into account that only γ-fraction of the space is occupied by tissue, the values , coincide with the concentration of T and I 13 measured in the bronchial tubes of healthy lung in [39], and coincides with value of TGF-β as computed in [40].

We also compute that E 0 = E* = 0.79 g/cm 3, f = f* = 4.75 × 10−3 g/cm 3, ρ = ρ* = 3.26 × 10−3 g/cm 3 and I 13 = 1.76 × 10−8 g/cm 3 at t = 0.

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 (t) the numerical approximation of C(x , y , z , t), and get the following ODE system by semi-discretization: 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.

Model simulation and validation

In this section, we simulate the model (1)-(17). The parameter values are listed in Tables 2 and 3 and the initial values are taken as explained above. The numerical simulation were carried out by finite difference scheme in spatial direction and Runge-Kutta method in time direction.

Table 2

Parameters’ description and value.

Parameter	Description	Value
D M	dispersion coefficient of macrophages	8.64 × 10−7 cm 2 day−1 [32]
D P	diffusion coefficient of MCP-1	1.728 × 10−1 cm 2 day−1 [32]
D G	diffusion coefficient of PDGF	8.64 × 10−2 cm 2 day−1 [32]
D Q	diffusion coefficient of MMP	4.32 × 10−2 cm 2 day−1 [32]
D Qr	diffusion coefficient for TIMP	4.32 × 10−2 cm 2 day−1 [32]
D Tβ	diffusion coefficient for TGF-β	4.32 × 10−2 cm 2 day−1 [32]
D Tα	diffusion coefficient for TNF-α	1.29 × 10−2 cm 2 day−1 [55]
D f	dispersion coefficient of fibroblasts	1.47 × 10−6 cm 2 day−1 [32]
D m	dispersion coefficient of myofibroblasts	1.47 × 10−5 cm 2 day−1 [32]
D I13	diffusion coefficient of IL-13	1.08 × 10−2 cm 2 day−1 [40]
D Tα	diffusion coefficient for TNF-α	1.29 × 10−2 cm 2 day−1 [40]
λ MT	transition rate of M2 to M1 macrophages by TNF-α	5 × 10−3 day−1 [56]
λ M1	polarization rate of M1 to M2 macrophages	9.02 × 10−6 day−1, estimated
λ E0	production rate of AEC	0.25 day−1 estimated
λ 1	repair rate of AEC	10−3 g/cm 3 day−1 estimated
λ EM	EMT rate of AEC	1.65 × 10−3 day−1 estimated
λ TβM	production rate of TGF-β by macrophages	1.5 × 10−2 day−1 [32]
λ Tβf	production rate of TGF-β by fibroblast	7.5 × 10−3 day−1 [32] & estimated
λ GM	production rate of PDGF by macrophages	2.4 × 10−5 day−1 [32]
λ QM	production rate of MMP by macrophages	3 × 10−4 day−1 [32]
λ QrM	production rate of TIMP by macrophages	6 × 10−5 day−1 [32]
λ PE	activation rate of MCP-1 due to AECs	1 × 10−8 day−1 [32]
λ ρf	activation rate of ECM due to fibroblasts	3 × 10−3 day−1 [32]
λ ρm	activation rate of ECM due to myofibroblasts	6 × 10−3 day−1 [32]
λ ρTβ	activation rate of ECM due to TGF-β	2 [32]
λ Ef	activation rate of fibroblasts due to bFGF and TGF-β	2.5 × 10−1 day−1 [32] & estimated
λ fE	production rate of fibroblasts	5 × 10−4 day−1 [32] & estimated
λ mfT	activation rate of myofibroblasts due to TGF-β	0.12 day−1 [32]
λ mfG	activation rate of myofibroblasts due to PDGF	0.12 day−1 [32]
λ TαM	activation rate of TNF-α due tomacrophage	1.39 × 10−5 day−1 [57]
λ TαE	activation rate of TNF-α due tomacrophage	6.9 × 10−6 day−1 [57] & estimated
λ I13	production rate of IL-13 by Th2 cells	2.12 × 10−7 g/ml day−1 [40]
λ I13M	production rate of IL-13 by macrophages	3.98 × 10−4 day−1 [40]
d M2	death rate of macrophages	0.015 day−1 [32]
d M1	death rate of macrophages	0.02 day−1 [56, 58]
d E	death rate of AECs	1.65 × 10−2 day−1 [32]
d E0	death rate of AECs	1.65 × 10−2 day−1 [32]
d E0T	death rate of AECs	1.65 × 10−3 day−1 [32]
δ	increased death rate of AECs by injury	1 × 10−3 day−1, estimated
d ρ	degradation rate of ECM	0.37 day−1 [32]
d P	degradation rate of MCP-1	1.73 day−1 [32]
d PM	internalization rate of MCP-1 by M1 macrophages	2.08 × 10−4 day−1 [32]
d G	degradation rate of PDGF	3.84 day−1 [32]
d QQr	binding rate of MMP to TIMP	4.98 × 108 cm 3 g −1 day−1 [32]
d QrQ	binding rate of TIMP to MMP	1.04 × 109 cm 3 g −1 day−1 [32]
d Q	degradation rate of MMP	4.32 day−1 [32]
d Qr	degradation rate of TIMP	21.6 day−1 [32]
d ρQ	degradation rate of ECM due to MMP	2.59 × 107 cm 3 g −1 day−1 [32]
d Tβ	degradation rate of TGF-β	3.33 × 102 day−1 [32]
d f	death rate of fibroblasts	1.66 × 10−2 day−1 [32]
d m	death rate of myofibroblasts	1.66 × 10−2 day−1 [32]
d Tα	degradation rate of TNF-α	55.45 day−1 [59]
d I13	degradation rate of IL-13	12.47 day−1 [40]

Table 3

Parameters’ description and value.

Parameter	Description	Value
χ P	chemotactic sensitivity parameter by MCP-1	10 cm 5 g −1 day−1 [32]
A E0	intrinsic AEC proliferation	8.27 × 10−3 g/cm 3 day−1 [32]
K G	PDGF saturation for activation of myofibroblasts	1.5 × 10−8 gcm −3 [32]
K Tβ	TGF-β saturation for apoptosis of AECs	1 × 10−10 gcm −3 [32]
K P	MCP-1 saturation for influx of macrophages	5 × 10−9 gcm −3 [32]
K Tα	TNF-α saturation	5 × 10−7 gcm −3 [40]
K I13	IL-13 saturation	2 × 10−7 g/cm 3 [40]
K E	AEC saturation	0.1 g/cm 3, estimated
ρ 0	ECM saturation	10−2 gcm −3 [32]
ρ*	ECM density in health	3.26 × 10−3 gcm −3 estimated
E*	TEC density in health	0.799 gcm −3 estimated
f*	fibroblast density in health	4.75 × 10−3 gcm −3 estimated
M 0	source/influx of macrophages from blood	5 × 10−5 gcm −3 [32]
β	influx rate of macrophages into interstitium	0.2 cm −1 [32]

Fig 3 shows the dynamics of the average densities of cells and concentrations of cytokines for 30 days.

Fig 3

The dynamics of the average concentrations of cells and cytokines in units of gm/cm 3 from homeostasis at day 0 to day 30.

I = 0.3 × 0.3 × 0.3 cm 3.

Fig 4 shows histogram of cells and cytokines in disease vs. homeostasis. The simulation results for MMP and TIMP shown in Fig 4 are in agreement with the experimental results, reported in [41] for protein concentration human lung tissue with IPF (n = 16 human subjects) and control (n = 6 human subjects). Indeed, although (in [41]) MMP 7 (for IPF) is nearly 4 times the level of MMP7 for control, all other MMPs (1,2,9,13) increased approximately twice or just a little more than twice, while the relatively small concentration of MMP8 decreased to 25% of the control level. The simulation results for TIMP shows an increase of 20% in the protein concentration for IPF vs. control, which is the same as in human lung tissues reported in [41] for TIMP-1,2,3. Levels of mRNA expression relate to levels of the translated proteins. The mRNA of TGF-β reported in [42] (which can also be deduced from [43]) shows increase by at least twice in IPF vs. control; this increase is the same for the TGF-β protein shown in Fig 4. However, we cannot make too much out of this comparison since TGF-β has to be activated post transcriptionally to be biologically active [20]. The mRNA expressions of TNF-α and PDGF reported in [43] show increased levels in IPF patients, which is in qualitative agreement with the increase in protein levels shown in Fig 4.

Fig 4

Comparison of cells and cytokines for IPF and healthy control at day 30 from the beginning of the disease (in fraction of healthy control).

Figs 5 and 6 are simulations of the disease for a larger period of 300 days. We see that the disease continue to grow but at slower rate.

Fig 5

The dynamics of the average concentrations of cells and cytokines in units of gm/cm 3 from homeostasis at day 0 to day 300.

I = 0.3 × 0.3 × 0.3 cm 3.

Fig 6

Comparison of cells and cytokines for IPF and healthy control at day 300 from the beginning of the disease (in fraction of healthy control).

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-α. Fig 7 displays the effect of treatment for mild case of IPF, namely I = 0.3 × 0.3 and λ = 2.5 × 10−3 day−1, and Fig 8 displays the effect of treatment for severe case of IPF, namely, I = 0.5 × 0.5 and λ = 3 × 10−3 day−1

Fig 7

Treatment studies for the mild case.

ECM is in units of gm/cm 3.

Fig 8

Treatment studies for the severe case.

ECM is in units of gm/cm 3.

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 λ in Eqs (1) (2) by λ /(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), λ by λ /(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 λ in Eq (13) by λ /(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 λ and λ by λ /(1 + A) and λ /(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.

The present model should be viewed as a first step in the development a more comprehensive study of IPF. Such a study should include altered DNA methylation [47, 48], epigenetic and environmental factors [49], gene mutation (e.g. of surfactant protein [50]), polymorphism (e.g. of IL-10 [51], IL-4 [52], Muc5B [53]), and telomerase mutations [54].