An In Silico Modeling Approach to Understanding the Dynamics of Sarcoidosis

Background Sarcoidosis is a polygenic disease with diverse phenotypic presentations characterized by an abnormal antigen-mediated Th1 type immune response. At present, progress towards understanding sarcoidosis disease mechanisms and the development of novel treatments is limited by constraints attendant to conducting human research in a rare disease in the absence of relevant animal models. We sought to develop a computational model to enhance our understanding of the pathological mechanisms of and predict potential treatments of sarcoidosis. Methodology/Results Based upon the literature, we developed a computational model of known interactions between essential immune cells (antigen-presenting macrophages, effector and regulatory T cells) and cytokine mediators (IL-2, TNFα, IFNγ) of granulomatous inflammation during sarcoidosis. The dynamics of these interactions are described by a set of ordinary differential equations. The model predicts bistable switching behavior which is consistent with normal (self-limited) and “sarcoidosis-like” (sustained) activation of the inflammatory components of the system following a single antigen challenge. By perturbing the influence of model components using inhibitors of the cytokine mediators, distinct clinically relevant disease phenotypes were represented. Finally, the model was shown to be useful for pre-clinical testing of therapies based upon molecular targets and dose-effect relationships. Conclusions/Significance Our work illustrates a dynamic computer simulation of granulomatous inflammation scenarios that is useful for the investigation of disease mechanisms and for pre-clinical therapeutic testing. In lieu of relevant in vitro or animal surrogates, our model may provide for the screening of potential therapies for specific sarcoidosis disease phenotypes in advance of expensive clinical trials.


Introduction
Sarcoidosis is a chronic granulomatous disease of unknown cause, for which relevant research models are lacking. Human research in sarcoidosis is hindered by the existence of diverse clinical phenotypes, presumably relating to genetic and environmental variables [1]. Genetic variability may also explain the unpredictable response to treatment among sarcoidosis patients [1]. Given the genetic diversity of the disease, environmental variables (e.g., antigen exposures) and the lack of relevant animal models, it would be necessary to recruit large numbers of patients, at a substantial cost, to represent all of the sarcoidosis phenotypes using conventional clinical research approaches. Alternatively, new generation, high-throughput genetic screening platforms provide an unprecedented opportunity to stratify the molecular basis of sarcoidosis disease phenotypes with the ultimate goal of individualizing therapy. To this end, it will be necessary to determine how genetic variability influences disease pathogenesis and treatment.
In this report, we focus on sarcoidosis phenotypes that are suspected to arise from defective antigen-dependent Th1 type immune responses associated with deregulated interactions among essential immune cells such as T effector cells, T regulatory cells, and antigen-presenting macrophages. The interactions among these cells are mediated by cytokines such as IL-2, IFNc, and TNFa. We hypothesized that this complex interaction network contained sufficient information for the investigation of ''normal'' and ''sarcoidosis-like'' Th1 responses to antigens. Thus, we developed a computational model to represent the dynamics of this interaction network and its responses to perturbations. Our results are the first demonstration of an in silico model of granulomatous inflammation with potential applications for mechanistic and therapeutic research relating to sarcoidosis and other related diseases.

A minimal model for Th1 activation
The hallmark of sarcoidosis is the preponderance of Th1 immune response to poorly characterized antigens. The differentiation of naive T cells upon antigen presentation and polarizing conditions has been the subject of previous mathematical modeling (see, for example, Yates et al. [2], Mariani et al. [3], and Mendoza et al. [4]). To cope with the complexity of the immune network, we first simplified these models to a minimal model that captures Th1 activation characterized by a bistable switch, as these previous models have shown. The following equation can describe this bistable switching behavior (without explicitly considering interactions with other immune cells for now): Note that we are interpreting T as Th1 activity or the density of CD4+ T effector cells (Teffs) with the Th1 phenotype. The first and last terms on the right-hand side of this equation set the net basal flux (influx rate, b T , minus the clearance rate, e T T) of Teffs in the system. These terms give the basal steady state level of T equal to b T /e T . The rate of Th1 activation is assumed to be given by the Hill-type function in the second term of Eqn 1. This sigmoidal function represents the autocatalytic increase in Th1 activity after antigen presentation and subsequent autocrine cytokine signaling which occurs for positive values of f T . As will be shown in the next subsection, the factor f T links Teffs to the influence of antigens (of density A) and cytokines (of density c); for now we use the following equation: where g 1 and g 2 are non-negative coefficients. Supplementary information, represented by Eqn S1 and Figure S1 (see Text S1), proves that Eqn 1 gives rise to a bistable range of f T ; in other words, within this range of f T , the system can be at low or at high Teff steady state activities, and which state is achieved depends on the initial conditions of the system. This bistability predicts that as f T is increased from a low value, there is a threshold value of f T where there occurs a discontinuous switch to high Th1 activity. In other words, sufficient increase in antigen and/or cytokine densities can bring about a sharp transition to high Teff steady state (T s ) levels. Alternatively, the system can switch to the upper branch of steady states within the bistable range by perturbing T so that its value crosses the middle branch of the curve. This is illustrated in Figure 1 where a square pulse of A is applied to the system, and a pulse of sufficient amplitude succeeds in switching the system to a larger T s .
Regulatory T cells (Tregs) are expected to increase the threshold for achieving high steady state T activity. As per Eqn 1, an increase in the value of the parameter h T in this simple model corresponds to an increase in Treg density in the detailed model discussed in the next subsection. This relationship is schematically represented in Figure S1 (see Text S1), wherein the effect of greater Treg activity (h T ) is to increase the value of f T where the switch to high T s activity occurs -that is, the antigen or cytokine threshold for promoting high-level Teff activity is increased.
A model involving essential components of granulomas (Teffs, Tregs, and macrophages) Here, we present a network model that involves the interactions of Teffs with macrophages and Tregs. The predisposition to sarcoidosis has implicated some dysfunction of naturally occurring, innate, Tregs [4]. Macrophages are also essential in the model because, first, they are known to be highly recruited by active Teffs to become part of granulomas (they represent the bulk of the granuloma by mass), macrophages synergize with Teffs in the secretion of IFNc and TNFa, and macrophages interact with and affect Treg activity [5][6][7][8][9]. The model immune network shown in Figure 2A minimally depicts the essential interactions among Teffs, Tregs, macrophages, and key cytokines secreted by these cells that mediate cell-cell interactions, namely, IL-2, IFNc and TNFa. As in the previous section, this model is phenomenological in that it does not include all known mechanistic details but rather captures the qualitative dynamics of the interactions among the players of the system -that is, how one player promotes or inhibits the densities or activities of other players in the network. Figure 2B summarizes the net interactions between cells in a so-called qualitative network (qNET), as we described previously [10]. As will be demonstrated in the next section, this qNET is useful to guide our intuition on the potential dynamics of the detailed model of Figure 2A.
For the model in Figure 2A, the rate of change of T is described by: where f T~k1a AMzk 1b dzk 1c c A, M, d and c are densities of antigens, macrophages, IL-2 and IFNc, respectively. Note that the product AM represents antigen presentation to Teffs by macrophages. We assume that macrophages are the only antigen-presenting cells, and increasing M leads to increased antigen presentation to Teffs. There are several mechanisms by which Tregs antagonize the activity of Teffs [9,11], but we subsume all these mechanisms in the inhibitory effect of Tregs on the proliferation of Teffs as implemented by including k 1d R in the denominator of the second term on the right-hand side of Eqn 3. Thus, the downregulation by Tregs of the production of IL-2, IFNc and TNFa is indirectly caused by the inhibition of Teff proliferation.
The following equations for the rates of change in the densities of Tregs (R) and macrophages (M) can be understood from the interactions shown in Figure 2: where f M~k6 czk 6c a Also, looking at Figure 2, the rates of change in the densities of the cytokines can be described by the following differential equations: Equations (3)-(8) form the complete set of ordinary differential equations that describe the dynamics of our model (see Figure S2).

Potential dynamics from the qNET structure
The qNET shown in Figure 2B summarizes the qualitative interactions between pairs of immune cells. This qNET can already give us some intuition as to what to expect for the longterm states of the system. Arrows mean ''activate'' or upregulate, while hammerheads mean ''inhibit'' or downregulate. Without Tregs, it is obvious that Teffs and macrophages synergize and activate each other, but only to a point where high levels of TNFa (a by-product of this synergy) are sufficient to suppress Teffs [this regulatory activity of macrophages towards Teffs has been suggested by Nagar et al [12] (see below)]. Without this macrophage-dependent suppression of Teffs, and without Tregs, a runaway immune response (unchecked activation and proliferation of Teffs and macrophages) is predicted.
Tregs exert suppressive activities on both Teffs and macrophages, as depicted through edges b and f, respectively. The solid black circle at the end of edge a is initially an arrow (activation) representing the IL-2-dependent activation of Tregs by Teffs [13]. However this arrow may turn into a hammerhead (suppression) if the negative effects of TNFa (whose levels are also induced by Teffs) begin to dominate [12]. With edge a being an arrow, the pair of edges a and b is a negative feedback loop that could lead to an oscillation between the numbers of Teffs and Tregs (this oscillation could be damped or sustained, as previous mathematical modeling, and experimental observations, have shown [14,15]. If edge a is a hammerhead, a toggle switch (mutual inhibition) could ensue; in other words, either Tregs annihilate Teffs or vice versa. Note that the interactions between Tregs and macrophages are also assumed to constitute a toggle switch (edges e and f), based upon established mechanisms [16,17]. . T in Eqn 1 is interpreted as Th1 activity. A square pulse of antigen with amplitude A1 = 3.5 (applied from t = 5 to t = 10) fails to increase T, while increasing the amplitude to A2 = 4 leads to a rapid increase in Th1 activity at the level T2. Parameters: hT = 1, bT = 0.02, eT = 1, c = 3, g 1 = g 2 = 1, initial T = 0.34. doi:10.1371/journal.pone.0019544.g001 What are the possible long-term state levels of Teffs, Tregs and macrophages as predicted by the qNET of Figure 2B? Using only a rough characterization of cellular levels (low and high, symbolized by 0 and +, respectively), not all of the 8 possible states can be realized due to the structure of the qNET, as shown in Table 1.
The 'healthy state'' (state 1) is characterized by very low or nonexistent levels of Teffs, Tregs and macrophages. Possible ''sarcoid states'' are states 5-8, although state 5 is predicted to be unlikely from the qNET structure alone (because macrophages go up with Teffs in the absence of Tregs which is the only source of inhibition according to the qNET). Note that states 6 and 7 involve simultaneous high activities of Th1 and Tregs which, as will be explained in the next section, explains the immune paradox discussed by Miyara et al [4].
The important question that modeling hopes to contribute is to predict the possible scenarios on how a healthy state switches to a sarcoid state. For this, we need to perform computer dynamic simulations of the full model. The base parameters used for the detailed model (depicted in Figure 2A) are provided in Table 2.

Thresholds of sarcoidosis (Th1 activation) related to amplitude of antigen exposure
The typical immune response to a single antigen exposure is self-limited. However, antigen exposures of sufficient duration or amplitude are capable of inducing a new steady state characterized by dramatic expansion of immune cell populations along with sustained increases in pro-inflammatory cytokines (Figure 3). This scenario would correspond to persistent granulomatous inflammation (active sarcoidosis).
Switching behavior of Treg levels due to sensitivity to IL-2 Tregs are considered to be important determinants of granulomatous inflammation in patients with sarcoidosis [5,18]. Based upon our modeling assumptions ( Figure 2A) and previous reports [15,19], it would appear that Treg activity is critically dependent upon IL-2. Figure 4 confirms that steady state level of Tregs sensitively depends on the parameter k 2a (which is the coefficient of the dependence of Treg proliferation on IL-2). Figure 4B corresponds to state 8 in Table 1, wherein Tregs are elevated along with Teff and macrophages. This represents the clinical scenario wherein active pulmonary sarcoidosis is paradoxically associated with anergy to new antigens (i.e., suppression by Tregs) [4].

Modeling therapeutic interventions
Assuming no change in Treg sensitivity, as modeled in Figure 4, we sought to determine the effects of suppression of either IL-2 or IFNc levels (corresponding to anti-IL-2 or anti-IFNc therapy). As shown in Figure 5, ''monotherapy'' with either anti-IL-2 or anti-IFNc is capable of suppressing Treg and macrophage activation, respectively, without affecting high steady-state activation of Teff cells. In contrast, Figure 6 shows that combinations of anti-IL-2 and anti-IFNc at appropriate ''doses'' (as represented by the amplitude of inhibition) maintains all cell lines in the low steady-state activation state following antigen  Table 1. Long-term states predicted by the qNET shown in Figure 2B. challenge (i.e., the normal response). Thus, assuming that the modeling conditions accurately depict the mechanisms regulating granuloma formation in the setting of sarcoidosis, careful titration of combination therapy with anti-IL-2 and anti-IFNc is predicted to effectively attenuate disease activity.

Anti-TNFa therapy can promote the sarcoidosis phenotype
Considering that sarcoidosis is a polygenic disease with diverse phenotypes, models must be adjusted to represent each disease phenotype. Whereas anti-TNFa therapy often suppresses sarcoid- Table 2. Base parameter values used for the detailed model in Fig. 2A (see also Fig. S2).   osis disease activity [20,21], one clinical variant of sarcoidosis is actually promoted by anti-TNFa therapy [22,23]. In this regard, it is interesting that the current model assumptions predict that monotherapy with anti-TNFa would either be ineffective or could paradoxically promote the sarcoidosis state (Figure 7).

Discussion
In silico modeling of lung disease is in its infancy, as reflected by the few published attempts to date, particularly relating to lung physiology [24,25]. However, leading scientific organizations like the National Institutes of Health (NIH) anticipate the need for in silico modeling to accommodate the exponential growth of information emerging from human genetic studies. With respect to strategic planning towards the goal of providing ''personalized medicine'' based upon genetic profiling, the NIH seeks to ''use in silico and statistical modeling/simulation strategies to design and inform genetic/pharmacogenetics clinical trials, so results can eventually be used for tailoring treatment or prevention'' (http:// www.nhlbi.nih.gov/about/dcvd/sp/dcvd-sp-goal-1.1.htm). The in silico sarcoidosis model presented here possesses relevant features, including representation of ''normal'' and ''disease'' phenotypes, and the capacity to perform preclinical therapeutic testing. As such, this model serves as a promising template for future sarcoidosis research.
Based upon putative interactions among critical components essential for granuloma formation, our model demonstrates that a biphasic immune response is expected depending on the The system is first allowed to reach the low steady states of all variables, and then exposed to a square pulse of antigen (with amplitude A = 100) from t = 20 to t = 70 as indicated by the arrows on the abscissa. All parameters as in Table 2, except k 2a which when increased above a specific threshold (compare Panel A vs Panel B) promotes a dramatic increase in Tregs. doi:10.1371/journal.pone.0019544.g004 nature of the antigenic challenge. As shown in Figure S1, suppression of Treg function, as was demonstrated in the lungs of sarcoidosis patients [4] presumably in response to TNFa [9], is sufficient to reduce the threshold for a sustained antigendependent Th1 response. In lieu of an altered immune response (i.e., normal host), the model further predicts that a sarcoidosislike Th1 response could be elicited by an antigen challenge of sufficient intensity. This model behavior is consistent with a number of clinical reports of previously healthy patients developing sarcoidosis following a single, massive antigen exposure, such as the case in first-responders to the 9/11 disaster in New York City [26,27].
To the extent that the computational model represents mechanisms predisposing to unregulated granuloma formation in the context of sarcoidosis, various ''treatments'' can be simulated by mathematically manipulating one or more of the extracellular immunomodulating molecules. This approach is illustrated in Figures 5, 6 and 7, wherein suppression of any single molecule (IL-2, IFNc or TNFa) fails to decrease Teff density in the sarcoidosis phenotype, whereas simultaneous suppression of IL-2 and IFNc restores the system to the baseline state (i.e., disease remission). These results are strongly dependent on the modeling assumptions, including the functional status of specific components of the network, as represented by changes in the kinetics or the character (inhibition, activation) of the interactions of one or more of the components in Eqns 3-8. A key feature of our model is that it conveniently allows for adjustments to accommodate diverse disease mechanisms. As a specific example, the sarcoidosisassociated BTNL2 mutation which results in a truncation mutation of the BTNL2 protein, a B7-family protein predicted to inhibit interactions between antigen-presenting cells and Teffs [28], would be represented by an increase in f T in Eqn 3. This would have the same effect as diminishing Treg function (h T , Eqn 3) in that the threshold for achieving the high Ts would be reduced, favoring the sarcoidosis phenotype. Figure 5. Single therapies using anti-IL2 (Id) or anti-IFNg (Ic). The system is first allowed to reach low steady states of all variables, and then exposed to a square pulse of antigen (with amplitude A = 100) from t = 20 to t = 70 to obtain a sarcoid state. As shown in (A) and (B), the model predicts that there is a threshold of anti-IL2 (between amplitudes Id = 29 and Id = 30) that will push the steady state of R to a low level, but not T and M. The anti-IL2 therapy is represented by asquare pulse of Id given between t = 100 and t = 120. And as shown in (C) and (D), the model also predicts that there is a threshold of anti-IFNc (between amplitudes Ic = 15 and Ic = 16) that will push the steady state of M to a low level, but not T and R. The anti-IFNc therapy is represented by a square pulse of Ic given between t = 100 and t = 150. All parameters as in Table 2 Interestingly, our model suggests a resolution of the clinical paradox discussed by Miyara et al [4], wherein excessive antigenmediated inflammation in the lungs coexists with impaired antigen responsiveness (anergy) elsewhere. Recent studies suggest that expansion of Tregs at the sites of active granulomatous tissue inflammation causes suppression of T cell proliferation elsewhere, leading to peripheral lymphopenia and anergy, which are welldocumented manifestations of active pulmonary sarcoidosis [29,30]. This scenario is represented by conditions 6 and 8 in Table 1, and Figures 4B, 5C, 6D, 6F, 7C. Admittedly, our model is incomplete in that detailed information relating to the actual kinetics and functions of the system Figure 6. Combination therapies using anti-IL2 and anti-IFNc. The system is first allowed to reach low steady states of all variables, and then exposed to a square pulse of antigen (with amplitude A = 100) from t = 20 to t = 70 to obtain a sarcoid state. For (A)-(C), k 2a = 2.0. Combinations of anti-IFNc (Ic) and anti-IL2 (Id) are applied as square pulses at the same time, between t = 100 and t = 120. The combination in (A) is not successful, while those of (B) and (C) are successful in pushing down the levels of T and M (although there is a slight increase in R compared to its level in (A), R is still very low). For (D)-(F), k 2a = 2.1. An increase in anti-IL2 (Id) from 29 to 40 brings down R, but not T and M. Note that a small addition of Ic as shown in (F) returns R to the high steady state level. All parameters as in Table 2, except k 2a . doi:10.1371/journal.pone.0019544.g006 Figure 7. Single therapy using anti-TNFa (Ia). The system is first allowed to reach low steady states of all variables, and then exposed to a square pulse of antigen with amplitude A = 100 from t = 20 to t = 70 (indicated by the first two black arrows from the left) to obtain a sarcoid state. (A) Without anti-TNFa (Ia = 0). (B) Anti-TNFa is introduced as a square pulse with amplitude Ia = 21.1 from t = 100 to t = 140 (indicated by the last pair of thick gray arrows); (C) Anti-TNFa is introduced as a square pulse with amplitude Ia = 21.2 from t = 100 to t = 140 (indicated by the last pair of thick gray arrows). All parameters as in Table 2 components are not well characterized, and many ancillary cells and molecules are not considered. For instance, recruitment of macrophages induced by IL-2 is not included in model [31] and many cells and cytokines involved in granuloma formation (e.g., natural killer T cells, a/d T cells, IL-12, IL-23 [8,32]) are not considered. Nonetheless, the qualitative framework of the model is adaptable as more details relating to mechanisms of granuloma formation and disease pathogenesis emerge, including the effects of genetic variability. Moreover, model simulations are convenient in that they can be run using inexpensive and readily-available computational facilities, as well as practical in that they accurately portray relevant clinical phenotypes, even in the simplified form presented herein.
In summary, we provide the first dynamic computer simulation of granuloma formation in sarcoidosis, a disease for which relevant in vitro or animal models do not currently exist. The qualitative network is shown to accurately model distinct disease mechanisms, including a sarcoidosis state resulting from either massive antigen exposure or altered immune (e.g., Treg) function, and specific sarcoidosis phenotypes (e.g., lymphopenic sarcoidosis). The model demonstrates the capacity for pre-clinical therapeutic testing, and adaptability based upon distinct disease immunophenotypes, such as those dictated by genetic variability.

Methods
To determine the steady states of the detailed model, the baseline values for each parameter was set to the values shown in Table 2 and further modified, as described for each experimental condition, and the corresponding systems of differential equations (Eqns 3-8, above) were solved numerically using the BerkeleyMadonna software (http://berkeleymadonna.com). Text S1 Supporting information demonstrating a bistable switch between low and high Teff steady states predicted by the model. (DOC)