1. Introduction

Fibrosis is the formation of excessive fibrous connective tissue in an organ or tissue, which occurs in reparative process or in response to inflammation. The excessive deposition disorganizes the architecture of the organ or tissue, and results in scars that disrupt the function of the organ or tissue. Fibrotic diseases are characterized by abnormal excessive deposition of fibrous proteins, such as collagen, and the disease is most commonly progressive, leading to organ disfunction and failure. Fibrotic diseases may become fatal when they develop in vital organs such as heart, lung, liver and kidney. Systemic sclerosis, an autoimmune disease, is another type of fibrotic disease whereby an area that usually has flexibility and movement thickens and hardens; examples include skin, blood vessels, and muscles which tighten under fibrous deposition. Myocardial infarction, commonly known as the heart attack, occurs when blood flow stops to a part of the heart, causing damage to the heart muscles. In this case a reparative response may initiate cardiac fibrosis. On the other hand, an autoimmune disease, such as Lupus Nephritis, begins with inflammation in the kidney, which may then lead to renal tubulointerstitial fibrosis. In general, no matter what initiates the disease, a reparative process or response to inflammation, fibrosis in an organ eventually involves both reparative process and response to inflammation.

The reparative process in fibrosis is similar to the process of wound healing. Anti-inflammatory macrophages (M2) secrete TGF-$\beta$ and PDGF which activate fibroblasts and myofibroblasts to excessively produce collagen. At the same time the M2 macrophages secrete MMP (and its antagonist TIMP) which disrupts the cross-linking in the collagen network and initiates the formation of a scar. On the other hand, in response to inflammation that develops in an organ, monocytes from the blood are induced to immigrate into the organ and differentiate into M1 macrophages. A polarization from M1 to M2 then takes place to promote tissue reparation. The heterogeneity of macrophages and the exchange of M1/M2 polarization have been reported in kidney fibrosis (Ricardo et al. [68], Duffield [18]), in liver fibrosis (Tacke et al. [78], Pellicoro et al. [64]), in cardiac fibrosis (Kong et al. [41]) and in idiopathic pulmonary fibrosis (IPF) (Hao et al. [31] and the references therein).

Proinflammatory macrophages M1 produce IL-12 which activates CD4+ T cells. The activities and heterogeneity of T cells have been reported in kidney fibrosis (Tapmeier et al. [81], Liu et al. [48], Nikolic-Paterson [63]), in liver fibrosis (Connolly et al. [13], Barron et al. [5], Hammerich et al. [27], Liedtke et al. [46]), in cardiac fibrosis (Wei et al. [89]), in pulmonary fibrosis (Shimizu et al. [74], Luzina et al. [54], Kikuchi et al. [39], Wei et al. [89], Lo Re et al. [51]), and in systemic sclerosis (Chizzolini [11]). The opposing effects of Th1/Th2 in fibrotic diseases was considered already in earlier work by Wynn [91].

A key protein in enhancing production of collagen by fibroblasts/myofibroblasts is TGF-$\beta$. The activity of TGF-$\beta$ ligand is mediated by a class of SMAD proteins which form complexes that enter into the nucleus of fibroblast/myofibroblast as transcription factors (Leask et al. [44], Kahn et al. [38], Rosenbloom et al. [69]). TGF-$\beta$ represents an attractive therapeutic target in the treatment of fibrotic diseases [69].

In this paper we focus on liver fibrosis. The liver is an organ that supports the body in many ways, as illustrated in Fig. 1. It produces bile, a substance needed to digest fats which, in particular, helps synthesize cholesterol. It stores sugar glucose and converts it to functional (glucose) sugar when the body sugar levels fall below normal. The liver detoxifies the blood; modifies ammonia into urea for excretion, and destroys old red blood cells. Hepatocytes are the cells of the main parenchymal tissue of the liver, making up 60% of the total liver cells; they perform most of the tasks of the liver. Kupffer cells are stellar macrophages, in direct contact with the blood, making up 15-20% of the liver cells. Hepatocyte stellate cells (HSCs) are quiescent cells, different from the portal fibroblasts of the liver [17]; they make up 5-8% of the liver cells. The liver walls are lined up with epithelial mesenchymal cells.

Damage to the liver may be caused by toxic drugs, alcohol abuse, hepatitis B and C, or autoimmune diseases. The damage may trigger reparative and inflammatory responses, and the onset of fibrosis. In homeostasis, the HSCs are quiescent cells, but after early liver damage, they become activated [65], producing hyaluronic acid (HA) [3,4,24,73]. HA promotes fibroblast activation and proliferation [23].

The gold standard for diagnosis and monitoring the progression of liver fibrosis is biopsy. But this procedure is invasive and incurs risk, and cannot be repeated frequently. Hence there is currently a great interest in identifying non-invasive markers for diagnosis of liver fibrosis that can detect the pathological progression of the disease [2,4,10,19,49,58,61]. The present paper develops for the first time a mathematical model of liver fibrosis. The model is a continuation and extension of the authors' articles [32,31] which dealt with kidney fibrosis and lung fibrosis. We use the model to describe the progression of the disease, and identify biomarkers that can be used to monitor treatment for liver fibrosis.

The model is based on the diagram shown in Fig. 1. In Section 2 we list all the variables and proceed to represent the network of Fig. 2 by a system of partial differential equations (PDEs) in a portion of the liver. In Section 3, we simulate the model, and explore potential anti-fibrotic drugs. We also use the mathematical model to quantify the scar density in terms of a combination of two serum biomarkers HA and TIMP, that have been observed in patients [2,4,19].

The conclusion of the paper is given in Section 4, and the parameter estimates used in the simulations are given in Section 5. In Section 5 we also perform sensitivity analysis, and in Section 6 we briefly describe the computational method used in the simulations.

---

2. Mathematical model

In this section we develop a mathematical model of liver fibrosis. The fibrosis takes place in some region $\Omega$ of the liver. The variables used in the model are given in Table 1. These variables satisfy a system of PDEs in $\Omega$.

Equation for macrophage density. The equation for the density of M1 macrophages is given by

$\begin{align*}&\frac{\partial M_1}{\partial t}-D_M{\nabla^2} M_1\\=&\underbrace{-\nabla\cdot({M_1}\chi_{P}\nabla P)}_{chemotaxis}+\underbrace{\lambda_{M_1}\frac{\varepsilon_1}{\varepsilon_1+\varepsilon_2}M_2}_{M_2\rightarrow M_1}-\underbrace{\lambda_{M_2}\frac{\varepsilon_2}{\varepsilon_1+\varepsilon_2}M_1}_{M_1\rightarrow M_2}\underbrace{-d_{M_1}M_1}_{death},\end{align*}$

where

$\varepsilon_1\;\;=\;\;\Big(\lambda_{MI_\gamma}\frac{I_\gamma}{K_{I_\gamma}+I_\gamma}+\lambda_{MT_\alpha}\frac{T_\alpha}{K_{T_\alpha}+T_\alpha}\Big)\frac{1}{1+I_{10}/K_{I_{10}}},\\ \varepsilon_2\;\;=\;\;\lambda_{MI_4}\frac{I_4}{K_{I_4}+I_4}+\lambda_{MI_{13}}\frac{I_{13}}{K_{I_{13}}+I_{13}}.$

The term $-\nabla\cdot(M_1\chi_{P}\nabla P)$ is the chemotactic effect {of} MCP-1 on M1 macrophages; $\chi_{P}$ is the chemotactic coefficient.

M2 macrophages may become M1 macrophages under the influence of IFN-$\gamma$ and TNF-$\alpha$, a process resisted by IL-10 [15,33], and M1 macrophage may become M2 macrophages under the influence of IL-4 and IL-13 [84,85]. Hence the transition between M1 and M2 macrophages depends on the ratio of $\varepsilon_1$ to $\varepsilon_2$: the transition $M_2\rightarrow M_1$ is at rate proportional to $\frac{\varepsilon_1}{\varepsilon_1+\varepsilon_2}$, while the transition $M_1\rightarrow M_2$ is at rate proportional to $\frac{\varepsilon_2}{\varepsilon_1+\varepsilon_2}$. These exchanges of polarization in M1/M2 are expressed by the second and third terms on the right-hand side of the M1 equation.

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 tissue [86,88]. Hence they satisfy the boundary condition

$D_M\frac{\partial M_1}{\partial n}+ \tilde{\beta}(P)M_0=0\hbox{ on the boundary of blood capillaries },$

where $\tilde{\beta}(P)$ depends on MCP-1 concentration, $P$. Here $M_0$ denotes the density of monocytes in the blood, i.e., the source of $M_1$ macrophages from the vascular system. As in [32] we replace the boundary conditions on the blood capillaries by a source term in the tissue, $\beta(P)M_0$.

Then the equation for M1 density in $\Omega$ satisfies the equation

$\frac{\partial M_1}{\partial t}-D_M{{\nabla^2}} M_1\;\;=\;\;\beta(P)M_0-\nabla\cdot(M_1\chi_{P}\nabla P)+\lambda_{M_1}\frac{\varepsilon_1}{\varepsilon_1+\varepsilon_2}M_2\\\;\;\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;-{\lambda_{M_2}\frac{\varepsilon_2}{\varepsilon_1+\varepsilon_2}M_1}- d_{M1}M_1. \label{Mac}$

(1)

We take $\beta(P)=\beta\frac{P}{K_P+P}$, where $\beta$ and $K_P$ are constants.

The M2 macrophage density satisfies the equation

$\frac{\partial M_2}{\partial t}-D_M\nabla^2 M_2\;\;=\;\;A_{M_2}+\underbrace{\lambda_{M_2}\frac{\varepsilon_2}{\varepsilon_1+\varepsilon_2}M_1}_{M_1\rightarrow M_2}\underbrace{-\lambda_{M_1}\frac{\varepsilon_1}{\varepsilon_1+\varepsilon_2}M_2}_{M_2\rightarrow M_1}\underbrace{-d_{M2}M_2}_{death}, \label{Mac1}$

(2)

where $A_{M_2}$ accounts for a source of Kupffer cells. The second and third terms on the right-hand side are complimentary to the corresponding terms in Eq. (1).

Equations for T cells density. Naive T cells, T0, are activated as either Th1 by contact with M1 macrophages in an IL-12 environment [37], a process down-regulated by IL-10 [15], or as Th2 cells by contact with M2 macrophages under an IL-4 environment [93]. IL-2 induces proliferation of Th1 cells [29]. Both activation and proliferation of Th1 cells are antagonized by IL-13 [16], while the activation of Th2 cells is antagonized by Th1 [6,55,93]. Hence, the equations of Th1 and Th2 densities are given as follows:

$\frac{\partial T_1}{\partial t}-D_T\Delta T_1\;\;=\;\; \Big(\lambda_{T_1M_1}T_0\frac{M_1}{K_{M_1}+M_1}\frac{I_{12}}{K_{I_{12}}+I_{12}}\frac{1}{1+I_{10}/K_{I_{10}}}\\ \;\; \;\;\;\;\;\;\;\;\;\;\;\;+\lambda_{TI_{2}}\frac{I_{2}}{K_{I_2}+I_2}T_1\Big)\frac{1}{1+I_{13}/K_{I_{13}}}\underbrace{-d_{T_1}T_1}_{death},\\ \frac{\partial T_2}{\partial t}-D_T\Delta T_2\;\;=\;\; \underbrace{\lambda_{T_2}T_0\frac{M_2}{K_{M_2}+M_2}\frac{I_{4}}{K_{I_{4}}+I_4}\frac{1}{1+T_1/K_{T_1}}}_{activation}$

(3)

$\underbrace{-d_{T_2}T_2}_{death},$

(4)

where $T_0$ is the density of T0.

Equation for TEC density ($E_0$ and $E$) The equation of the inactivated TEC ($E_0$) density is given by

$\begin{align} \frac{d E_0}{d t}=A_{E_0}\Big(1+\underbrace{\frac{\lambda_1E_0I_D}{K_E+E_0I_D}}_{repair}\Big)\underbrace{-\Big(d_{E_0}E_0+\delta+d_{E_0T}\frac{T_\beta}{K_{T_\beta}+T_\beta}\Big)}_{apoptosis}\underbrace{-\lambda_{E_0}E_0I_D}_{E_0\rightarrow E},\label{E0}\end{align}$

(5)

and the equation for the activated TEC (E) is

$\frac{d E}{d t}\;\;=\;\;\underbrace{\lambda_{E_0}E_0I_D}_{activation}{\underbrace{-\lambda_{EM}EI_D}_{EMT}}\underbrace{-d_EE}_{death}.$

(6)

In homeostasis, the production of $E_0$ is represented by the term $A_{E_0}$, and the death rate is represented by $d_{E0}E_0$. $I_D=0$, $\delta=0$ and activated TGF-$\beta$ concentration is very small. The injury to the epithelium is expressed in two ways: (i) by activation of TEC, which is represented by term $\lambda_{E_0}E_0I_D$, where $D$ is the damaged region and $I_D=1$ on $D$, $I_D=0$ elsewhere, and (ii) by increased apoptosis caused by oxidative stress [14,40] (the term with $\delta$) and by TGF-$\beta$ [63,70]. The damaged epithelium is partially repaired by fibrocytes, and this is expressed by the term $\frac{\lambda_1E_0I_D}{K_D+E_0I_D}$ [34]. The second term of the right-hand side in Eq. (6) accounts for epithelial mesenchymal transition (EMT) due to injury [63].

Equations for fibroblast concentration ($f$) and myofibroblast concentration ($m$) The fibroblasts and myofibroblasts equations are given by:

$\frac{\partial f}{\partial t}-D_f\nabla^2 f\;\;=\;\; {\underbrace{\lambda_{Ef}E_0}_{source}}+\lambda_{fH_A}\frac{H_A}{K_{H_A}+H_A}f\\ \;\;\ \;\;\;\;\;\;\;\;\;\;+{\underbrace{\lambda_{fE}\Big(\frac{T_\beta}{K_{T_\beta}+T_\beta}+\frac{I_{13}}{K_{I_{13}}+I_{13}}\Big)\frac{E}{K_E+E}f}_{production}}\\ \;\;\;\;\;\;\;\;\;\;\ \;\;\underbrace{-\Big(\lambda_{mfT}\frac{T_\beta}{K_{T_\beta}+T_\beta} +\lambda_{mfG}\frac{G}{K_G+G}\Big)f}_{f\rightarrow m}\underbrace{-d_{f}f}_{death},$

(7)

$\frac{\partial m}{\partial t}-D_m\nabla^2 m\;\;=\;\;\underbrace{\Big(\lambda_{mfT}\frac{T_\beta}{K_{T_\beta}+T_\beta} +\lambda_{mfG}\frac{G}{K_G+G}\Big)f}_{f\rightarrow m}\underbrace{-d_M m}_{death}.$

(8)

The first term on the right-hand side of Eq. (7) is a source from $E_0$-derived fibroblast growth factor (bFGF), which for simplicity we take to be in the form $\lambda_{Ef}E_0$ [9,70]. The second term represents the activation and proliferation of fibroblasts by hyaluronic acid [23]. The third term on the right-hand side of Eq. (7) accounts for the fact that TGF-$\beta$ and IL-13, combined with E-derived bFGF, increase proliferation of fibroblasts [9,12,21,36,52]. For simplicity, we do not include bFGF specifically in the model, and instead represent it by {$E$}. As in [31,32], TGF-$\beta$ and PDGF transform fibroblasts into myofibroblasts [20,52,59,66,77,82,92] (the fourth term on the right-hand side of Eq. (7)).

Equation for HSC ($H$) HSC is activated by PDGF and TGF-$\beta$ [50]. Hence,

$\frac{\partial H}{\partial t}-D_H\nabla^2 H\;\;=\;\;A_H+\underbrace{\Big(\lambda_{HG}\frac{G}{K_G+G}+\lambda_{H T_\beta}\frac{T_\beta}{K_{T_\beta}+T_\beta}\Big)H}_{production}\underbrace{-d_HH}_{death}.$

(9)

Equation for HA ($H_A$) HA is produced by HSCs and degraded by sinusoidal epithelial cells [3,4,19,24,73]. Hence,

$\frac{\partial H_A}{\partial t}-D_{H_A}\nabla^2 H_A\;\;=\;\;\underbrace{\lambda_{H_A}H}_{production\;}\underbrace{-d_{H_A}H_A}_{degradation}.$

(10)

Equation for ECM density ($\rho$) and scar ($S$). The ECM consists primarily of fibrillar collagen and elastin, but includes also fibronectin, lamina and nitrogen that support the matrix network by connecting or linking collagen (Lu et al. [53]). For simplicity, we represent the ECM by the density of collagen. ECM is produced by fibroblasts, myofibroblasts [52,59,60,76,82,92] and by HSCs whose production rate of ECM is similar to that of myofibroblast[22], and TGF-$\beta$ enhances the production of ECM by myofibroblasts [43,47,52,59,82,92]. MMP ($Q$), in fibrosis, degrades collagen by cutting the protein into small fragments (Veidal et al. [83]); we assume that the loss of collagen is proportional to $Q\rho$. The equation for the density of ECM is then given by:

$\frac{\partial \rho}{\partial t}\;\;=\;\;\underbrace{\lambda_{\rho f}f\big(1-\frac{\rho}{\rho_0}\big)^++\lambda_{\rho m}\Big(1+\lambda_{\rho T_\beta}\frac{T_\beta}{K_{T_\beta}+T_\beta}\Big)(m+H)}_{production}\\ \;\;\ \;\;\;\;\;\;\;\;-d_{\rho Q}Q\rho\underbrace{-d_\rho\rho}_{death},$

(11)

where $\big(1-\frac{\rho}{\rho_0}\big)^+=1-\frac{\rho}{\rho_0}$ if $\rho<\rho_0$, $\big(1-\frac{\rho}{\rho_0}\big)^+=0$ if $\rho\geq\rho_0$. There are several computational models of collagen network under various biological conditions (Lee et al. [45]) and under strain-dependent degradation (Hadi et al. [25]). But the parameters used in these models are not helpful in determining how scar develops by excess of ECM while under the effect of MMP. Since MMP increases scarring in cases of excessive collagen concentration, we shall use the following simple formula to describe the the growth of a scar:

$S=\lambda_S(\rho-\rho^*)^+\Big(1+\lambda_{SQ}\frac{Q}{K_Q+Q}\Big),$

(12)

where $\rho^*$ is the concentration of collagen in normal healthy tissue and $\lambda_S$, $\lambda_{SQ}$ and $K_Q$ are positive parameters.

Equation for MCP-1 ($P$) The MCP-1 equation is given by

$\frac{\partial P}{\partial t}-D_{P}\nabla^2 P\;\;=\;\;\underbrace{{\lambda_{PE}E}}_{production}\underbrace{-d_{PM}\frac{P}{K_P+P}M_1-d_PP}_{degradation},$

(13)

where $\lambda_{PE}$ represents the growth rate by activated TEC following damage to the endothelium [32,34,35,71,72,75,88,87]. The second term on the right-hand side accounts for the internalization of MCP-1 by macrophage, which is limited due to the limited rate of receptor recycling.

Equations for concentrations of PDGF ($G$), MMP ($Q$), and TIMP ($Q_r$). These cytokines are produced by macrophages [86,94] and, as in [32], the following sets of diffusion equations hold for $G$, $Q$ and $Q_r$:

$\frac{\partial G}{\partial t}-D_{G}\nabla^2 G\;\;=\;\; \underbrace{\lambda_{GM}M_2}_{production\;}\underbrace{-d_{G}G}_{degradation},$

(14)

$\frac{\partial Q}{\partial t}-D_{Q}\nabla^2 Q\;\;=\;\; \underbrace{{\lambda_{QM}M_2}}_{production}\underbrace{-d_{QQ_r}Q_rQ}_{depletion}\underbrace{-d_{Q} Q}_{degradation},$

(15)

$\frac{\partial Q_r}{\partial t}-D_{Q_r}\nabla^2 Q_r\;\;=\;\;\underbrace{{\lambda_{Q_rM}M_2}}_{production}\underbrace{-d_{Q_rQ}QQ_r}_{depletion}\underbrace{-d_{Q_r}Q_r}_{degradation}.$

(16)

Note that in Eq. (15), MMP is lost by binding with TIMP (second term), a process which also depletes TIMP in Eq. (16).

Equations for concentrations of TGF-$\beta$ ($T_\beta$) and TNF-$\alpha$ ($T_\alpha$). TGF-$\beta$ is produced and is activated by M2 macrophages, a process jointly enhanced by IL-13 [12,21,36,80]; in addition, TGF-$\beta$ is produced and is activated by TECs and fibroblasts [8,70]. Hence $T_\beta$ satisfies the equation:

$\begin{align}\frac{\partial T_\beta}{\partial t}-D_{T_\beta}\nabla^2 T_\beta= \underbrace{{\lambda_{T_\beta M}M_2\Big(1+\lambda_{T_\beta I_{13}}\frac{I_{13}}{I_{13}+K_{I_{13}}}\Big)}+{\lambda_{T_\beta f}f\frac{E}{E+K_E}}}_{production}\underbrace{-d_{T_\beta}T_\beta}_{degradation}. \label{15}\end{align}$

(17)

TNF-$\alpha$ is produced by M1 macrophages [67], and by TEC [8,57], and is depleted when it combines with receptors on M2 in the process which produces phenotype exchange $M_2\rightarrow M_1$:

$\frac{\partial T_\alpha}{\partial t}-D_{T_\alpha}\nabla^2 T_\alpha\;\;=\;\; \underbrace{\lambda_{T_\alpha M}M_1+{\lambda_{T_\alpha E}E}}_{production}\underbrace{-\lambda_{MT_\alpha}\frac{T_\alpha}{K_{T_\alpha}+T_\alpha}M_2}_{M_2\rightarrow M_1}\underbrace{-d_{T_\alpha}T_\alpha}_{degradation}.\label{16}$

(18)

Equation for interleukins: IL-2 is produced by Th1 cells [29]:

$\frac{\partial I_{2}}{\partial t} -D_{I_{2}}\Delta I_{2}\;\;=\;\;\underbrace{\lambda_{I_{2} T_1}T_1}_{production}\underbrace{-d_{I_{2}}I_{2}}_{degradation}.\label{I2}$

(19)

IL-4 is produced by Th2 cells and M2 macrophages [6,55], hence

$\frac{\partial I_{4}}{\partial t} -D_{I_{4}}\Delta I_{4}\;\;=\;\;\underbrace{\lambda_{I_{4}T_2 }T_2+\lambda_{I_{4}M_2}M_2}_{production}\underbrace{-d_{I_{4}}I_{4}}_{degradation}.\label{I4}$

(20)

IL-10 is produced primarily by M2 macrophages, while IL-12 is produced primarily by M1 macrophages in a process that is s antagonized by IL-10 [15] and IL-13 [1]. Hence IL-10 and IL-12 satisfy the equations:

$\frac{\partial I_{10}}{\partial t}-D_{I_{10}}\Delta I_{10}\;\;=\;\;\underbrace{\lambda_{I_{10}M_{2}}M_{2}}_{production}\underbrace{-d_{I_{10}}I_{10}}_{degradation},\\ $

(21)

$\frac{\partial I_{12}}{\partial t}-D_{I_{12}}\Delta I_{12}\;\;=\;\;\underbrace{\lambda_{I_{12}M_{1}}M_{1}\frac{1}{1+I_{10}/K_{I_{10}}}\frac{1}{1+I_{13}/K_{I_{13}}}}_{production}\underbrace{-d_{I_{12}}I_{12}}_{degradation}.\\ $

(22)

IL-13 is produced by M2 macrophages [28,68] and by Th2 cells [79,84], so that

$\frac{\partial I_{13}}{\partial t} -D_{I_{13}}\Delta I_{13}\;\;=\;\;\underbrace{\lambda_{I_{13}T_2 }T_2+\lambda_{I_{13}M}M_2}_{production}\underbrace{-d_{I_{13}}I_{13}}_{degradation}.\label{I13}$

(23)

---

2.1. Equations for IFN-$\gamma$

($I_\gamma$): IFN-$\gamma$ is produced by Th1 cells [29]:

$\frac{\partial I_\gamma}{\partial t}-D_{I_\gamma}\Delta I_\gamma\;\;=\;\;\underbrace{\lambda_{I_\gamma T_1}T_1}_{production}-d_{I_\gamma}I_\gamma.$

(24)

The parameters which appear in Eqs (1)-(21) are listed in Tables 2-4 of Section 5 together with their dimensional values.

---

3. Results

A model of of renal fibrosis was introduced by Hao et al. [32]. The model combines M1 and M2 macrophages into one variable $M$, and does not include T cells. Based on patients' data, the model suggests that urine measurements of (TGF-$\beta$, MCP-1) could serve as biomarkers to determine the severity of the disease. The lung contains many tiny alveoli, where oxygen is absorbed. The tissue surrounding them is the lung interstitium. Idiopathic pulmonary fibrosis (IPF) is a fibrosis of the interstitium whose etiology is unknown. A model of IPF was developed by Hao et al. [31] takes into account of the unique alveoli structure of the lung. The model includes the alveolar macrophages (M2) and the proinflammatory macrophages (M1) but, it does include T cells.

The model developed in the present paper is more comprehensive than the models developed in [31,32] since it includes T cells, and liver-specific cells, namely HSCs, as well as the hyaluronic acid (HA), produced by HSCs; both HSCs and HA play important roles in liver fibrosis. By including T cells and HSCs we can explore potential anti-fibrotic drugs such as injection of IFN-$\gamma$, and potential serum biomarkers such as HA.

Boundary conditions. We assume that fibrosis occurs only within the region $\Omega$, hence:

$\rm{all\;\;the\;\;variables\;\;satisfy\;\;non-flux\;\;condition\;\;on\;\;the\;\;boundary\;\;of\;\;}\Omega.$

(25)

Initial conditions. We take the following initial conditions (mostly from [32,33]):

$\small M_1=3.73\times10^{-5}~ g/ml, ~M_2=3.38\times10^{-5}~ g/ml,~T_1=4.83\times10^{-5} ~ g/ml,\\ \;\;\;\;\;\;T_2=2.37\times10^{-5}~ g/ml, ~E_0=0.1~ g/ml, ~E=1\times10^{-6}~ g/ml, \\ \;\;\;\;\;\;f=1.2\times10^{-2}~g/ml, ~m=7.1\times10^{-6}~ g/ml,~\rho=0.002~ g/ml\\ \;\; ~P=5.59\times10^{-8}~ g/ml, ~G=3.07\times10^{-10}~ g/ml, ~Q=2.29\times10^{-6}~ g/ml, \\ \;\;\;\;\;\;Q_r=10^{-6}~ g/ml, ~T_\beta=1.52\times10^{-9}~ g/ml, ~T_\alpha=1.47\times10^{-9}~ g/ml,\\ \;\; I_2= 2.49\times10^{-8}~ g/ml,~I_4=3.22\times10^{-12}~ g/ml,~I_{13}=1.13\times10^{-9}~ g/ml,\\ I_{10}=7.66\times10^{-12}~ g/ml, ~I_{12}=1.64\times10^{-8}~ g/ml, \hbox{~and~}I_\gamma=1.82\times10^{-11}~ g/ml.$

(26)

We also assume initial homeostasis with a small amount of inflammation represented by the term $\lambda_{PE}E$ in Eq. (13):

$\lambda_{E_0}E_0I_D=0,~~\lambda_{PE}E=\varepsilon_0, \varepsilon_0\hbox{ is small}.$

(27)

Finally, we take at $t=0$ homeostasis values of HA (from [7,19]) and H (from Sec. 5, under Eq. (9)):

$HA=10^{-4}~ g/ml ,~ H=0.001~ g/ml.$

(28)

In the following simulations the parameter values are taken from Tables 2-4. For simplicity, we simulate the model for a 2-d domain $\Omega$, taking

$\Omega \hbox{ a square of side 1 cm, and~} D \hbox{ a concentric square of side 0.3 cm.}$

(29)

Fig. 3 shows the dynamics of the average concentrations of cells, cytokines and ECM for the first 200 days. The parameters are taken from Tables 2-4.

We see that most densities/concenrations nearly stabilize by day 100; however TEC density and fibroblast concentration continue to decrease while the myofibroblasts concentration increases. We note that ECM increases up to 6 times its initial value for healthy case, in agreement with [19].

We are mostly interested in scar formation, hence in scar density $S$. The parameters in Eq. (12) are unknown, and for illustration we take $\lambda_S=100$, $\lambda_{SQ}=1$ and $K_Q=5\times10^{-6}$ in Eq. (12). The profile of the scar density $S(t)$ for the first 200 days is shown by the blue curve in Fig. 4. We see that $S(t)$ grows initially fast, but the growth rate gradually decreases. Other choices of the parameters $\lambda_S$, $\lambda_{SQ}$ and $K_Q$ show the same qualitative behavior.

Treatment studies. We can use the model to explore potential drugs. We express the effect of a drug indirectly by either reducing some of the parameters in the relevant equations by factors such as $\frac{1}{1+A}$, $\frac{1}{1+B}$, $\theta$, or by adding a constant term $c$ in the relevant equation during the treatment period. The choices of $A$, $B$ and $c$ are somewhat arbitrary, since they depend on the actual amount of dozing. Such drugs could be, for instance, anti-TGF-$\beta$, NOX inhibitor or IFN-$\gamma$ injection. Fig. 4 displays the effect of treatment when the drug is administered at day 100 for 100 days, continuously. We determine the efficacy of the drug by how much it blocks or reduces the scar density.

Anti-TGF-$\beta$. We first consider an anti-TGF$\beta$ drug, such as Pirfenidone which was recently approved in the United States. In our model we need to replace $\lambda_{T_\beta M}$ and $\lambda_{T_\beta f}$ by $\lambda_{T_\beta M}/(1+A)$ in Eq. (17) and $\lambda_{T_\beta f}/(1+A)$, and $T_\beta$ by $T_\beta/(1+B)$ in all terms where $T_\beta$ acts to promote fibrosis. The green curve in Fig. 4 shows the effect of the drug for $A=B=0.1$. We see that in terms of scar, the drug is initially effective in decreasing the scar density, but in the long term its effect diminishes.

NOX inhibitor. One of suggested novel drugs for treatment of hepatic fibrosis is NOX inhibitor [42]. NOX are membrane proteins that activate HSCs [65]. The effect of anti NOX drug is to decrease $\lambda_{HG}$ and $\lambda_{HT_\beta}$ in Eq. (9) by a factor of $\theta\in(0,1]$; $\theta=1$ when no drug is applied. The black curve in Fig. 4 shows the dynamics of the scar for $\theta=0.5$ and suggests that the drug may be very effective as anti-fibrotic drug. Micro RNA-21 (miR-21) modulates ERK1 signaling in HSCs activation and is overexpressed in hepatic fibrosis [95]. Anti miR-21 is a potential anti-fibrotic drug which, in our model, has the same effect as NOX inhibitor.

Injection of IFN-$\gamma$. It was suggested by Weng et al. [90] that IFN-$\gamma$ treatment may reduce liver fibrosis. Such a treatment means that we need to add in our model a source term $c$ in Eq. (24) to represent the injection of IFN-$\gamma$. Taking $c=10^{-9}$ g/ml/day, we see, in Fig. 4, that the drug initially promotes fibrosis but later on it reverses fibrosis. This behavior may be attributed to the fact that fibrosis has both non-inflammatory reparative aspect and proinflammatory aspect. Hence IFN-$\gamma$ injection can affect the disease in either negative and positive ways.

Biomarkers. Patients with liver fibrosis have higher concentration of HA and TIMP in the liver [2,4,19]. We can use the mathematical model to develop a diagnostic tool to determine the state of the disease based on combined measurements of HA and TIMP. We do not know when the disease of an individual patient began, or equivalently, what was the damaged area $D$ of an individual patient at time $t=0$ when the biomarkers were measured. Hence we take $D$ to be a rectangle with variable side $\lambda$, where $\lambda\in[0.2,0.9]$ depends on the individual patient.

For each $\lambda$ we simulate the model for time $0\leq t\leq 200$ days and determine the quantities $HA(t,\lambda)$, $TIMP(t,\lambda)$ and the scar density $S(t,\lambda)$. As $\lambda$ increases, the curves $\Gamma(\lambda)=\{HA(t,\lambda),TIMP(t,\lambda), 0\leq t\leq 200\}$ increase and span a region in the $HA-TIMP$ plane shown in Fig. 5. We associate to each point in this region the corresponding value of $S(t,\lambda)$, using color from the color column in Fig. 5. Fig. 5 can then be used to determine, for any individual, based on his/her concentrations of HA and TIMP in the liver, what the scar density is.

As reported in [26,61,62], the serum biomarkers of HA and TIMP reflect the disease state, and thus roughly the tissue levels of HA and TIMP. As the correlation between tissue and serum concentrations of HA and of TIMP become more precise, Fig. 5 could then provide a quantitative non-invasive diagnostic tool for liver fibrosis.

We note however that some of the parameters in the model equations may not be sufficiently precise; there are also variations from one person to another. Sensitivity analysis (such as that carried in Sec 5.1) shows that the scar density varies in a continuous way when parameters are changed continuously within a limited range. Hence Fig. 5 should be viewed as just one possible prediction map; similar maps could be produced with other parameters. When new experimental and clinical data become available, some of the parameters, especially these under "estimated" in Tables 2-4, could be modified to make simulations better fit the data. Sensitivity analysis could be used in order to modify collectively a group of parameters.

---

4. Conclusion.

Fibrosis in an organ is characterized by excessive deposition of fibrous connective tissue. It disorganizes the architecture of the organ, leading to the formation of scars and eventual disfunction and failure of the organ. There are currently no drugs that can appreciably reverse the progress of the disease. The present paper focuses on liver fibrosis. The gold standard for diagnosis and monitoring the pathological progression of liver fibrosis is biopsy. But this procedure is invasive and incurs risk, and cannot repeated frequently. For this reason we developed, for the first time, a mathematical model that describes the progression of the disease and the effect of drug treatment, and we used the model to construct a diagnostic map based on a combination of biomarkers. The model is represented by a system of 24 partial differential equations for the concentrations of cells and cytokines. The cells are macrophages M1 and M2, T cells Th1 and Th2, fibroblasts, myofibroblasts, HSCs, and tissue epithelial cells. The cytokines are either produced by these cells, or affect the activities of the cells. The mathematical model builds on the models developed in [31,32], but it also includes HSCs and CD4+ T cells: Th1 and Th2. This extended model enables us to explore new potential drugs. For example, we tested with our model the efficacy of treatment by injection of IFN-$\gamma$, a suggestion made in [90]. We found, interestingly, that the drug initially increases fibrosis but later on decreases it.

We used the model to explore the efficacy of other potential drugs aimed to block liver fibrosis. Currently, most of the available data on anti-fibrotic drugs are obtained from mice experiments. As more clinical data become available, our model could be refined (by modifying some of the parameters) and validated, and it could then serve as as useful tool in exploring the efficacy of anti-fibrotic drugs for the treatment of liver fibrosis in human patients.

There is currently a great interest in determing reliable serum biomarkers for diagnosis and prognosis of liver fibrosis [2,4,10,19,49,58,61]. Our mathematical model can be used as diagnostic and prognostic tool by using a combination of two biomarkers. Thus, in Fig. 5 we quantified the dependence of scar density in liver fibrosis in terms of concentrations of TIMP and HA in the fibrotic tissue; these two concentrations are overexpressed in serum of patients with liver fibrosis [19,61]. Our model can be used to explore other combinations of biomarkers in liver fibrosis as more experimental and clinical data become available.

---

5. Parameters

The parameters of the model are listed in Tables 2-4. Most of the parameter are taken from previous works [29,30,30,32,33]. The remaining parameters are estimated below.

Eq. (3). In the blood of a healthy adult there are 5000,000-750,000 T cells per ml, which translates into an average of approximately $5\times10^{-4}$g/ml. We assume that the density of naive CD4+ T cells in the tissue is significantly smaller, taking $T_0=3\times10^{-5}$ g/ml.

Eq. (4). The production of $T_1$ by $M_1$ under $I_{12}$ environment is $\lambda_{T_1M_1}=10$/day [33]. We assume that the production of $T_2$ by $M_2$ under $I_4$ environment is much smaller, taking $\lambda_{T_2}=0.75$/day.

Eq. (5). We assume that the density of $E_0$ of the inactivated epithelial cells in the liver is 0.1 g/ml. The repair term $A_{E_0}$ was estimated in [32] to be $8.27\times10^{-3}$ g/ml/day. We assume that the repair is 5 times slower in the liver, taking $A_{E_0}=1.65\times10^{-3}$ g/ml/day.

Eq. (7). The production rate of fibroblasts by activated TEC is $5\times10^{-4}$/day [32]. We assume that the production of fibroblasts by HSC-produced HA is five times larger, taking, $\lambda_{fH_A}=2.5\times10^{-3}$/day. The transition rate from fibroblasts to myofibroblasts, $\lambda_{mfT}$ and $\lambda_{mfG}$, in the lung were estimated in [32] by 0.12/day. We assume these rates are larger in the liver than in the lung, taking $\lambda_{mfT}=\lambda_{mfG}=0.3$/day.

Eq. (9). HSCs make up 5-8% of all the liver tissue. Accordingly we take $H=0.02$ g/ml. We assume that in homeostasis only 5% HSCs are activated, that is, $H=0.001$ g/ml. The death rate of HSCs is assumed to be the same as for fibroblast, $d_H=d_f=1.66\times10^{-2}$ day$^{-1}$. From the steady state equation $A_H-d_HH=0$, we then get, $A_H=1.66\times10^{-5}$ g/ml/day.

We take $\lambda_{HG}=\lambda_{HT_\beta}$ and assume that, when activated, the number of HSC increased by 25%. We account for this by taking $\lambda_{GH}=\lambda_{T_\beta H}=0.2d_H=3.32\times10^{-3}$ day$^{-1}$.

Eq. (10). We assume the degradation rate of HA by sinusoidal epithelial cells to be $d_{H_A}=0.29$/day [7]. In health the concentration of HA is $10^{-4}$ g/ml [7,19]. From the steady state equation

$\lambda_{H_A}H-d_{H_A}H_A=0,$

with $H=0.001$, we then get $\lambda_{H_A}=2.9\times10^{-2}$ day$^{-1}$.

Eq. (15) The production of MMP and TIMP by M2 macrophage in the lung was taken in [33] to be $3\times10^{-4}$/day and $6\times10^{-5}$/day, respectively. We assume that the production rate is larger in the liver, taking $\lambda_{QM}=3\times10^{-3}$/day and $\lambda_{Q_rM}=6\times10^{-4}$/day.

---

5.1. Sensitivity analysis

We performed sensitivity analysis on some of the production parameters of the system (1)-(17). Following the method in [56], we performed Latin hypercube sampling and generated 1000 samples to calculate the partial rank correlation (PRCC) and the p-values with respect to the scar concentration at day 200. The results are shown in Fig. 6 (The p-value was $< 0.01$).

Scar density grows if $\rho$ and $Q$ are increased (Eq. (12) and $\rho$ increases with increase in $f$, $m$, $H$ and $T_\beta$ (Eqs. (11), (12)); $f$ is increased by $H_A$ which is produced by $H$. These observations explain why the parameters $\lambda_{HG}$, $\lambda_{HT_\beta}$, $\lambda_{H_A}$, $\lambda_{T_\beta M}$ and $\lambda_{T_\beta I_{13}}$ are positively correlated. We next observe that $T_\beta$ is produced by $M_2$ (and $f$), hence the transition $M_1\rightarrow M_2$ positively affects scar growth. This transition is increased if $T_\alpha$ is decreased while $I_{10}$ and $I_{13}$ are increased (see the form of $\varepsilon_1$, $\varepsilon_2$ which appear in Eqs. (1), (2)). Hence $\lambda_{I_{10}M_2}$, $\lambda_{I_{13}T_2}$ and $\lambda_{M T_\alpha}$ are positively correlated while $\lambda_{T_\alpha M}$ is negatively correlated (see Eq. (18)). Since MCP-1 attracts macrophages to the liver, the production rate of MCP-1 by TEC, $\lambda_{PE}$, is positively correlated. The sensitivity analysis can be carried out in a similar way for the remaining production parameters.

---

6. Computational method

In order to illustrate our numerical method, we consider the following diffusion equation:

$\frac{\partial X}{\partial t}-D_X\nabla^2 X=F_X \hbox{ in }\Omega,\label{NS}$

(30)

where the right-hand side accounts for all the 'active' terms. Let $X_{i,j}^n$ denote a numerical approximation of $X(ih_x,jh_y,n\tau)$, where $h_x$ and $h_y$ are the stepsize in the $x$ and $y$ directions respectively, and $\tau$ is the time stepsize. Then a discretization is derived by the explicit Euler five-point finite difference scheme, i.e.,

$\begin{align}&\frac{X_{ij}^{n+1}-X_{ij}^N}{\tau}-D_X\Big(\frac{X_{i+1,j}^n+X_{i-1,j}^n-2X_{i,j}^n}{h_x^2}+\frac{X_{i,j+1}^n+X_{i,j-1}^n-2X_{i,j}^n}{h_y^2}\Big)\nonumber\\=&F_X(X_{i,j}^n) \hbox{ in }\Omega.\label{NS1}\end{align}$

(31)

In order to make the scheme stable, we take $\tau\leq \frac{h^2}{4D_X}$, namely $\tau= 0.1\frac{h^2}{D_X}$, where $h=h_x=h_y$.

---

Acknowledgments

This work has been supported by the Mathematical Biosciences Institute and the National Science Foundation under Grant DMS 0931642.

---