Inflammatory Bowel Disease: How Effective Is TNF-α Suppression?

Crohn’s Disease (CD) results from inappropriate response toward commensal flora. Earlier studies described CD as a Th1 mediated disease. Current models view both phenotypes as a continuum of various permutations between Th1, Th2 and Th17 pathways compounded by a range of Treg disfunctions. In the present paper, we develop a mathematical model, by a system of differential equations, which describe the dynamic relations among these T cells and their cytokines. The model identities four groups of CD patients according to up/down regulation of Th1 and Th2. The model simulations show that immunosuppression by TNF-α blockage benefits the group with Th1High/Th2Low while, by contrast, the group with Th1Low/Th2High will benefit from immune activation.

Patients with inflammatory bowel diseases have elevated levels of circulating and gut mucosal cytokines [1]. Downstream signaling from these inflammatory mediators, through Janus kinase (JAK) and signal transducers and activators of transcription (STAT) proteins, activate transcription factors T-bet, GATA3, Foxp3 and RORγt [3]. A complex network of regulatory feedback loops involving these cytokines and their targets, are responsible for the polarization of naïve T cells into specific T helper cells: Th1, Th2, Th17 and Treg [3,4].
Earlier animal and human studies described CD as a Th1 mediated [5]. Current models view both phenotypes as a continuum of various permutations between the Th1, Th2, and Th17 pathways compounded by a range of Treg dysfunctions [4].
The development of current biological therapies in CD, mirrors our understanding of immune regulatory pathways. Unfortunately the clinical triage of CD patients, based on Montreal classification, cannot identify relevant immune targets in a given patient [6]. Thus our treatments are inherently a trial and error approach. Pretreatment knowledge of the relevant gut mucosal immune dysfunction in a given patient would significantly improve the risk/benefit balance and open the way toward personalized medical care. We have previously developed a theoretical model that described the relationship between Th1, Th2 and Treg circuits in patients with CD [7]. Our current study expands this model to include the Th17 pathway which plays an important role in inhibiting Treg cell differentiation that associated with autoimmune disorders and inflammation [8,9]. The levels of Th17 cell transcription factor and related cytokines have been collected in our clinical data to provide better predictions of disease outcome through our new model. Furthermore, based on patients' data, it allowed us to simulate the effect of TNF-α suppression in a cohort of patients with CD, and thus identify distinct groups which will benefit from TNF-α suppression and groups which may benefit from immune activation rather than immune suppression.

Mathematical model
A system of differential equations was developed based on the network that incorporates cytokines and transcriptions factors relevant to Th1, Th2, Th17 and regulatory T cells pathway as shown in Fig 1. The relative concentrations of either immune cells or inflammatory mediators were defined in g/cm 3 and were based on a theoretical density in a cm 3 of tissue. The variables (concentrations) included in the model are listed in Table 1.
Differential equations for cytokines. IFN-γ is produced by Th1 cells and activates M1 macrophages [10]. In order to derive an equation which expresses the dynamics of I γ , we denote by dI γ the change in I γ over a small time interval dt, and compute the quotient dI γ /dt (the derivative of I γ with respect to t, if dt is infinitesimally small). If I γ is produced by Th1 cells and the production rate coefficient for IFN-γ by Th1 cells is ν γ1 , then, at a given Th1 cell density T 1 , we can write dI γ /dt = ν γ1 T 1 . Similarly, the production of this cytokine by activated macrophages would be described by the differential equation dI γ /dt = ν γM M 1 . We finally have to account for the fact that IFN-γ tissue concentration decays at some known rate δ γ . Therefore the final formula for the dynamics of I γ is given by the differential equation: Similar differential equations were established for IL-2, IL-4, IL-6, IL-21, TGF-β and TNF-α. IL-2 is secreted by Th1 cells [11,12], so we have Th2 cells and M2 macrophages produce IL-4 during inflammation [4,10,13]; hence IL-21 is produced by Th17 cells [14,15], so that Inflammatory Bowel Disease: TNF-α Blockage IL-6 and TNF-α are produced by M1 [14,15] and Th1 cells also produce TNF-α [16]; hence, the equations of I 6 and I α are TGF-β is mainly expressed from M2 macrophages and Treg cells [4,10,17]; hence we have Cytokine IL-10 can be the product of both M2 macrophages and Treg cells [4,10,17], and IL-2 further increases IL-10 production by Treg cells [18]. IL-2 binds its cognate receptors and it is rapidly internalized. Allowing for the receptor recycling rate (n 2r ) and maximum amount of IL-2 bound at any given time, we can write an equation that links this rate to the tissue concentration of IL-2 (I 2 ) and Treg (T r ): n 2r (I 2 /(z 2 + I 2 )) T r where z 2 is a constant. Therefore the differential equation that incorporates IL-10 secretion from macrophages, Treg cells and IL-2 contribution can be written as: IL-12 is produced by M1 macrophages, and IL-10 can suppress its expression [19]. The IL-10 mediated suppression can be described using the expression (1+ I 10 /z 10 ) where z 10 is a constant specific to IL-10 activation rate. Thus the differential equation for the net IL-12 production becomes: Differential equations for macrophages. The specific polarization of naïve T cells toward Th1, Th2, Th17 and Treg phenotypes is supported by cytokines produced by activated macrophages [4]. Activated macrophages fall into two major groups, M1 and M2. Increased tissue concentration of TNF-α in IBD patients can lead undifferentiated M0 macrophages to polarize to M1 macrophages, classically activated macrophages [20]. M1 macrophages can also be derived from M2 in the presence of high levels of IFN-γ, a proinflammatory cytokine [21]. In contrast, M2 (alternatively activated macrophages) are induced by anti-inflammatory cytokines such as IL-10 and TGF-β [21]. The following factors were considered when obtaining a differential equation for M1 macrophage phenotype formation: concentration of tissue factors (f 1 ) that contribute to macrophage polarization; tissue concentrations of TNF-α (I α ), IFN-γ (I γ ) and TGF-β (I β ); cytokine specific rates of macrophage activation (σ Mα, σ Mβ, σ Mγ ) and transition between M1 and M2 phenotypes (σ Mβ , σ Mγ ), and macrophage apoptosis (decay) rate (μ M ). The differential equation for M1 is given by Similarly, undifferentiated M0 macrophages polarize to M2 macrophages under the influence of IL-10 [22], while M1 convert to M2 macrophages by TGF-β [21], so that the complete formula for dM 2 /dt, after including death rate (μ M ) is: Differential equations for T cells. We finally turn to the dynamics of the T cells. The initiation of the inflammatory process involves macrophage activation. Macrophages produce several types of cytokines including IL-12, IFN-γ, IL-4, IL-6, TGF-β, IL-10 and TNF-α, and all these cytokines activate T cells upon binding to their specific receptors. The inflammatory network in IBD involves several auto-activation loops among the four types of T cells: Th1 cells: As shown in Fig 1(A), the loop of Th1 auto-activation includes IL-12, IFN-γ, and macrophages. IL-12 induces Th1 activation while IFN-γ, produced by Th1 cells, activates macrophages to produce more IL-12 [10,13].
Th2 cells: The loop of Th2 auto-activation includes IL-4 which leads to Th2 activation while Th2 cells further drive the expression of IL-4 [4,10,13]. This positive feedback loop is shown in Fig 1(B).
Macrophages activation after encounter of microbial antigen induces Th1 development through direct contact and IL-12 production [10,13]. Inhibition of IL-12 production by macrophages may explain the ability of IL-10 to suppress Th1 development [27]. Mathematically this is expressed by the following equation: IL-2 secretion by Th1 cells can further enhance their activation [28]. On the other hand, negative regulatory signals are provided by Th2 and Treg cells [24,25]. Finally the net rate of Th1 cell production needs to incorporate the activation-induced apoptosis. Hence our original differential equation becomes: Similarly, for Th2 we accounted for IL-4 signaling, reciprocal Th2 inhibition by Th1 and Treg programs as well as activation-induced death of Th2 cells [4,10,13,24]: IL-6 in conjunction with TGF-β activates Th17 cells [8,14,15]. IL-21 is also a potent inducer of Th17 differentiation [8,14,15]; downregulation of IL-21 expression decreases Th17 cell infiltration in intestinal mucosa of IBD patients. Th17 activation is resisted by both IFN-γ and IL-4 [29]. Since Treg cells are also an inhibitor of Th17 activation [9], the final equation for T 17 is: Treg cells number is regulated by a dynamic homeostatic process that balances high rates of cell division with apoptosis. These cells have high affinity for IL-2 which is a potent negative regulator of pro-apoptotic signals in Treg cells [23]. Interaction with innate immune cells like macrophages as well as cytokines like TGF-β and IL-10 promotes survival and immunosuppressive activity of this cell population. IBD patients have increased gut production of the proinflammatory cytokine TNF-α which was shown to increase Treg cell apoptosis, and TNFα blockade can reverse this process, which is also observed in rheumatoid arthritis (RA) patients [30]. Proinflammatory cytokines that drive a Th17 response negatively regulate Treg cell development. Increased Th17 cell density in the gut mucosa of IBD patients corresponds to a reciprocal impairment in the Treg population. Taking in consideration all these positive and negative regulators of Treg cells, we derive the following equation: Parameter estimation Macrophages. We assumed that under normal conditions majority of gut macrophages display an M2 like phenotype. The parameters σ Mγ , σ M10 , σ Mβ and σ Mα represent the macrophage activation related coefficients under the influence of IFN-γ, IL-10, TGF-β and TNF-α cytokines, respectively. Given the assumption that under normal conditions M2 concentration is usually higher than M1, the activation rate of M2 is higher than that of M1; we arbitrarily take σ M10 = 10σ Mα . We also considered that the transition rates between M1 and M2 phenotypes are the same as the activation rate of M2, and take σ Mβ = σ M10 and σ Mγ = σ M10 .
Cytokines. Th1 cells in the inflamed gut likely produce much more IFN-γ than macrophages. Thus we empirically considered the IFN-γ production rate in Th1 cells compared to that of macrophages, to be ν γ1 = 5ν γM . We also assumed a similar production rate of IL-2 (ν 21 ) and IFN-γ (ν γ1 ) by Th1 cells, so that ν 21 = ν γ1 .
Mouse models of colitis have indicated that T cell derived IL-10 plays a more important role than the innate immune pool. We assumed that Treg cells produce more IL-10 than macrophages, and we took the IL-10 production rates to be ν 10M = 3.72 × 10 −4 week −1 [29,33] and ν 10r = 3ν 10M .
T cells. Maintenance of Th1 cells pool would require the activation and decay coefficients to be of similar magnitude. The coefficient σ 2 (representing the IL-2 maximal signal output) was estimated to be 1.23 week −1 while μ 1 (Th1 degradation rate coefficient) was 1.4 week -1 .
According to [34], we assume that Th17 activation rates by IL-6 is the same as that by IL-21, so that σ 6 = σ 21 . Similarly, IL-10 and TGF-β contributions to Treg activation/survival rate were taken to be of equal magnitude, so that σ 10 = σ β .
From the data described in the section Data collection and analysis we obtained, in particular, the mRNA concentrations of TNF-α, IL-6 and IL-10 for healthy individuals as well as the master regulators of T cells for healthy individuals. The cytokine concentrations of TNF-α, IL-6 and IL-10 are assumed to be proportional to those of the mRNA concentration, with proportionality parameter λ c . In [29], the cytokine concentrations are within 10 −5 -10 -9 g/cm 3 ; we take the IL-6 concentration in healthy tissue to be 8.00 × 10 −6 g/cm 3 . We can then compute λ c and then also the concentrations of TNF-α, IL-10 in healthy tissue. We accordingly get the steadystate concentration of TNF-α to be 9.75 × 10 −6 g/cm 3 and of IL-10 to be 1.54 × 10 −6 g/cm 3 .
The densities of Th1, Th2, Th17 and Treg in healthy tissue are assumed to be proportional to the concentrations of their master regulators with proportionality parameter λ T which needs to be determined; the master regulators are obtained from the mRNA analysis described in the section Data collection and analysis. In addition to λ T , there are still six important parameters that need to be determined, namely, v 6M , v α1 , σ 12 , σ 4 , σ 21 and σ β We determine these 7 parameters as follows: We assume that each half-saturation constant is the same as the steady state and solve the 15 steady-state equations of the system (1)-(15) for the following 15 variables: the 7 parameters defined above, and the 8 steady-state concentrations of IL-2, IL-4, IL-12, IL-21, TGF-β, IFN-γ, M1 macrophages and M2 macrophages.
Solving the steady-state system of 15 algebraic equations, we get the values of the parameters v 6M , v α1 , σ 12 , σ 4 , σ 21 and σ β (as shown in Tables 2 and 3 under "estimated from data"), and the steady-state concentrations of IL-2, IL-4, IL-12, IL-21, TGF-β, IFN-γ, M1 macrophages and M2 macrophages (shown in Table 4). Here we use the fact that the steady-state concentrations of λ T Th1, λ T Th2, λ T Th17 and λ T Treg are known from the analysis of the section Data collection and analysis. The half-saturation parameters are listed at the end of Table 3. Table 4 all the steady-state concentrations of cytokines, macrophages and T cells; this list partially overlaps with Table 3.

Parameter sensitivity analysis
Since the model is highly complicated with a lot of parameters and some parameters are estimated roughly from data, we performed sensitivity analysis to determine the robustness of the simulation results and effect of the parameters on the concentrations of Th1 and Th2 cells which are used to determine the types of IBD patients. In our parameter analysis, we focused on nine parameters, σ Mα , σ M10 , σ 12 , σ 2 , σ 4 σ 21 , σ 6 , σ β and σ 10 , which are the activation rates of macrophages and T cells by different types of cytokines and are the most significant in the disease dynamics. We applied the method of Partial Rank Correlation Coefficient (PRCC) [44] for our sensitivity analysis. We ran 5000 simulations in which all the nine parameters are varied according to Latin hypercube sampling with the range of ±20% perturbation around the parameter values obtained for healthy individuals. The results of the sensitivity analysis are summarized in Table 5. Fig 2 shows the scatter plots of rank transformed T 1 (Fig 2A) and T 2 (Fig 2B) after 100 weeks (solution close to steady state) versus the rank transformed parameters with significant correlation (|PRCC|>0.5) and p-value (p<0.01); the title of each subplot shows its PRCC value.
All the nine parameters have significant PRCC values either with T 1 or with T 2 . Among the nine parameters, we find that the statistically significant PRCC values with both T 1 and T 2 are σ 12 , σ β and σ 10 . The parameters σ β and σ 10 are negatively correlated to T 1 and T 2 : the two Table 3

Parameter Definition
Value References parameters are the activation rates of Treg cells which inhibit Th1 and Th2 activations. This result is consistant with the observation that the abnormal levels of T cells usually arise from abnormal regulation of Th1 and Th2 cells by Treg cells [7]. On the other hand, σ 12 is positively correlated to T 1 and T 2 : as the parameter increases, the amount of Th1 cells increases and thus Treg cell regulation on Th1 and Th2 cells decreases through increased TNF-α inhibtion produced by Th1 cells. Note that although Th1 and Th2 cells inhibit each other, reduced Treg inhibtion has a larger effect than the mutual inihbition between Th1 and Th2 cells. The roles of these three parameters (σ 12 , σ β and σ 10 ) on different types of IBD will be discussed in the Result section later.

Data collection and analysis
Human subjects and tissue biopsy collection. Colon biopsies from healthy controls and patients with Inflammatory Bowel Diseases were obtained from a de-identified tissue bank. The tissue bank is Host Modifiers of Inflammatory Bowel Diseases Phenotypes and IRB protocol number is 06-0270-F1V. The Institutional Review Board at the University of Kentucky approved the tissue bank protocol. Banked, normal biopsies were obtained from healthy controls that underwent screening colonoscopy. IBD patients had an established diagnosis of Crohn's disease involving the colon and were off immunosuppressant treatment at the time of the endoscopic procedure. Biopsies were obtained inflamed but not ulcerated colonic mucosa. All measures were derived from mucosal biopsies. No blood samples were utilized in our experiment. Patient characteristics are shown in Table 6. Cytokine profile data tends to be affected by the contamination of epithelial cells that cannot produce typical Th1/Th2 cytokines. Nevertheless all biopsies are expected to have a similar content of epithelial cells: we used the same standard biopsies forceps, applied en-face orientation, and for all IBD patients the biopsy was performed at the area abutting ulceration. If anything, biopsies from IBD patients would have been expected to contain less epithelial cells, due to surface erosion. Since control and IBD biopsies contained epithelial cells the differences in immune markers are expected to be representative of the inflammatory process. mRNA analysis. Frozen biopsies were placed in 600μl lysis buffer (MagNA Pure Compact RNA Isolation Kit) at room temperature for 10 min. Biopsies were disrupted using a MagNA Lyser instrument (Roche, Basel Switzerland) and beads: 1 st spin @ 7000x for 40 sec and 2 nd spin @6000x for 30 sec. Tubes were placed in the instrument cooling block for 15 min and then centrifuged briefly to pellet the debris. The supernatant containing total RNA was then purified by MagNA Pure Compact RNA Isolation Kit (Roche) protocol with elution volume of 50ul. For control RNA quality and quantity were determined with a spectrophotometer (NanoDrop Technologices, Wilmington, DE).
Gene expression analysis. Total mRNA (15μl) were analyzed using the nCounter Master Kit from NanoString Technologies (NanoString Technologies, Inc Seattle, WA). Gene Expression CodeSets of our interest for the nCounter Analysis System were preordered from Nano-String Technologies, Inc (www.nanostring.com). The nCounter Analysis System is an integrated system comprised of a fully-automated Prep Station, a Digital Analyzer, and the CodeSet (molecular barcodes) reader.

Results
In the sequel we shall use the following abbreviations: We obtained gut mucosal samples from healthy individuals and IBD-Crohn's patients and measured the mRNA expression for cytokines IL-6, IL-10, TNF-α and transcription factors Tbet, GATA3, RORγt and Foxp3. Based on the observed values we were able to divide the patients into four categories: For example, if an IBD patient had a T-bet level higher than T 1ssh and a GATA-3 level lower than T 2ssh , this patient is classified to Type 1 group. There are totally 58 patients: 7 in Type 1, 18 in Type 2, 17 in Type 3 and 16 in Type 4. The deviations of mRNA fold expression of IBD patients from healthy individuals are shown in Table 7.
Our next step was to determine the values of coefficients that define the production rate of individual cytokines and activation of master regulators. These theoretical parameter variations (cytokine production, T cell activation), in each of the four types, were determined using the steady-state solution of our mathematical model and are shown in Table 7. We assumed that the blood sample data reflect the tissue data (by the same proportionality parameter), and that the disease state occurred at the steady state when some of the production/activation rates related to T cells were deregulated, either increased or decreased. Here we consider the variations of the following parameters, ν g1 , ν α1 , ν βr , ν 10r , σ 12 , σ 4 , σ 21 , σ 6 , σ β and σ 10 ( Table 8).
In type 1 group the production rate of Th1 cytokines, IFN-γ and TNF-α , was decreased while the anti-inflammatory IL-10 was increased. Nevertheless the signal strength of IL-12, IL-21 and IL-6 favored activation of Th1 and Th17 pathways. On the other hand, decreased activation of Treg by TGF-β and IL-10 was able to offset their increased production rate.
In type 2 group, despite an increase in TNF-α production rate and IL-12 signal strength, the lack of IL-21 and IL-6 signal strength coupled with increased IL-4 signal strength tipped the balance in favor of Th2 response.
Type 3 group was the only one with a robust IFN-γ production rate. Interestingly we no longer observed the reciprocity between the Th1, Th2, Th17 and Treg activation rates, as the signal strength for IL-12, IL-4, IL-21, IL-6, TGF-β and IL-10 were increased across the board.
Type 4 group was characterized by a global decrease in signal strength and thus lack of T cell activation even though this group had the highest production rate of TNF-α Table 7. Fold changes of the cytokine and T cell concentrations obtained from the clinical data and the simulations in different types of diseases.

Conclusions and Discussion
Inflammatory Bowel Diseases are characterized by a deregulated immune response toward gut microbiota. A simplistic characterization of disease phenotype divides patients into either Crohn's Disease or Ulcerative Colitis. Nevertheless, if we take into account the individual variations in genetic background, microbiome, environmental exposure, nutrition, and life-style, we can foresee multiple immune phenotypes. The development of current biological therapies in CD mirrors our understanding of immune regulatory pathways. Current Montreal classification distinguishes phenotypes based on anatomical and clinical criteria. Unfortunately this will not help identify individual immune phenotype relevant to a particular biologic treatment. Pretreatment knowledge of the relevant gut mucosal immune dysfunction in a given patient would significantly improve the risk/benefit balance and open the way toward personalized medical care.
Patients with inflammatory bowel diseases have elevated levels of circulating and gut mucosal cytokines [1]. Downstream signaling from these inflammatory mediators, activate transcription factors T-bet, GATA3, Foxp3 and RORγt [3]. A complex network of regulatory feedback loops involving these cytokines and their targets, are responsible for the polarization of naïve T cells into specific T helper cells: Th1, Th2, Th17 and Treg [3,4].
In this study we developed a mathematical model of differential equations to describe immune regulatory pathways in patients with CD. The output of these differential equations was dependent on the relative mRNA expression of specific immune system targets.
Patient stratification based on calculated Th1 and Th2 immune cell activation, relative to healthy controls, identified 4 distinct immune phenotypes. These were further characterized based on cytokine production as well as the mRNA expression of the Th1/Th2/Th17 and Treg master regulators.
While it is tempting to analyze individual cytokines when predicting response to a particular biological treatment, the overall immune response is a resultant of opposing factors. The relative magnitude of the end state compared to healthy controls or pre/post treatment status may offer a way to predict response to a particular biologic treatment. We performed an in-silico simulation of TNF-α blockade for each of the 4 immune phenotypes identified by our mathematical algorithm. Globally the largest impact was seen in type 3 group (Th1 high /Th2 high ) where a diminished Th1 and Th17 activation was predicted along with a reciprocal increase in Treg and TGF-β. Therefore this group appeared to have the best immunologic outcome by blocking the TNF-α. Our findings through mathematical modeling mirror the observations in patients with IBD, where Infliximab treatment down-regulates IL-17 expression in the gut mucosa and promotes healing [45]. On the other hand, type 4 was more consistent with a global immunosuppressed status, and anti-TNF blockade further magnified this changes relative to healthy controls. Our study included samples from symptomatic and/or with evidence of disease activity (abnormal CRP, Calprotectin, anemia) CD patients. All patients were off immunosuppressive therapy at the time of endoscopy for at least 3 month. 8 (50%) of type 4 group patients had previously experienced a failure to at least one biologic treatment. Thus type 4 group (Th1 Low /Th2 Low ) should be a candidate for immune stimulatory therapy rather than further immunosuppression. Also, importantly, the paradoxical presence of clinical and endoscopic activity in this group underscores the importance of global analysis rather than individual immune mediator analysis. In conclusion we propose that mathematical modeling of immune regulatory networks in patients with IBD may allow identifying distinct groups with high relevance to biological therapies. Based on the results of the in-silico simulation of TNF blockade, a multi-pronged approach aimed at several distinct pathways might be superior to biologic monotherapy. In certain cases (type 4) immune activation rather than suppression might be recommended.
Although parameters (coefficients for differential equations) were obtained from other models, they are representative of a chronic inflammatory process and thus relevant to IBD. Sarcoidosis and Lupus (models referenced in our study) also share genetic susceptibility loci with Crohn's disease. Ultimately, validation and refinement of the mathematical parameters will require a prospective approach in IBD patients. Our current model is supposed to provide a basis for such future studies.
We acknowledge that more cytokines could have been included in the present study. Coefficients utilized in the present mathematical model were also derived from prior models of inflammation, which did not include IL-5, IL-13, IL-17 and IL-22. Naïve T helper cells can be induced to differentiate towards T helper 1 (Th1), Th2, Th17 and regulatory (Treg) phenotypes according to the local cytokine milieu. Stable expression of specific transcription factors, T-bet (Th1), GATA-3 (Th2), Foxp3 (Tregs) and RORγt (Th17) will ultimately represent the net effect of the combine cytokine production. Thus the inclusion of all the master regulators relevant to the T cell polarization is expected to reflect the tissue status of these pathways. This potential extension of the modeling with inclusion of the extra cytokines should be considered in future work when more patients data become available.
Our present study represents a snapshot of CD patients evaluated for assessment of disease activity. We intentionally excluded patients with recent exposure to immunosuppressant therapy to avoid a direct influence of the chosen markers. Nevertheless the purpose of the current study was to establish a mathematical model based on IBD patient samples that will serve as a basis for prospective studies.