Research Article

Mathematical model of sarcoidosis

See allHide authors and affiliations

PNAS November 11, 2014 111 (45) 16065-16070; first published October 27, 2014; https://doi.org/10.1073/pnas.1417789111

  1. Contributed by Avner Friedman, September 17, 2014 (sent for review August 8, 2014; reviewed by Marc Judson and Reinhard C. Laubenbacher)

Significance

Sarcoidosis is a disease involving abnormal collection of granulomas that develop in the lung and other organs. The origin of the disease is unknown, clinical data are very limited, and there is no current effective treatment. This paper develops a mathematical model with simulations that are validated by the available clinical data. The model is then used to explore potential treatments of the disease, to suggest therapeutic targets that may reduce the disease activity, and thus to predict treatment responses in preclinical settings.

Abstract

Sarcoidosis is a disease involving abnormal collection of inflammatory cells forming nodules, called granulomas. Such granulomas occur in the lung and the mediastinal lymph nodes, in the heart, and in other vital and nonvital organs. The origin of the disease is unknown, and there are only limited clinical data on lung tissue of patients. No current model of sarcoidosis exists. In this paper we develop a mathematical model on the dynamics of the disease in the lung and use patients’ lung tissue data to validate the model. The model is used to explore potential treatments.

Sarcoidosis is a disease involving an abnormal collection of inflammatory cells that can interact to form nodules, called granulomas, which are capable of altering the functions of affected tissues and organ systems. The granulomas contain macrophages, T lymphocytes whose functions are regulated by inflammatory mediators such as TNF-α, IFN-γ, IL-2, IL-10, IL-12 and TGF-β. The sarcoidosis granulomas are most commonly detected in the lungs and mediastinal lymph nodes; however, recent technological advances have improved disease detection in other organs, such as the heart, which is now recognized to be involved in one-third of cases, and other vital and nonvital organs (1⇓⇓–4). Sarcoidosis of the lungs and heart contributes to disability and increased mortality in these patients.

The primary cause of sarcoidosis remains a mystery, and progress has been limited by the lack of relevant disease models. It is unknown to what extent genetic predisposition or specific environmental exposures (e.g., antigens derived from infectious organisms) trigger the inflammatory immune response. It is reasonable, however, to assume that, in the lung, inflammation is initiated following inhalation of an environmental antigen, which leads to a typical Th1 immune response that is initiated by macrophages. Activated macrophages secrete proinflammatory cytokines such as IL-12 (5) and TNF-α (6, 7) and anti-inflammatory cytokine IL-10 (8) and IL-13 (9); they and Th17 cells secrete chemokine (C-C motif) ligand 20 (CCL20) CCL20 (10, 11). The CD4+ T cells in sarcoidosis are primarily Th1, Th17, and Treg. Th1 is activated by IL-12, and activated Th1 cells produce IFN-γ, which further activates macrophages; these processes are inhibited by IL-10 (12, 13). Cytokine CCL20 chemoattracts both Treg and Th17 cells (14) into the granuloma. Treg and Th17 are both activated by TGF-β (15). IL-2 secreted by Th1 (16) increases the proliferation of Th1 cells (16), blocks the proliferation of Th17 cells (17), and enhances the activation of Treg by TGF-β (18, 19); TGF-β is secreted by activated macrophages and Treg (20, 21).

A detailed diagram of the network involving the interactions among all these immune cells and cytokines, including GM-CSF production by macrophages, is shown schematically in Fig. 1.

Fig. 1.
Fig. 1.

Schematic network of sarcoidosis. Arrowhead means production or activation, rectangle means inhibition, and oval means chemoattraction.

No current models of sarcoidosis exist. In this paper we develop a mathematical model of sarcoidosis based on the diagram shown in Fig. 1. The model is represented by a system of partial differential equations. Parameters are determined by using the clinical data on cytokine levels in healthy lung tissue as reported in Crouser et al. (22). The model is then validated by data on cytokine levels in lung tissue of patients (22).

We use the model to explore the effect of anti–TNF-α (currently in use) and potential drugs, anti–IL-12, anti–IFN-γ, and TGF-β enhancement, in decreasing the size of sarcoid granulomas.

Mathematical Model

The variables of the model are listed in Table 1. Here we assume that the granuloma occupies a region that varies in time and that macrophages and T cells are in movement with velocity u within the granuloma. The need to use a spatial model arises from the facts that granulomas are regions that evolve in time, and chemotaxis by chemokine CCL20 plays an important role in attracting both Th17 and Treg cells.

View this table:
Table 1.

The variables of the model

We assume that all species are dispersing or diffusing, in the granuloma, with appropriate diffusion coefficients. The equation for each species Xi (1≤i≤k) has the form∂Xi∂t−θ∇⋅(uXi)−DXi∇2Xi=FXi(X1,⋯,Xk),[1]where ∇2 is the Laplace operator, DXi is the diffusion coefficient, and FXi is a function that depends on all of the species and expresses the result of their interactions on the growth of Xi. The term DXi∇2Xi for cells means dispersion, which decreases crowding, and for cytokines it means diffusion. Because cells are much larger than cytokines, their dispersion coefficients are much smaller than the diffusion coefficients of cytokines. The term θ∇⋅(uXi) represent movement by advection: θ=1 for cells but θ=0 for cytokines (cytokines are very small and move only by diffusion). In determining the structure of FXi, we use, for simplicity, the linear mass conservation law. For instance, if Xj+Xk→Xi, then FXi=mXjXk, where m is the production rate of Xi. However, this law will apply only when Xj and Xk are unlimited. If, for example, Xj represents cells and Xk represents molecules that are bound and internalized by Xj, then the internalization of Xk may be limited due to the limited rate of receptor recycling. In this case, we use the Michaelis–Menten law FXi=mXj(Xk/(K+Xk)); we do not use Hill’s law FXi=mXj(Xkn/(K+Xkn)) (n≥2) because we want to keep the linear conservation law for small concentrations.

Equation for Activated Macrophages (MA).

Alveolar macrophages are M2 macrophages (23). They are activated by IFN-γ (24), GM-CSF (25), and TNF-α (26).

The density of the activated macrophages follows the equationdMAdt−∇⋅(uMA)︸transport−DMA∇2MA︸dispersion= (f︸inflammation+ λMIγIγIγ+KIγ︸activation by IFN-γ+λMGGG+KG︸activation by GM-CSF+ λMTαTαTα+KTα︸activation by TNF-α)M0−dMAMA,︸death[2]where f indicates the inflammation. As mentioned above, the use of the Michaelis–Menten law, for instance, λMIγ(Iγ/(Iγ+KIγ)), expresses the fact that internalization of Iγ may be limited due to the limited rate of receptor recycling. The Michaelis–Menten law is similarly used throughout this paper.

Equation for Th1 Cells (T1).

The density of Th1 cells satisfies the equation∂T1∂t−∇⋅(uT1)−DT∇2T1= (λTI12I12KI10T1+I10MA︸activation by macrophage+λTI2I2KI2+I2T1︸proliferation by IL-2)11+Tr+KTr︸inhibition by Treg−dTT︸death.[3]

In the first term on the right-hand side, Th1 cells are activated by IL-12 and direct contact with MHCII of activated macrophages and inhibited by IL-10 (12, 13). IL-2 increases the proliferation of Th1 cells (16). These processes are inhibited by Treg (27). Internalization of IL-12 is much smaller than the inhibition by IL-10, and therefore the term I12/KI12 was neglected in the denominator of the first term of the right-hand side (i.e., KI12/KI10T1=750, in SI Text).

Equation for Treg cells (Tr).

The density of Treg cells satisfies the equation∂Tr∂t−∇⋅(uTr)−DR∇2Tr =−∇(χCR∇C)︸chemotaxis−dTrTr︸death+ λTrI2I2KI2+I2Tβ1+Tα/KTαMAMA+KMA︸activation.[4]Treg cells are activated by TGF-β in the presence of IL-2 (18, 19) and direct contact with activated macrophages and chemoattracted by CCL20 (14). TNF-α resists activation of Treg (28), and this resistance is assumed to be significantly higher than the internalization of TGF-β (which is neglected here).

Equation for Th17 Cells (T17).

Th17, in direct contact with activated macrophages, is activated by TGF-β (15) and other cytokines including IL-6, IL-21, and IL-23 (29). For simplicity we include only TGF-β in our model and accordingly adjust its activation rate λT17. This activation is resisted by Treg (via IL-6) (30) as well as by IL-2 (17). Th17 is also chemoattracted by CCL20 (14). Hence the equation for the density of Th17 cells is given by∂T17∂t−∇⋅(uT17)−DH∇2T17 =−∇(χCH∇C)︸chemotaxis−dT17T17︸death+λT17Tβ(1+Tr/KTr)(1+I2/KI2)MAMA+KMA︸activation by TGF-β resisted by IFN-γ and IL-2.[5]

Equation for IFN-γ (Iγ).

The concentration of IFN-γ is modeled by∂Iγ∂t−DIγ∇2Iγ=λIγT1T1+λIγMMA︸production−dIγIγ︸degradation.[6]IFN-γ is secreted by Th1 lymphocytes and by activated macrophages (24).

Equation for TGF-β (Tβ).

TGF-β is secreted by Treg (21), activated macrophages (20), and Th1 lymphocytes (21). Hence the TGF-β concentration satisfies the equation∂Tβ∂t−DTβ∇2Tβ=λTβTrTr+λTβMMA+λTβT1T1︸production−dTβTβ︸degradation.[7]

Equation for IL-12 (I12).

IL-12 comes in two forms, IL-12 p40 and IL-12 p70 (5). Both forms are produced by activated macrophages and inhibited by IL-10 (31). This process is enhanced by IFN-γ, but the production of IL-12 p70 is negligible without the participation of IFN-γ (5). Hence the equations for the concentration of IL-12 are given by∂I1240∂t−DI12∇2I1240=λI1240MAMA1+I10/KI10︸production by macrophage1+(λI1240IγIγKIγ+Iγ︸enhanced by IFN-γ)−dI12I1240︸degradation,∂I1270∂t−DI12∇2I1270=λI1270MAMA1+I10/KI10IγKIγ+Iγ︸production−dI12I1270︸degradation,I12=I1240+I1270.[8]

Equation for TNF-α (Tα).

The concentration of TNF-α evolves according to the equation∂Tα∂t−DTα∇2Tα=λTαMMA1+I13/KI13(1+λTαIγIγIγ+KIγ)︸production−dTαTα︸death.[9]TNF-α is secreted by activated macrophages (6, 7), a process enhanced by IFN-γ (32) and inhibited by IL-13 (33).

Equation for IL-2 (I2).

IL-2 is produced by Th1 cells (16):∂I2∂t−DI2∇2I2=λI2T1T1︸production−dI2I2︸degradation.[10]

Equations for Other Cytokines: GM-CSF (G), IL-13 (I13), IL-10 (I10), and CCL20 (C).

Activated alveolar macrophages produce GM-CSF (34), IL-10 (8), IL-13 (9), and CCL20 (10). Hence,∂G∂t−DG∇2G=λGMMA︸production−dGG︸degradation,[11]

∂I10∂t−DI10∇2I10=λI10MMA︸production−dI10I10︸degradation− dI10MAMAI10I10+KI10I12I12+KI12︸absorption,[12]∂I13∂t−DI13∇2I13=λI13+λI13MMA︸production−dI13I13︸degradation,[13]

∂C∂t−DC∇2C=λCMMA+λCT17T17︸production−dCC︸degradation− (dCRTr+dCHT17)CKC+C︸loss by chemotaxis.[14]

The absorption term in Eq. 12 is based on the fact that IL-10 enters macrophages to block production of IL-12 at the transcription level (31). The rate of absorption depends on the level of IL-12 in the microenvironment and is taken to be I12/(I12+KI12).

Eq. 14 includes a loss due to the chemotaxis by CCL20, which is bound and internalized by Treg and Th17 that are chemoattracted by CCL20.

Equations for the Velocity u.

We assume that the cells are distributed uniformly throughout the granuloma, and their total density is 0.1 g/mL (35), so thatMA+T1+T17+Tr=0.1.We also assume that all cells have approximately the same volume and surface area, so that the diffusion coefficients of all cell types are the same. By adding Eqs. 2–5, we get∇⋅u=(f+λMIγIγIγ+KIγ+λMGGG+KG+λMTαTαTα+KTα)M0+(λTI12I12KI10+I10MA+λTI2I2KI2+I2T1)11+Tr/KTr−dMAMA−dTT+λTrI2I2KI2+I2Tβ1+Tα/KTα−∇(χCR∇C)+λT17Tβ(1+Tr/KTr)(1+I2/KI2)−∇(χCH∇C)− dTrTr−dT17T17.[15]

Results

The parameter values of the system of Eqs. 1–14 are given SI Text. In this section we simulate the mathematical model developed in the previous section. For computational simplicity, we assume that the granuloma is a sphere with radius r=R(t) and all of the variables are radially symmetric [that is, functions of (r,t)], satisfying a no-flux boundary condition. The boundary is moving by the kinematic lawdR(t)dt=u(r,t),[16]where u is the component of the velocity u in the radial direction. We take the initial state to be that of a healthy individual. Some of the unknown parameters of the model were chosen so that the cytokine levels in healthy lung tissue coincide with the levels reported in ref. 22.

To validate the model, we assume that the initial stage of the disease triggers macrophage inflammatory reaction. We then simulated the development of a granuloma by taking a spherical tissue with initial radius R(0)=0.1 cm and an inflammation level f=0.4 in Eq. 2. Fig. 2 shows the profiles of all of the cells and cytokines for the first 100 d, at which time the disease reached steady state; note that the granuloma radius increased from R(0)=0.1 cm to R(100)=0.13 cm. Patient data (n = 11) of cytokine concentration in lung tissue are reported in ref. 22. We can compare our results (at day 100) with the data in ref. 22. Fig. 3 shows a good fit of our simulation results with the patient data.

Fig. 2.
Fig. 2.

Simulation results during the first 100 d since the start of the disease.

Fig. 3.
Fig. 3.

Comparison with clinical data.

The slight discrepancy in the level of CCL20 can be attributed to the fact that in ref. 22 it was chemokine MIP-1α that was measured rather than chemokine CCL20.

Treatment.

The most commonly used agents in the treatment of pulmonary sarcoidosis are corticosteroids; taken orally they provide relief of symptoms and control potentially disabling respiratory impairments (36). However, the exact mechanism of the action of the drugs is unknown, and they do not cure the disease (37). Infliximab, an anti–TNF-α drug, is used for chronic resistance sarcoidosis, but it has serious side effects and its effectiveness is uncertain (36).

From Fig. 2 we see that, starting from heathy state, the sarcoid granuloma radius will increase in 100 d from radius R(0)=0.1 cm to radius R(100)=0.13 cm. A complete recovery, in our model, will reduce the granuloma radius to 0.1 cm. Clinical tests for sarcoidosis use expression of FVC, cytokine expression in bronchoalveolar lavage (BAL), or BAL in fluid (BALF). Forced vital capacity (FVC) tests are easy to take; however, they do not provide a good representation of the state of the disease and may also be inconsistent. Data for anti–TNF-α reported in ref. 38 by FVC test show inconsistency in the effect of the drug over time, particularly as the level of the drug injection is increased. Using BAL measurements, it is reported in ref. 39 that TGF-β is a regulator of the inflammatory process in sarcoidosis. Expression of cytokine mRNA in BALF of patients with sarcoidosis and healthy controls is reported in ref. 40, including response to reduced Th1 cells.

Here we use our model to explore the efficacy of several drugs in terms of how they reduce the radius of the granuloma.

In the simulation of sarcoidosis in Fig. 3, we have taken the inflammation (in Eq. 1) to be f=0.4. By day 100 the granuloma stabilized and all of the concentrations became constant. In the sequel, starting with these constant values, treatment will begin at day 100, it will be continuous, and we expect the granuloma radius to continuously decrease until it reaches some steady-state value. We begin with the anti–TNF-α drug and assume that treatment with anti–TNF-α reduces the growth rate λTαM in Eq. 8 by a factor 1/(1+0.5). Fig. 4 shows how the radius R(t) decreases, reaching steady state around week 20.

Fig. 4.
Fig. 4.

The profile of the radius under anti–TNF-α treatment administered after week 15.

We proceed to use the model to explore other potential drugs: anti IL-12, anti IFN-γ, and injection of TGF-β.

We represented the effect of anti–IL-12 by reducing by half the production rates λI1240MA and λI1270MA in Eq. 8. Fig. 5 shows a decrease in the granuloma radius. A steady state is reached after ∼20 wk.

Fig. 5.
Fig. 5.

The profile of the radius under anti–IL-12 treatment administered after week 15.

Next we consider anti–IFN-γ and represent its effect by reducing the production rates λIγT1, λIγM by half in Eq. 6. Fig. 6 shows how radius R(t) changes under treatment.

Fig. 6.
Fig. 6.

The profile of the radius under anti–IFN-γ treatment administered after week 15.

Finally we consider injection of TGF-β and represent its effect by introducing a source term (10−11 g⋅mL−1⋅d−1) in Eq. 7. Fig. 7 shows the reduction of R(t) under TGF-β injection.

Fig. 7.
Fig. 7.

The profile of the radius under injection with TGF-β treatment administered after week 15.

We note that the steady states of R(t) in Figs. 6 and 7 are all different, and they depend on the “amount” of drug that was administered.

Discussion.

Sarcoidosis is a disease whose origin remains a mystery. Pulmonary chronic sarcoidosis is currently treated by drugs that are generally known to reduce inflammation, but not curative. Among these drugs, infliximab is perhaps the most specific, an anti–TNF-α drug. In an attempt to explore the progression of the disease, we developed in this paper a mathematical model based on patient data (22). The model is represented by a system of partial differential equations within a granuloma of varying radius R(t). The variable quantities in the model are cells (macrophages and T cells) and cytokines. We assume that the disease is associated with initial inflammation, and we then use the model to simulate the growth/decrease of each of the variables and of the radius of the granuloma. The resulting model closely approximates the sustained multicellular, multicytokine inflammatory response (for the first 100 d) following antigen stimulation, which is characteristic of the Th1 immune response causing sarcoidosis. In particular, we have simulated the effect of infliximab on the decrease in the granuloma radius R(t). We also explored the effect on R(t) of other potential treatments: anti IL-12, anti IFN-γ, and injection of TGF-β, which are viable therapeutic targets based upon clinical evidence implicating higher levels of IL-12 and IFN-γ (24, 41) and lower levels of TGF-β (39, 42) with chronic disease activity. Whereas no treatments directed specifically at these targets have been assessed in the clinical setting in terms of altering sarcoidosis disease burden (e.g., granuloma prevalence in the lung), the model suggests that these predicted therapeutic targets may reduce disease activity. Thus, the model could be used to predict treatment responses in the preclinical setting.

Our model of sarcoid granuloma was based on Fig. 1. A similar network can be used to describe granuloma in tuberculosis (TB), although in that case one has to consider both classically activated and alternatively activated macrophages (43). However, the source of inflammation in TB arises from the TB antigen (i.e., the Mycobacterium tuberculosis); hence some of the model parameters will have to be changed in the TB case, leading to different conclusions. Indeed, BAL measurements in pulmonary sarcoidosis and pulmonary TB show differences in the expression of cytokines (44, 45). Big differences may occur when some of the cytokines are overexpressed or underexpressed. For instance, a genetic variant associated with excessive IFN-γ production in response to TB antigens may predispose the lung to sustained granuloma formation (sarcoidosis) while protecting against TB (because it will kill the bacteria and hence limit the inflammation). On the other hand, if the genetic defect were associated with impaired IFN-γ production to the same antigenic challenge, this condition would favor the development of latent TB (inability to kill/clear the organism) but would protect against sarcoidosis (inflammation would be self-limited and would readily resolve).

The present work is a step toward a more comprehensive study of sarcoidosis and its treatment. As more data become available, the model could be further refined. It would be important to include in this refined model adverse side effects of drugs.

Methods

All of the computations used to solve the PDE system apply second-order finite difference discretization on the radial direction and a forward Euler method on the time direction.

To support the robustness of the simulation results, we ran sensitivity analysis on parameters that appear in the differential equations and in the boundary conditions; details are in SI Text.

The experimental results displayed in Fig. 3 were obtained by the following procedure: Gene expression analysis was performed on tissues obtained from patients with sarcoidosis at the time of diagnosis compared with normal lung tissue. Expression of select genes was further confirmed in lung tissue from a second series of patients with sarcoidosis and disease-free control subjects by semiquantitative RT-PCR. The expression of proteins corresponding to selected overexpressed genes was determined using fluorokine multiplex analysis, and immunohistochemistry.

Footnotes

  • ↵1To whom correspondence should be addressed. Email: afriedman{at}mbi.osu.edu.

  • Author contributions: W.H., E.D.C., and A.F. designed research, performed research, analyzed data, and wrote the paper.

  • Reviewers: M.J., Albany Medical College; and R.C.L., University of Connecticut Health Center.

  • The authors declare no conflict of interest.

  • This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1417789111/-/DCSupplemental.

Back to top
Article Alerts
Email Article
Citation Tools
Request Permissions
Mathematical model of sarcoidosis
Wenrui Hao, Elliott D. Crouser, Avner Friedman
Proceedings of the National Academy of Sciences Nov 2014, 111 (45) 16065-16070; DOI: 10.1073/pnas.1417789111
del.icio.us logo Digg logo Reddit logo Twitter logo CiteULike logo Facebook logo Google logo Mendeley logo

Article Classifications

Proceedings of the National Academy of Sciences: 111 (45)
Table of Contents

Submit

Jump to section

You May Also be Interested in

Similar Articles