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.

Introduction

Inflammatory Bowel Diseases (IBD), Crohn’s Disease (CD) and Ulcerative Colitis (UC) result from an inappropriate immune response toward commensal flora [1]. Genome wide association studies (GWAS) indicated that the majority of IBD susceptibility loci belong to immunoregulatory networks [2, 3].

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.

Materials and Methods

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/cm3 and were based on a theoretical density in a cm3 of tissue. The variables (concentrations) included in the model are listed in Table 1.

Fig 1

Schematic diagram of immune system with application to inflammatory bowel disease (IBD).

Table 1

List of the variables included in the model.

M1: activated M1 macrophages	T1: Th1 cells
M2: activated M2 macrophages	T2: Th2 cells
M: total macrophages (M1 + M2)	T17: Th17 cells
Iα: TNF-α	Iβ: TGF-β
Iγ: IFN-γ	I2: IL-2
I4: IL-4	I6: IL-6
I10: IL-10	I12: IL-12
I21: IL-21

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 ν, then, at a given Th1 cell density T, we can write dI/dt = νT. Similarly, the production of this cytokine by activated macrophages would be described by the differential equation dI/dt = νM. 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

IL-6 and TNF-α are produced by M1 [14, 15] and Th1 cells also produce TNF-α [16]; hence, the equations of I 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) 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) and Treg (T): n (I/(ζ + I)) T where ζ 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/ζ) where ζ 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 (f1) that contribute to macrophage polarization; tissue concentrations of TNF-α (I), IFN-γ (I) and TGF-β (I); cytokine specific rates of macrophage activation (σ σ σ) and transition between M1 and M2 phenotypes (σ, σ), and macrophage apoptosis (decay) rate (μ). 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/dt, after including death rate (μ) 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).

Th17 cells: Th17 auto-activation involves IL-6, IL-21 and TGF-β. Th17 cells induce production of IL-21 which, together with IL-6 and TGF-β, may further stimulate the activation of Th17 cells [8, 14, 15]. This is shown in Fig 1(C).

Treg cells: The loop of Treg auto-activation includes TGF-β, IL-10 and IL-2. As shown in Fig 1(D), TGF-β and IL-10 trigger Treg activation while Treg cells express TGF-β and IL-10, forming a positive feedback loop in Treg differentiation [4, 10, 17]. We assume that IL-2 is critical for the activation of Treg cells [23].

Aside from activation, the system is regulated by several inhibitory processes (red lines in Fig 1). Treg inhibits the activation of Th1 and Th2 cells [24] while the pairs, Th1-Th2 and Th17-Treg, are mutually antagonistic [9, 25]. IFN-γ produced by Th1 and IL-4 produced by Th2 inhibits Th17 activation [26].

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 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 σ, σ, σ and σ 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 σ = 10σ. We also considered that the transition rates between M1 and M2 phenotypes are the same as the activation rate of M2, and take σ = σ and σ = σ.

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 ν = 5ν. We also assumed a similar production rate of IL-2 (ν) and IFN-γ (ν) by Th1 cells, so that ν = ν.

It was previously shown that Th2 cells at a density of 3 × 10−2 g/cm3 produce an IL-4 concentration of 15 × 10−9 g/cm3 [31]. Given δ = 349.37 week−1 [7, 32] and using the steady-state equation (this is expressed as νT − δ I = 0), we get ν = 1.75 × 10−4 week−1. We further assumed that Th2 cells produce more IL-4 than macrophages, and take ν = ν/3 = 5.83 × 10−5 week−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 ν = 3.72 × 10−4 week−1 [29, 33] and ν = 3ν.

T cells

Maintenance of Th1 cells pool would require the activation and decay coefficients to be of similar magnitude. The coefficient σ (representing the IL-2 maximal signal output) was estimated to be 1.23 week−1 while μ (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 σ = σ. Similarly, IL-10 and TGF-β contributions to Treg activation/survival rate were taken to be of equal magnitude, so that σ = σ.

Other parameters

We take the inhibition of Th2 by Th1 to be γ = 1.83 × 10−1 g/cm3, the inhibition of Treg by Th17 to be γ = 3.37 × 10−1 g/cm3, and assume that the inhibition of Th1 by Th2 is larger, taking γ = 5.35 × 10−2 g/cm3. We also assume that Treg inhibition effect on the other T cells is still somewhat larger, taking γ = γ = γ = 6.06 × 10−2 g/cm3.

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-9g/cm3; we take the IL-6 concentration in healthy tissue to be 8.00 × 10−6 g/cm3. We can then compute λc and then also the concentrations of TNF-α, IL-10 in healthy tissue. We accordingly get the steady-state concentration of TNF-α to be 9.75 × 10−6 g/cm3 and of IL-10 to be 1.54 × 10−6 g/cm3.

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, v, σ, σ, σ 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, v, σ, σ, σ 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 λTTh1, λTTh2, λTTh17 and λTTreg 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 2

Parameters used in the equations of macrophages and cytokines.

Parameter	Definition	Value	References
Macrophages
f1	Basal activation rate of M1	1 week-1	[29] & estimated
f2	Basal activation rate of M2	10 week-1	[29] & estimated
σMα	Activation rate of macrophages by TNF-α	2.4 week-1	[29]
σMγ	Transition rate from M2 to M1	24 week-1	[29] & estimated
σM10	Activation rate of macrophages by IL-10	24 week-1	[29] & estimated
σMβ	Transition rate from M1 to M2	24 week-1	[29] & estimated
M0	Density of inactivated macrophages	1.5 x 10−3 g/cm3	[29, 35]
μM	Decay rate of macrophages	1 week-1	[36]
Cytokines
νγM	Production rate of IFN-γ by macrophages	8.2 x 10−6 week-1	[29, 36]
νγ1	Production rate of IFN-γ by Th1	4.1 x10-5 week-1	[29, 36] & estimated
ν4M	Production rate of IL-4 by macrophages	5.83 x10-5 week-1	[31] & estimated
ν42	Production rate of IL-4 by Th2	1.75 x 10−4 week-1	[31] & estimated
ν6M	Production rate of IL-6 by macrophages	3.63 x 10−3 week-1	Estimated from data
ν10M	Production rate of IL-10 by macrophages	3.72 x 10−4 week-1	[29, 33]
ν10r	Production rate of IL-10 by Treg	1.12 x 10−3 week-1	[33] & estimated
ν12M	Production rate of IL-12 by macrophages	2.65 x 10−4 week-1	[29, 36]
ν2117	Production rate of IL-21 by Th17	8.05 x 10−4 week-1	[15]
ναM	Production rate of TNF-α by macrophages	2.10 x 10−2 week-1	[29]
να1	Production rate of TNF-α by Th1	7.35 x 10−2 week-1	Estimated from data
νβM	Production rate of TGF-β by macrophages	5.6 x 10−9 week-1	[29]
νβr	Production rate of TGF-β by Treg	3.90 x 10−8 week-1	[37]
ν21	Production rate of IL-2 by Th1	4.1 x 10−5 week-1	[29, 36] & estimated
n2r	Production ratio of IL-10 by IL-2	3	Estimated
δγ	Degradation rate of IFN-γ	29.12 week-1	[7, 18]
δ4	Degradation rate of IL-4	349.37 week-1	[7, 32]
δ6	Degradation rate of IL-6	29.11 week-1	[38]
δ10	Degradation rate of IL-10	116.48 week-1	[39]
δ12	Degradation rate of IL-12	8.33 week-1	[29]
δ21	Degradation rate of IL-21	63.98 week-1	[40]
δα	Degradation rate of TNF-α	388.15 week-1	[41]
δβ	Degradation rate of TGF-β	349.37 week-1	[7, 42]
δ2	Degradation rate of IL-2	537.46 week-1	[7, 43]

Table 3

Parameters used in the equations of T cells and a list of half-saturation constants.

Parameter	Definition	Value	References
T cells
μ1	Decay rate of Th1	1.4 week-1	[29, 36]
μ2	Decay rate of Th2	1.4 week-1	[29, 36]
μ17	Decay rate of Th17	1.4 week-1	[29, 36]
μr	Decay rate of Treg	1.4 week-1	[29, 36]
σ12	Maximum rate of IL-12 signaling	10.93 week-1	Estimated from data
σ2	Maximum rate of IL-2 signaling	1.23 week-1	[29, 36] & estimated
σ4	Maximum rate of IL-4 signaling	1.94 week-1	Estimated from data
σ21	Maximum rate of IL-21 signaling	156.17 week-1	[34] & estimated from data
σ6	Maximum rate of IL-6 signaling	156.17 week-1	[34] & estimated from data
σβ	Maximum rate of TGF-β signaling	14.02 week-1	Estimated from data
σ10	Maximum rate of IL-10 signaling	14.02 week-1	Estimated from data
γ1	Constants for Th1 inhibition of Th2	1.83 x 10−1 g/cm3	Estimated
γ2	Constants for Th2 inhibition of Th1	5.35 x10-2 g/cm3	Estimated
γr1, γr2, γr17	Constants for Treg inhibition	6.06 x10-2 g/cm3	Estimated
γ17	Constants for Th17 inhibition	3.37 x10-1 g/cm3	Estimated
Half-saturation constants
ζ6	Constant for IL-6	8 x 10−6 g/cm3	The constants are the same as the steady-state concentrations of the cytokines in healthy individuals; in solving for the steady-state equations of (1)-(15)
ζα	Constant for TNF-α	9.75 x 10−6 g/cm3
ζ10	Constant for Il-10	1.54 x 10−6 g/cm3
ζ2	Constant for IL-2	6.86 x 10−8 g/cm3
ζγ	Constant for IFN-γ	2.58 x 10−6 g/cm3
ζ12	Constant for IL-12	4.90 x 10−8 g/cm3
ζ4	Constant for IL-4	9.70 x 10−9 g/cm3
ζ21	Constant for IL-21	4.25 x 10−6 g/cm3
ζβ	Constant for TGF-β	6.77 x 10−12 g/cm3

Table 4

Steady-state concentrations of cytokines, macrophages and T cells in a healthy individual.

	Value (g/cm3)
M1 macrophage	2.34 x 10−2
M2 macrophage	2.50 x 10−2
Th1 cell	4.58 x 10−2
Th2 cell	1.34 x 10−2
Th17 cell	3.37 x 10−1
Treg cell	6.06 x 10−2
IFN-γ	2.58 x 10−6
IL-2	6.86 x 10−8
IL-4	9.70 x 10−9
IL-21	4.25 x 10−6
IL-6	8.00 x 10−6
TNF-α	9.75 x 10−6
IL-10	1.54 x 10−6
TGF-β	6.77x 10−12
IL-12	4.90 x 10−8

For completeness we listed in 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, σ, σ, σ, σ, σ σ, σ, σ and σ, 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 T1 (Fig 2A) and T2 (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.

Table 5

PRCC values of some parameters for the outputs T1 and T2.

	PRCC for T1	PRCC for T2
σMα	0.3298*	-0.0190
σM10	0.6264*	-0.1448*
σ12	0.6206*	0.8420*
σ2	0.2257*	0.4404*
σ4	-0.3123*	-0.4336*
σ21	0.2756*	0.3192*
σ6	0.4181*	0.4954*
σβ	-0.6690*	-0.7636*
σ10	-0.6866*	-0.7858*

* Denotes significant PRCC values (p-value <0.01).

Fig 2

Scatter plots of rank transformed (A) and (B) after 100 weeks versus some rank transformed parameters with statistically significant correlation (|PRCC|>0.5 and p-value <0.01). The title of each subplot shows its PRCC value.

All the nine parameters have significant PRCC values either with T1 or with T2. Among the nine parameters, we find that the statistically significant PRCC values with both T1 and T2 are σ, σ and σ. The parameters σ and σ are negatively correlated to T1 and T2: the two 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, σ is positively correlated to T1 and T2: 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 (σ, σ and σ) 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.

Table 6

Baseline characteristics for patients with IBD-Crohn’s Disease.

Total number	N = 58
Gender
Female	27(46.55%)
Male	31(53.45%)
Location
Upper GI	0(0%)
Ileal	0(0%)
Ileo-colonic	47(81.03%)
Colonic	11(18.97%)
Perianal	6(10.34%)
Behavior
Inflammatory	41(70.68%)
Stricturing	14(24.13%)
Fistulizing	3(5.17%)
Current immunosuppressant treatment	0(0%)
Former use of immunomodulator	27(46.55%)
Former use of biologic therapy	19(32.75%)

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: 1st spin @ 7000x for 40 sec and 2nd 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 NanoString 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:

T: the mean value of the Th1 master regulator T-bet of healthy individuals

T: the value of the Th1 master regulator T-bet of IBD patient

T: the mean value of the Th2 master regulator GATA-3 of healthy individuals

T: the value of the Th2 master regulator GATA-3 of IBD patient

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 T-bet, GATA3, RORγt and Foxp3. Based on the observed values we were able to divide the patients into four categories:

T > T, T < T

T < T, T > T

T > T, T > T

T < T, T < T

For example, if an IBD patient had a T-bet level higher than T and a GATA-3 level lower than T, 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.

Table 7

Fold changes of the cytokine and T cell concentrations obtained from the clinical data and the simulations in different types of diseases.

Type of diseases	# of cases	IL-6	IL-10	TNF-α	T-Bet	Gata3	RORγt	Foxp3	TGF-β*	IFN-γ*
Th1↑ Th2↓	7	↓-93%	↓-51%	↓-31%	↑+50%	↓-22%	↑+50%	↓-71%	↓-56%	↓-99%
Th1↓Th2↑	18	↓-74%	↓-15%	↓-31%	↓-29%	↑+50%	↓-27%	↑+1%	↓-6%	↓-92%
Th1↓Th2↑	18	↓-74%	↓-15%	↓-31%	↓-29%	↑+50%	↓-27%	↑+1%	↓-6%	↓-92%
Th1↑Th2↑	17	↑+35%	↑+44%	↑+47%	↑+110%	↑+75%	↓-14%	↑+70%	↑+19%	↑+481%
Th1↓Th2↓	16	↓-83%	↓-34%	↓-45%	↓-46%	↓-40%	↓-11%	↓-42%	↓-20%	↓-95%

* The values of TGF-β and IFN-γ are predicted by the model.

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, ν, ν, ν, ν, σ, σ, σ, σ, σ and σ (Table 8).

Table 8

Simulation results: Parameter variations in different types of diseases.

Type of diseases	νγ1	να1	νβρ	ν10ρ	σ12	σ4	σ21	σ6	σβ	σ10
	coefficient of IFN-γ by Th1	coefficient of TNF-α by Th1	coefficient of TGF-β by Treg	coefficient of IL-10 by Treg	activation of Th1 cells	activation of Th2 cells	activation of Th17 cells	activation of Th17 cells	activation of Treg cells	activation of Treg cells
Th1↑Th2↓	↓-99.6%	↓-50%	↑+35%	↑+35%	↑+298%	↓-32%	↑+49%	↑+49%	↓-50%	↓-50%
Th1↓Th2↑	↓-90%	↑+4%	↓-10%	↓-10%	↑+81%	↑+25%	↓-13%	↓-13%	-0%	-0%
Th1↑Th2↑	↑+194%	↓-29%	↓-30%	↓-30%	↑+255%	↑+122%	↑+298%	↑+298%	↑+17%	↑+17%
Th1↓Th2↓	↓-92%	↑+11%	↑+32%	↑+32%	↑+5%	↓-48%	↓-28%	↓-28%	↓-18%	↓-18%

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-α

Currently, anti-TNF antibodies represent the main group of biologic agents for the treatment of IBD. Therefore we simulated the TNF-α blockade in our mathematical model. For modeling the TNF-α blockade, we fix the concentration of TNF-α (Iα) to be zero and calculate the change of other variables at the steady-state concentrations. Fig 3 shows the predicted mRNA fold change in cytokines and T cell master regulators. In type 1 group, TNF-α blockade had the most profound effect on Th1 (93%) and Th17 (40%) master regulators. In the third group our model predicted a significant decrease in Th1 (105%) and increase in Treg (108%). In the second and fourth group the decrease in both Th1 (35%, 34%) and Th17 (19%, 21%) master regulators was modest whereas Treg was virtually unaffected (1%, 2%).

Fig 3

Simulation results: fold changes of the cytokine and T cell concentrations when TNF-α is completely blocked in different types of diseases.

Blue bars represent the results of pre-treatment; red bars represent the results of post-treatment with TNF-α blockage.

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 (Th1high/Th2high) 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 (Th1Low/Th2Low) 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.