Stochastic Effects in Autoimmune Dynamics

Among various possible causes of autoimmune disease, an important role is played by infections that can result in a breakdown of immune tolerance, primarily through the mechanism of “molecular mimicry”. In this paper we propose and analyse a stochastic model of immune response to a viral infection and subsequent autoimmunity, with account for the populations of T cells with different activation thresholds, regulatory T cells, and cytokines. We show analytically and numerically how stochasticity can result in sustained oscillations around deterministically stable steady states, and we also investigate stochastic dynamics in the regime of bi-stability. These results provide a possible explanation for experimentally observed variations in the progression of autoimmune disease. Computations of the variance of stochastic fluctuations provide practically important insights into how the size of these fluctuations depends on various biological parameters, and this also gives a headway for comparison with experimental data on variation in the observed numbers of T cells and organ cells affected by infection.

1. Introduction

Breakdown of immune tolerance and the resulting autoimmune disease occur when the immune system fails to distinguish the host's own healthy cells from the cells affected by the infection, thus triggering an immune response that also targets healthy cells. Autoimmune disease is usually focused in a specific organ or part of the body, such as retina in the case of uveitis, central nervous system in multiple sclerosis, or pancreatic β-cells in insulin-dependent diabetes mellitus type-1 (Prat and Martin, 2002; Kerr et al., 2008; Santamaria, 2010). Whilst it is close to impossible to pinpoint precise causes of autoimmunity in each individual case, it can usually be attributed to a number of factors, which can include the genetic predisposition, age, previous immune challenges, exposure to pathogens etc. A number of distinct mechanisms have been identified for how an infection of the host with a pathogen can result in the subsequent onset of autoimmune disease, and these include bystander activation (Fujinami, 2011) and molecular mimicry (von Herrath and Oldstone, 1996; Ercolini and Miller, 2008), which is particularly important in the context of autoimmunity caused by viral infections.

Over the years, a number of mathematical models have investigated various origins and aspects of immune response, with an emphasis on the onset and the development of autoimmune disease. Some of the earlier models studied interactions between regulatory and effector T cells without looking at causes of autoimmunity but instead focusing on T cell vaccination (Segel et al., 1995). Borghans et al. (Borghans and De Boer, 1995; Borghans et al., 1998) looked into this process in more detail and showed the onset of autoimmunity, which was defined as oscillations in the number of autoreactive cells that exceeded a certain threshold. León et al. (2000, 2003, 2004) and Carneiro et al. (2005) have analyzed interactions between different T cells and their effect on regulation of immune response and control of autoimmunity. More recently, Iwami et al. (2007, 2009) considered a model of immune response to a viral infection, in which they explicitly included the dynamics of a virus population. Although this model is able to demonstrate an emergence of autoimmunity, it fails to produce a regime of normal viral clearance. Alexander and Wahl (2011) have focused on the role of professional antigen-presenting cells (APCs) and their interactions with regulatory and effector effector cells for the purposes of controlling autoimmune response. Burroughs et al. (2011a,b) have analyzed the autoimmunity through the mechanism of cytokine-mediated bystander activation. A special issue on “Theories and modeling of autoimmunity” provides an excellent overview of recent research in the area of mathematical modeling of various aspects of onset and development of autoimmune disease (Root-Bernstein, 2015).

These are several different frameworks for modeling the role of T cells in controlling autoimmune response. Alexander and Wahl (2011) and Burroughs et al. (2011a,b) have explicitly included a separate compartment for regulatory T cells that are activated by autoantigens and suppress the activity of autoreactive T cells. Another modeling approach is to consider the possibility of the same T cells performing different immune functions through having different or tunable activation thresholds, which allows T cells to adjust their response to T cell antigen receptor stimulation by autoantigens. This methodology was originally proposed theoretically to study peripheral and central T cell activation (Grossman and Paul, 1992, 2000; Grossman and Singer, 1996), and has been subsequently used to analyse differences in activation/response thresholds that are dependent on the activation state of the T cell (Altan-Bonnet and Germain, 2005). van den Berg and Rand (2004) and Scherer et al. (2004) have studied stochastic tuning of activation thresholds. Interestingly, the need for T cells to have tunable activation can be shown to emerge from the fundamental principles of of signal detection theory (Noest, 2000). A number of murine and human experiments have confirmed that activation of T cells can indeed dynamically change during their circulation (Nicholson et al., 2000; Bitmansour et al., 2002; Stefanova et al., 2002; Römer et al., 2011), thus supporting the theory developed in earlier papers.

Since immune response is known to be a complex multi-factor process (Perelson and Weisbuch, 1997), a number of studies have looked into various stochastic aspects of immune dynamics, such as T cell selection and proliferation. Deenick et al. (2003) have analyzed stochastic effects of interleukin-2 (IL-2) on T cell proliferation from precursors. Blattman et al. (2000) have shown that repertoires of the CTL (cytotoxic T cell lymphocyte) populations during primary response to a viral infection and in the memory pool are similar, thus providing further support to the theory of stochastic selection for the memory pool. Detours and Perelson (2000) have explored the distribution of possible outcomes during T cell selection with account for variable affinity between T cell receptors and MHC-peptide complexes. Chao et al. (2004) analyzed a detailed stochastic model of T cell recruitment during immune response to a viral infection. Stirk et al. (2010a,b) have developed a stochastic model for T cell repertoire and investigated the role of competitive exclusion between different clonotypes. Using the methodology of continuous-time Markov processes, the authors computed extinction times, a limited multivariate probability distribution, as well as the size of fluctuations around the deterministic steady states. Reynolds et al. (2012) have used a similar methodology to investigate an important question of asymmetric cell division and its impact on the extinction of different T cell populations and the expected lifetimes of naïve T cell clones. With regards to modeling autoimmune dynamics, Alexander and Wahl (2011) have studied the stochastic model of immune response with an emphasis on professional APCs to show that the probability of developing a chronic autoimmune response increases with the initial exposure to self-antigen or autoreactive effector T cells. An important aspect of stochastic dynamics that has to be accounted for in the models is the so-called stochastic amplification (Alonso et al., 2007; Kuske et al., 2007), which denotes a situation where periodic solutions with decaying amplitudes in the deterministic model can result in sustained stochastic periodic oscillations in individual realizations of the same model. This suggests that whilst on average the behavior may show decaying-amplitude oscillations, individual realizations represented by stochastic oscillations can explain relapses/remissions in clinical manifestations of the disease as caused by endogenous stochasticity of the immune processes.

Blyuss and Nicholson (2012, 2015) have proposed and analyzed a mathematical model of immune response to a viral infection that explicitly takes into account the populations of two types of T cells with different activation thresholds and also allows for infection and autoimmune response to occur in different organs. This model supports the regimes of normal viral clearance, a chronic infection, and an autoimmune state represented by exogenous oscillations in cell populations, associated with episodes of high viral production followed by long periods of quiescence. Such behavior, that in the clinical observation could be associated with relapses and remissions, has been observed in a number of autoimmune diseases, such as MS, autoimmune thyroid disease and uveitis (Ben Ezra and Forrester, 1995; Davies et al., 1997; Nylander and Hafler, 2012). Despite its successes, this model has a limitation that the periodic oscillations are only possible when the amount of free virus and the number of infected cells are also exhibiting oscillations, while in laboratory and clinical situations, one rather observes a situation where the initial infection is completely cleared, and this is then followed by the onset of autoimmune reaction. To overcome this limitation, Fatehi et al. (unpublished) have recently extended the model of Blyuss and Nicholson to also include the population of regulatory T cells and the cytokine mediating T cell activity.

In this paper we analyse the effects of stochasticity on the dynamics of immune response in a model with the populations of T cells with different activation thresholds, regulatory cells and cytokines, as presented in Methods. Starting with a system of ordinary differential equations, we apply the methodology of continuous-time Markov chains (CTMC) to derive a Kolmogorov, or chemical master equation, describing the dynamics of a probability distribution of finding the system in a particular state. To make further analytical and numerical progress, we derive an Itô stochastic differential equation, whose solutions provide similar stochastic paths to those of the CTMC models. This then allows us to numerically study the stationary multivariate probability distributions for the states in the model, explore stochastic amplification, determine how the magnitude of stochastic fluctuations around deterministic steady states depends on various parameters, and investigate the effects of initial conditions on the outcome in the case of bi-stability between different dynamical states. These results suggest that the experimentally observed variation in the progression of autoimmune disease can be attributed to stochastic amplification, and they also provide insights into how the variance of fluctuations depends on parameters, which can guide new laboratory experiments.

2. Methods

2.1. Continuous-Time Markov Chain Model of Immune Dynamics

In a recent paper we introduced and analyzed a deterministic model for autoimmune dynamics with account for the populations of T cells with different activation thresholds and cytokines (Fatehi et al. unpublished). The analysis showed that depending on parameters and initial conditions, the model can support the regimes of normal disease clearance, where an initial infection is cleared without further consequences for immune dynamics, chronic infection characterized by a persistent presence of infected cells in the body, and the state of autoimmune behavior where after clearance of initial infection, the immune system supports stable endogenous oscillations in the number of autoreactive T cells, which can be interpreted in the clinical practice of autoimmune disease as periods of relapses and remissions. This work extended earlier results on modeling the effects of tunable activation thresholds (Blyuss and Nicholson, 2012, 2015) by including regulatory T cells, as well as the cytokine mediating proliferation and activity of different types of T cells.

A deterministic model for immune response to a viral infection, as illustrated in a diagram shown in Figure 1, has the form

$\begin{array}{l}\frac{dS}{dt}=rS\left(1-\frac{S}{N}\right)-\beta SF-{\mu }_a{T}_{aut}S,\\ \frac{dF}{dt}=\beta SF-d_{F}F-{\mu }_F{T}_{nor}F-{\mu }_a{T}_{aut}F,\\ \frac{d{T}_{in}}{dt}={\lambda }_{in}-d_{in}{T}_{in}-\alpha T_{in}F,\\ \frac{d{T}_{reg}}{dt}={\lambda }_r-d_r{T}_{reg}+p_1\alpha T_{in}F+{\rho }_{1}I{T}_{reg},\\ \frac{d{T}_{nor}}{dt}=p_2\alpha T_{in}F-d_n{T}_{nor}+{\rho }_{2}I{T}_{nor},\\ \frac{d{T}_{aut}}{dt}=\left(1-p_1-p_2\right)\alpha T_{in}F-d_a{T}_{aut}-\delta T_{reg}{T}_{aut}+{\rho }_{3}I{T}_{aut},\\ \frac{dI}{dt}={\sigma }_1{T}_{nor}+{\sigma }_2{T}_{aut}-d_{i}I,& (1)\end{array}$

where S(t) is the number of susceptible organ cells, F(t) is the number of infected cells, T_in(t) is the number of naïve T cells, T_reg(t) is the number of regulatory T cells, T_nor(t) is the number of activated T cells which recognize foreign antigen and destroy infected cells, T_aut(t) is the number of autoreactive T cells which destroy cells presenting both foreign and self-antigen, and I(t) is the amount of interleukin 2 (IL-2) cytokine. In this model, it is assumed that in the absence of infection, organ cells in the host reproduce logistically with a linear growth rate r and carrying capacity N, and they can become infected at rate β by already infected cells that are producing new virus particles. Unlike earlier models (Blyuss and Nicholson, 2012, 2015; Fatehi et al. unpublished), we consider the situation where the process of producing virions by infected cells is quite fast, hence, we do not explicitly incorporate a separate compartment for free virus. Regarding immune response, we assume that naïve T cells remain in homeostasis, and upon activation at rate α by a signal from infected cells, a proportion p₁ of them will develop into regulatory T cells, a proportion p₂ will become normal activated T cells able to destroy infected cells at rate μ_F, and the remaining T cells will become autoreactive, in which case their threshold for activation by susceptible cells is reduced, and hence, they will be destroying both infected and susceptible host cells at rate μ_a. The effect of regulatory T cells is in reducing the number of autoreactive T cells at rate δ, and regulatory T cells are themselves assumed to be in a state of homeostasis. Finally, normal and autoreactive T cells produce IL-2 at rates σ₁ and σ₂, and IL-2 in turn facilitates proliferation of regular, normal and autoreactive T cells at rates ρ₁, ρ₂, and ρ₃, respectively. One should note that in light of experimental evidence suggesting the possibility of autoimmunity in the absence of B cells (Wolf et al., 1996) and the fact that the development of antibodies can itself depend on prior T cell activation with bacteria (Wu et al., 2010), the above model does not take into account antibody response, but rather focuses on T cell dynamics.

FIGURE 1

Figure 1. A schematic diagram of immune response to infection. Blue indicates host cells (susceptible and infected), red denotes T cells (naïve, regulatory, normal activated, and autoreactive), yellow shows cytokines (interleukin-2).

As a first step in the analysis of stochastic effects in immune dynamics, we construct a CTMC model based on the ODE model (1) using the methodology developed earlier in the context of modeling stochastic effects in epidemic and immunological models (Brauer et al., 2008; Stirk et al., 2010a; Allen, 2011). To this end, we introduce variables X₁(t), …, X₇(t) ∈ {0, 1, 2, …} as discrete random variables representing the number of uninfected cells, infected cells, naïve T cells, regulatory T cells, normal activated T cells, autoreactive T cells, and interleukin-2 at time t, respectively. Let the initial condition be fixed as

$\begin{array}{l}{\text{X}}_0=\left(X_1\left(0\right),\dots ,X_7\left(0\right)\right)=\left(n_{10},n_{20},n_{30},n_{40},n_{50},n_{60},n_{70}\right).\end{array}$

The probability of finding the system in the state n = (n₁, n₂, n₃, n₄, n₅, n₆, n₇) with n_i ∈ {0, 1, 2, …} at time t can be defined as

$\begin{array}{l}P\left(\text{n},t\right)=\mathrm{\text{Prob}}\left\{\text{X}\left(t\right)=\text{n}|\text{X}\left(0\right)={\text{X}}_0\right\}.\end{array}$

Let Δt be sufficiently small such that ΔX_i(t) = X_i(t + Δt)−X_i(t) ∈ {−1, 0, 1} for 1 ≤ i ≤ 7. The CTMC can then be formulated as a birth and death process in each of the variables (Allen, 2011). The infinitesimal transition probabilities corresponding to Figure 1 are as follows,

$\begin{array}{ll}\mathrm{\text{Prob}}\left(\Delta \text{X}=\text{i}|\text{X}=\text{n}\right)=\{\begin{array}{cc}{q}_1\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(1,0,0,0,0,0,0\right),\hfill \\ q_2\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(-1,0,0,0,0,0,0\right),\hfill \\ q_3\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(-1,1,0,0,0,0,0\right),\hfill \\ q_4\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,1,0,0,0,0\right),\hfill \\ q_5\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,-1,0,0,0,0\right),\hfill \\ q_6\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,-1,0,1,0,0\right),\hfill \\ q_7\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,-1,0,1,0,0\right),\hfill \\ q_8\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,-1,0,0,1,0\right),\hfill \\ q_9\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,-1,0,0,0,0,0\right),\hfill \\ q_{10}\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,0,1,0,0,0\right),\hfill \\ q_{11}\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,0,-1,0,0,0\right),\hfill \\ q_{12}\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,0,0,1,0,0\right),\hfill \\ q_{13}\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,0,0,-1,0,0\right),\hfill \\ q_{14}\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,0,0,0,1,0\right),\hfill \\ q_{15}\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,0,0,0,-1,0\right),\hfill \\ q_{16}\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,0,0,0,0,1\right),\hfill \\ q_{17}\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,0,0,0,0,-1\right),\hfill \\ 1-{\displaystyle \sum _{i=1}^{17}}{q}_i\Delta t+o\left(\Delta t\right),\hfill & \text{i}=\left(0,0,0,0,0,0,0\right),\hfill \\ o\left(\Delta t\right),\hfill & \mathrm{\text{otherwise}},\hfill \end{array}\hfill & (2)\end{array}$

where

$\begin{array}{l}{q}_1=b_1{n}_1+b_2{n}_1^2,q_2=d_1{n}_1+d_2{n}_1^2+{\mu }_a{n}_1{n}_6,q_3=\beta n_1{n}_2,q_4={\lambda }_{in},\\ q_5=d_{in}{n}_3,q_6=p_1\alpha n_2{n}_3,q_7=p_2\alpha n_2{n}_3,q_8=\left(1-p_1-p_2\right)\alpha n_2{n}_3,\\ q_9=\left(d_F+{\mu }_F{n}_5+{\mu }_a{n}_6\right)n_2,q_{10}={\lambda }_r+{\rho }_1{n}_4{n}_7,q_{11}=d_r{n}_4,q_{12}={\rho }_2{n}_5{n}_7,\\ q_{13}=d_n{n}_5,q_{14}={\rho }_3{n}_6{n}_7,q_{15}=\left(d_a+\delta n_4\right)n_6,q_{16}={\sigma }_1{n}_5+{\sigma }_2{n}_6,q_{17}=d_i{n}_7.\end{array}$

Here, $b_1{n}_1+b_2{n}_1^2$ and $d_1{n}_1+d_2{n}_1^2$ are natural birth and death rates for uninfected cells with b₁ − d₁ = r and d₂ − b₂ = r/N (Allen, 2011).

The probabilities P(n, t) satisfy the following master equation (forward Kolmogorov equation)

$\begin{array}{l}\frac{dP(n,t)}{dt}=\left\{\right({\epsilon }_1^{-}-1)q_1+({\epsilon }_1^{+}-1)q_2+({\epsilon }_1^{+}{\epsilon }_2^{-}-1)q_3+({\epsilon }_3^{-}-1)q_4+({\epsilon }_3^{+}-1)q_5\\ +({\epsilon }_3^{+}{\epsilon }_4^{-}-1)q_6+({\epsilon }_3^{+}{\epsilon }_5^{-}-1)q_7+({\epsilon }_3^{+}{\epsilon }_6^{-}-1)q_8+({\epsilon }_2^{+}-1)q_9+({\epsilon }_4^{-}-1)q_{10}\\ +({\epsilon }_4^{+}-1)q_{11}+({\epsilon }_5^{-}-1)q_{12}+({\epsilon }_5^{+}-1)q_{13}+({\epsilon }_6^{-}-1)q_{14}+({\epsilon }_6^{+}-1)q_{15}\\ +({\epsilon }_7^{-}-1)q_{16}+({\epsilon }_7^{+}-1)q_{17}\left\}P\right(n,t).& (3)\end{array}$

where the operators ${\epsilon }_i^{\pm }$ are defined as follows,

$\begin{array}{l}{\epsilon }_i^{\pm }f\left(n_1,n_2,n_3,n_4,n_5,n_6,n_7,t\right)=f\left(n_1,\dots ,n_i\pm 1,\dots ,n_7,t\right),\end{array}$

for each 1 ≤ i ≤ 7, and if n_i < 0 for any 1 ≤ i ≤ 7, then P(n, t) = 0.

By solving this master equation, one can find the probability density function for this model. However, since this is a high-dimensional difference-differential equation, solving it is a very challenging task. Normally, the number of events occurring in a small time step in the CTMC model is extremely large, hence using the CTMC model for plotting stochastic trajectories is very computationally intensive (Mandal et al., 2014). A much more computationally efficient approach is to use chemical Langevin equations (Gillespie, 2000, 2002), also known as Itô stochastic differential equation (SDE) models, which provide very similar sample paths to those of the CTMC models (Mandal et al., 2014). While both Itô and Stratonovich interpretations of stochastic calculus can be applied (Øksendal, 2000), in biological applications Itô formulation is more frequently used due to its non-anticipatory nature and a closer connection to numerical implementation (Allen, 2007, 2011; Braumann, 2007).

2.2. Stochastic Differential Equation Model

To derive Itô SDE model, let Y(t) = (Y₁(t), Y₂(t), Y₃(t), Y₄(t), Y₅(t), Y₆(t), Y₇(t)) be a continuous random vector for the sizes of various cell compartments at time t. Similar to the CTMC model, we assume that Δt is small enough so that during this time interval at most one change can occur in state variables. These changes together with their probabilities are listed in Table 1, which is again based on Figure 1 and transitions in the CTMC model (2). Using this table of possible state changes, one can compute the expectation vector and covariance matrix of ΔY for sufficiently small Δt (Allen et al., 2008; Mandal et al., 2014). The expectation vector to order Δt is given by

$\begin{array}{l}𝔼\left(\Delta \text{Y}\right)\approx {\displaystyle \sum _{i=1}^{17}}{P}_i{\left(\Delta \text{Y}\right)}_i\Delta t=\mu \Delta t,\end{array}$

where

$\begin{array}{l}\mu =\left(\begin{array}{c}{P}_1-P_2-P_3\\ P_3-P_9\\ P_4-P_5-P_6-P_7-P_8\\ P_6+P_{10}-P_{11}\\ P_7+P_{12}-P_{13}\\ P_8+P_{14}-P_{15}\\ P_{16}-P_{17}\end{array}\right)\end{array}$

is the drift vector, which can be easily seen to be identical to the right-hand side of the deterministic model Equation (1). The covariance matrix is obtained by keeping terms of order Δt only, i.e.,

$\begin{array}{l}\begin{array}{c}\mathrm{\text{cov}}\left(\Delta \text{Y}\right)=𝔼\left[\left(\Delta \text{Y}\right){\left(\Delta \text{Y}\right)}^T\right]\hfill \\ -𝔼\left[\Delta \text{Y}\right]{\left(𝔼\left[\Delta \text{Y}\right]\right)}^T\approx 𝔼\left[\left(\Delta \text{Y}\right){\left(\Delta \text{Y}\right)}^T\right]\hfill \\ ={\displaystyle \sum _{i=1}^{17}}{P}_i{\left(\Delta \text{Y}\right)}_i{\left(\Delta {\text{Y}}_i\right)}^T\Delta t=\Sigma \Delta t,\hfill \end{array}\hfill \end{array}$

where

$\begin{array}{l}\Sigma =\left(\begin{array}{ccccccc}{P}_1+P_2+P_3& -P_3& 0& 0& 0& 0& 0\\ -P_3& P_3+P_9& 0& 0& 0& 0& 0\\ 0& 0& P_4+P_5+P_6+P_7+P_8& -P_6& -P_7& -P_8& 0\\ 0& 0& -P_6& P_6+P_{10}+P_{11}& 0& 0& 0\\ 0& 0& -P_7& 0& P_7+P_{12}+P_{13}& 0& 0\\ 0& 0& -P_8& 0& 0& P_8+P_{14}+P_{15}& 0\\ 0& 0& 0& 0& 0& 0& P_{16}+P_{17}\end{array}\right)\end{array}$

is a 7 × 7 covariance matrix. To derive Itô SDE model, we need to find a diffusion matrix H defined according to HH^T = Σ. Although this matrix is not unique, different forms of this matrix give equivalent systems (Allen, 2007; Allen et al., 2008).

TABLE 1

Table 1. Possible state changes ΔY during a small time interval Δt.

If one rewrites the covariance matrix Σ in the form

$\begin{array}{l}\Sigma =\left(\begin{array}{ccc}U& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& W& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& Z\end{array}\right),\end{array}$

with

$\begin{array}{l}U=\left(\begin{array}{cc}{P}_1+P_2+P_3& -P_3\\ -P_3& P_3+P_9\end{array}\right),Z=P_{16}+P_{17},\end{array}$

and

$\begin{array}{l}W=\left(\begin{array}{cccc}{P}_4+P_5+P_6+P_7+P_8& -P_6& -P_7& -P_8\\ -P_6& P_6+P_{10}+P_{11}& 0& 0\\ -P_7& 0& P_7+P_{12}+P_{13}& 0\\ -P_8& 0& 0& P_8+P_{14}+P_{15}\end{array}\right),\end{array}$

we can define three matrices H₁, H₂ and H₃ as follows,

$\begin{array}{l}{H}_1=\left(\begin{array}{ccc}\sqrt{P_1+P_2}& -\sqrt{P_3}& 0\\ 0& \sqrt{P_3}& \sqrt{P_9}\end{array}\right),H_3=\sqrt{P_{16}+P_{17}},\\ H_2=\left(\begin{array}{ccccccc}\sqrt{P_4+P_5}& -\sqrt{P_6}& -\sqrt{P_7}& -\sqrt{P_8}& 0& 0& 0\\ 0& \sqrt{P_6}& 0& 0& \sqrt{P_{10}+P_{11}}& 0& 0\\ 0& 0& \sqrt{P_7}& 0& 0& \sqrt{P_{12}+P_{13}}& 0\\ 0& 0& 0& \sqrt{P_8}& 0& 0& \sqrt{P_{14}+P_{15}}\end{array}\right).\end{array}$

Now if we consider

$\begin{array}{l}H=\left(\begin{array}{ccc}{H}_1& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& H_2& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& H_3\end{array}\right),\end{array}$

then HH^T = Σ, where H is a 7 × 11 matrix. The Itô SDE model now has the form

$\begin{array}{ll}\{\begin{array}{c}d\text{Y}\left(t\right)=\mu dt+Hd\text{W}\left(t\right),\hfill \\ \text{Y}\left(0\right)={\left(A\left(0\right),F\left(0\right),T_{in}\left(0\right),T_{reg}\left(0\right),T_{nor}\left(0\right),T_{aut}\left(0\right),I\left(0\right)\right)}^T,\hfill \end{array}\hfill & (4)\end{array}$

and $\mathbf{\text{W}}(t)={\left[W_1(t),W_2(t),\dots ,W_{11}(t)\right]}^T$ is a vector of 11 independent Wiener processes (Allen, 2007).

In order to make further analytical progress, we find an approximate probability density function for the model (4) as given by an approximate solution of the master equation (van Kampen, 1981; Allen, 2007). Let P(Y, t) be the probability density function of the model (4). Then P(Y, t) satisfies the following Fokker-Planck equation (Gardiner, 2004; Allen, 2007) which is an approximation of the master equation

$\begin{array}{l}\{\begin{array}{c}\frac{\partial P\left(\text{Y},t\right)}{\partial t}=-{\displaystyle \sum _{i=1}^7}\frac{\partial }{\partial y_i}\left[{\mu }_{i}P\left(\text{Y},t\right)\right]+\frac{1}{2}{\displaystyle \sum _{i=1}^7}{\displaystyle \sum _{j=1}^7}\frac{{\partial }^2}{\partial y_i\partial y_j}\left[{\Sigma }_{ij}P\left(\text{Y},t\right)\right],\hfill \\ \hfill \\ P\left(\text{Y},0\right)={\delta }_7\left(\text{Y}-{\text{Y}}_0\right).\hfill \end{array}\hfill \end{array}$

By solving this PDE, one can find the probability density function of our model, but since this equation is high-dimensional and nonlinear, solving it analytically is impossible. Hence, we use another approach, a so-called system size expansion or van Kampen's Ω-expansion (van Kampen, 1981), which is a method for constructing a continuous approximation to a discrete stochastic model (Stirk et al., 2010a,b), which allows one to study stochastic fluctuations around deterministic attractors (Black et al., 2009).

2.3. System Size Expansion

In order to apply the van Kampen's approach, we consider fluctuations within a systematic expansion of the master equation for a large system size Ω. Specifically, we write each n_i(t) as a deterministic part of order Ω plus a fluctuation of order Ω^½ as follows,

$\begin{array}{ll}{n}_i\left(t\right)=\Omega x_i\left(t\right)+{\Omega }^{½}{\zeta }_i\left(t\right),i=1,\dots ,7,& (5)\end{array}$

where x_i(t) and ζ_i(t) are two continuous variables, and Ωx_i(t) = 𝔼[n_i(t)]. The probability density P(n, t) satisfying the master Equation (3) is now represented by the probability density Π(ζ, t), i.e., Π(ζ, t) = P(n, t) = P(Ωx + Ω^½ ζ, t), which implies

$\begin{array}{ll}\frac{dP\left(\text{n},t\right)}{dt}=\frac{\partial \Pi }{\partial t}-{\displaystyle \sum _{i=1}^7}{\Omega }^{½}\frac{d{x}_i}{dt}\frac{\partial \Pi }{\partial {\zeta }_i}.& (6)\end{array}$

To expand the master equation (3) in a power series in Ω^−½, we use the following expansion for the step operators

$\begin{array}{ll}{\epsilon }_i^{\pm }=1\pm {\Omega }^{-½}\frac{\partial }{\partial {\zeta }_i}+\frac{1}{2}{\Omega }^{-1}\frac{{\partial }^2}{\partial {\zeta }_i^2}\pm \cdots .& (7)\end{array}$

Substituting expressions (6, 7) into the master equation (see Supplementary Material for details) and collecting terms of order Ω^½ yields the following deterministic model for macroscopic behavior

$\begin{array}{ll}\begin{array}{c}\frac{d{x}_1}{dt}=b_1{x}_1+{\tilde{b}}_2{x}_1^2-d_1{x}_1-{\tilde{d}}_2{x}_1^2-\tilde{\beta }{x}_1{x}_2-{\tilde{\mu }}_a{x}_1{x}_6,\hfill \\ \frac{d{x}_2}{dt}=\tilde{\beta }{x}_1{x}_2-d_F{x}_2-{\tilde{\mu }}_F{x}_2{x}_5-{\tilde{\mu }}_a{x}_2{x}_6,\hfill \\ \frac{d{x}_3}{dt}={\tilde{\lambda }}_{in}-d_{in}{x}_3-\tilde{\alpha }{x}_2{x}_3,\hfill \\ \frac{d{x}_4}{dt}={\tilde{\lambda }}_r-d_r{x}_4+p_1\tilde{\alpha }{x}_2{x}_3+{\tilde{\rho }}_1{x}_4{x}_7,\hfill \\ \frac{d{x}_5}{dt}=p_2\tilde{\alpha }{x}_2{x}_3-d_n{x}_5+{\tilde{\rho }}_2{x}_5{x}_7,\hfill \\ \frac{d{x}_6}{dt}=\left(1-p_1-p_2\right)\tilde{\alpha }{x}_2{x}_3-d_a{x}_6-\tilde{\delta }{x}_4{x}_6+{\tilde{\rho }}_3{x}_6{x}_7,\hfill \\ \frac{d{x}_7}{dt}={\sigma }_1{x}_5+{\sigma }_2{x}_6-d_i{x}_7,\hfill \end{array}\hfill & (8)\end{array}$

where

$\begin{array}{l}{b}_2=\frac{{\tilde{b}}_2}{\Omega },d_2=\frac{{\tilde{d}}_2}{\Omega },\beta =\frac{\tilde{\beta }}{\Omega },{\mu }_a=\frac{{\tilde{\mu }}_a}{\Omega },\\ {\mu }_F=\frac{{\tilde{\mu }}_F}{\Omega },\alpha =\frac{\tilde{\alpha }}{\Omega },\delta =\frac{\tilde{\delta }}{\Omega },\\ {\rho }_i=\frac{{\tilde{\rho }}_i}{\Omega },i=1,2,3,{\lambda }_{in}={\tilde{\lambda }}_{in}\Omega ,{\lambda }_r={\tilde{\lambda }}_r\Omega .\end{array}$

Model (8) has been analyzed in Fatehi et al. (unpublished), and it can have at most four biologically feasible steady states. The first one, a disease-free steady state, is given by

$\begin{array}{l}{S}_1^{*}=\left(\frac{b_1-d_1}{{\tilde{d}}_2-{\tilde{b}}_2},0,\frac{{\tilde{\lambda }}_{in}}{d_{in}},\frac{{\tilde{\lambda }}_r}{d_r},0,0,0\right),\end{array}$

and it is stable if $d_F>\tilde{\beta }$. The second and third steady states can be found as

$\begin{array}{l}{S}_2^{*}=\left(0,0,\frac{{\tilde{\lambda }}_{in}}{d_{in}},x_4^{*},0,\frac{d_i\left(d_a+\tilde{\delta }{x}_4^{*}\right)}{{\tilde{\rho }}_3{\sigma }_2},\frac{d_a+\tilde{\delta }{x}_4^{*}}{{\tilde{\rho }}_3}\right),\end{array}$

and

$\begin{array}{l}{S}_3^{\ast }=(\frac{{\tilde{\rho }}_3{\sigma }_2(b_1-d_1)-{\tilde{\mu }}_a{d}_i\left(d_a+\tilde{\delta }{x}_4^{\ast }\right)}{{\tilde{\rho }}_3{\sigma }_2\left({\tilde{d}}_2-{\tilde{b}}_2\right)},\\ 0,\frac{{\tilde{\lambda }}_{in}}{d_{in}},x_4^{\ast },0,\frac{d_i\left(d_a+\tilde{\delta }{x}_4^{\ast }\right)}{{\tilde{\rho }}_3{\sigma }_2},\frac{d_a+\tilde{\delta }{x}_4^{\ast }}{{\tilde{\rho }}_3}),\end{array}$

where $x_4^{*}$ satisfies the following quadratic equation

$\begin{array}{ll}{\tilde{\rho }}_1\tilde{\delta }{\left(x_4^{*}\right)}^2+\left({\tilde{\rho }}_1{d}_a-{\tilde{\rho }}_3{d}_r\right)x_4^{*}+{\tilde{\rho }}_3{\tilde{\lambda }}_r=0.& (9)\end{array}$

These steady states are stable, provided

$\begin{array}{l}\begin{array}{c}\frac{{\sigma }_2}{{\tilde{\mu }}_a{d}_i}K<\frac{d_a+\tilde{\delta }{x}_4^{*}}{{\tilde{\rho }}_3}<\frac{d_n}{{\tilde{\rho }}_2},\tilde{\delta }{\tilde{\rho }}_1{\left(x_4^{*}\right)}^2>{\tilde{\lambda }}_r{\tilde{\rho }}_3,\hfill \\ {\tilde{\rho }}_3{\tilde{\lambda }}_r^2+{\tilde{\rho }}_3{d}_i{\tilde{\lambda }}_r{x}_4^{*}-{\tilde{\rho }}_3{d}_i{d}_a{\left(x_4^{*}\right)}^2-\tilde{\delta }\left({\tilde{\rho }}_1{d}_a+{\tilde{\rho }}_3{d}_i\right){\left(x_4^{*}\right)}^3\hfill \\ -{\tilde{\rho }}_1{\tilde{\delta }}^2{\left(x_4^{*}\right)}^4>0,\hfill \end{array}\hfill \end{array}$

where K = 1 for $S_2^{*}$, and $K=(\tilde{\beta }-d_F)/(1+\tilde{\beta })$ for $S_3^{*}$. Biologically, the steady state $S_2^{*}$ represents the death of organ cells, while $S_3^{*}$ corresponds to an autoimmune regime.

The last steady state $S_4^{*}$ has all of its components positive and corresponds to the state of chronic infection.

At the next order, stochastic fluctuations are determined by linear stochastic processes, hence, this is known as a linear noise approximation (van Kampen, 1981; Wallace et al., 2012). The dynamics of these fluctuations is described by the following linear Fokker-Planck equation

$\begin{array}{ll}\frac{\partial \Pi \left(\zeta ,t\right)}{\partial t}=-{\displaystyle \sum _{i,j}}{A}_{ij}\frac{\partial }{\partial {\zeta }_i}\left({\zeta }_j\Pi \right)+\frac{1}{2}{\displaystyle \sum _{i,j}}{B}_{ij}\frac{{\partial }^2\Pi }{\partial {\zeta }_i\partial {\zeta }_j},& (10)\end{array}$

where A is the Jacobian matrix of system (8)

$\begin{array}{l}A=\left(\begin{array}{ccccccc}{b}_1+2{\tilde{b}}_2{x}_1-d_1-2{\tilde{d}}_2{x}_1-{\tilde{\mu }}_a{x}_6-\tilde{\beta }{x}_2& -\tilde{\beta }{x}_1& 0& 0& 0& -{\tilde{\mu }}_a{x}_1& 0\\ \tilde{\beta }{x}_2& \tilde{\beta }{x}_1-d_F-{\tilde{\mu }}_F{x}_5-{\tilde{\mu }}_a{x}_6& 0& 0& -{\tilde{\mu }}_F{x}_2& -{\tilde{\mu }}_a{x}_2& 0\\ 0& -\tilde{\alpha }{x}_3& -d_{in}-\tilde{\alpha }{x}_2& 0& 0& 0& 0\\ 0& p_1\tilde{\alpha }{x}_3& p_1\tilde{\alpha }{x}_2& {\tilde{\rho }}_1{x}_7-d_r& 0& 0& {\tilde{\rho }}_1{x}_4\\ 0& p_2\tilde{\alpha }{x}_3& p_2\tilde{\alpha }{x}_2& 0& {\tilde{\rho }}_2{x}_7-d_n& 0& {\tilde{\rho }}_2{x}_5\\ 0& \left(1-p_1-p_2\right)\tilde{\alpha }{x}_3& \left(1-p_1-p_2\right)\tilde{\alpha }{x}_2& -\tilde{\delta }{x}_6& 0& {\tilde{\rho }}_3{x}_7-d_a-\tilde{\delta }{x}_4& {\tilde{\rho }}_3{x}_6\\ 0& 0& 0& 0& {\sigma }_1& {\sigma }_2& -d_i\end{array}\right),\end{array}$

and B is a 7 × 7 symmetric matrix given by

$\begin{array}{l}{B}_{ij}=\{\begin{array}{ll}{b}_1{x}_1+{\tilde{b}}_2{x}_1^2+d_1{x}_1+{\tilde{d}}_2{x}_1^2+\tilde{\beta }{x}_1{x}_2+{\tilde{\mu }}_a{x}_1{x}_6,\hfill & \text{if }(i,j)=(1,1),\hfill \\ \tilde{\beta }{x}_1{x}_2+d_F{x}_2+{\tilde{\mu }}_F{x}_2{x}_5+{\tilde{\mu }}_a{x}_2{x}_6,\hfill & \text{if }(i,j)=(2,2),\hfill \\ {\tilde{\lambda }}_{in}+d_{in}{x}_3+\tilde{\alpha }{x}_2{x}_3,\hfill & \text{if }(i,j)=(3,3),\hfill \\ {\tilde{\lambda }}_r+d_r{x}_4+p_1\tilde{\alpha }{x}_2{x}_3+{\tilde{\rho }}_1{x}_4{x}_7,\hfill & \text{if }(i,j)=(4,4),\hfill \\ p_2\tilde{\alpha }{x}_2{x}_3+d_n{x}_5+{\tilde{\rho }}_2{x}_5{x}_7,\hfill & \text{if }(i,j)=(5,5),\hfill \\ (1-p_1-p_2)\tilde{\alpha }{x}_2{x}_3+d_a{x}_6+\tilde{\delta }{x}_4{x}_6+{\tilde{\rho }}_3{x}_6{x}_7,\hfill & \text{if }(i,j)=(6,6),\hfill \\ {\sigma }_1{x}_5+{\sigma }_2{x}_6+d_i{x}_7,\hfill & \text{if }(i,j)=(7,7),\hfill \\ -\tilde{\beta }{x}_1{x}_2,\hfill & \text{if }(i,j)=(1,2)\text{ or }(2,1),\hfill \\ -p_1\tilde{\alpha }{x}_2{x}_3,\hfill & \text{if }(i,j)=(3,4)\text{ or }(4,3),\hfill \\ -p_2\tilde{\alpha }{x}_2{x}_3,\hfill & \text{if }(i,j)=(3,5)\text{ or }(5,3),\hfill \\ -(1-p_1-p_2)\tilde{\alpha }{x}_2{x}_3,\hfill & \text{if }(i,j)=(3,6)\text{ or }(6,3),\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\end{array}$

Since the Fokker-Planck Equation (10) is linear, the probability density Π(ζ, t) is Gaussian, and hence, just the first two moments are enough to characterize it (Hayot and Jayaprakash, 2004; Pahle et al., 2012). Due to the way the system size expansion was introduced in (Equation 5), the mean values of fluctuations for all variables are zero, i.e., 〈ζ_i(t)〉 = 0 for all 1 ≤ i ≤ 7, while the covariance matrix Ξ with Ξ_ij = 〈ζ_i(t)ζ_j(t)〉 − 〈ζ_i(t)〉〈ζ_j(t)〉 = 〈ζ_i(t)ζ_j(t)〉 satisfies the following equation (van Kampen, 1981; Pahle et al., 2012)

$\begin{array}{ll}{\partial }_t\Xi =A\Xi +\Xi A^T+B,& (11)\end{array}$

where A^T is the transpose of A.

We are mainly interested in the dynamics of fluctuations when the oscillations of the deterministic model have died out, and the system is in a stationary state, i.e., the fluctuations take place around the steady states (Black et al., 2009). If the model (8) tends to a steady state as t → ∞, then in the equation (10) one can substitute the values of x_i's with the corresponding constant components of that steady state to study the fluctuations around it, as described by the linear Fokker-Planck equation. At any steady state, the covariance matrix Ξ is independent of time, and the fluctuations are described by a Gaussian distribution with the zero mean and the stationary covariance satisfying the equation

$\begin{array}{l}A\Xi +\Xi A^T+B=0.\end{array}$

In order to be able to relate the results of this analysis to simulations, it is convenient to express the covariance matrix in terms of actual numbers of cells in each compartment, rather than deviations from stationary values. To this end, we instead use the covariance matrix C defined as C_ij = 〈(n_i − 〈n_i〉)(n_j − 〈n_j〉)〉, which, in light of the relation C_ij = ΩΞ_ij, satisfies the following Lyapunov equation (Pahle et al., 2012)

$\begin{array}{ll}AC+C{A}^T+\Omega B=0.& (12)\end{array}$

This equation can be solved numerically for each of the stable steady states to determine the variance of fluctuations around that steady state depending on system parameters.

3. Results

To simulate the dynamics of the model, we solve the system Equation (4) numerically using the Euler-Maruyama method with parameter values given in Table 2, and Ω = 1000. The initial condition is chosen to be of the form

$\begin{array}{l}\left(x_1\left(0\right),x_2\left(0\right),x_3\left(0\right),x_4\left(0\right),x_5\left(0\right),x_6\left(0\right),x_7\left(0\right)\right)\\ =\left(18,2,7.2,6.3,0,0,0\right),& (13)\end{array}$

which corresponds to a small number of host cells being initially infected.

TABLE 2

Table 2. Table of parameters.

Figure 2 shows the results of 20,000 simulations with the initial condition (13) and σ₂ = 1. In the deterministic model (8), for σ₂ = 1 both steady states $S_1^{*}$ (disease-free) and $S_3^{*}$ (autoimmune state) are stable, but with the initial condition (13) the system is in the basin of attraction of $S_3^{*}$. In the stochastic model, the majority of trajectories also enter the attraction region of $S_3^{*}$, but a small proportion of them went into the basin of attraction of $S_1^{*}$. This figure illustrates a single stochastic path around $S_1^{*}$, and a single stochastic path around $S_3^{*}$, together with the deterministic trajectory. These individual solutions indicate that whilst deterministically, the system exhibits decaying oscillations around $S_3^{*}$, the same behavior is observed in the stochastic simulations only upon taking an average of a very large number of simulations. At the same time, individual realizations exhibit sustained stochastic oscillations in a manner similar to that observed in models of stochastic amplification in epidemics (Alonso et al., 2007; Kuske et al., 2007). Figure 2 also illustrates the size of areas of one standard deviation from the mean for trajectories in the basins of attraction $S_1^{*}$ and $S_3^{*}$, in which individual stochastic trajectories may exhibit stochastic oscillations (Conway and Coombs, 2011; Reynolds et al., 2012).

FIGURE 2

Figure 2. Numerical simulation of the model (4) with parameter values from Table 2, σ₂ = 1, and the initial condition (13). Red curves are two sample paths that have entered the basins of attraction of $S_1^{*}$ or $S_3^{*}$, black curve is the deterministic trajectory from (1), and the shaded areas indicate the regions of one standard deviation from the mean.

Figures 3A,B show temporal evolution of the probability distribution in the case of bi-stability between the steady states $S_1^{*}$ and $S_3^{*}$, as illustrated in Figure 2. They indicate that after some initial transient, the system reaches a stationary bimodal normal distribution. The width of the probability distribution around each stable steady state, as described by its variance or standard deviation, gives the size of fluctuations around this steady state observed in individual stochastic realizations, as is shown in Figure 2. Similar behavior has been observed in stochastic realizations of other deterministic models with bi-stability (Earnest et al., 2013; Bruna et al., 2014; Hufton et al., 2016). For the parameter values given in Table 2, the deterministic system exhibits a bi-stability between $S_1^{*}$ and $S_2^{*}$, and with the initial condition

$\begin{array}{l}\left(x_1\left(0\right),x_2\left(0\right),x_3\left(0\right),x_4\left(0\right),x_5\left(0\right),x_6\left(0\right),x_7\left(0\right)\right)\\ =\left(18,9,7.2,6.3,0,0,0\right),& (14)\end{array}$

it is in the basin of attraction of $S_2^{*}$. Due to stochasticity, the stationary probability distribution in this case is also bimodal, with the majority of solutions being distributed around $S_2^{*}$, and a very small number being centered around $S_1^{*}$, as can be seen in Figures 3C,D. Increasing the system size Ω is known to result in the bimodal distribution becoming unimodal due to the size of fluctuations scaling as Ω^−1/2, which results in a reduced variability in trajectories (Black and McKane, 2012; Hufton et al., 2016), and the same conclusion holds for the system (4).

FIGURE 3

Figure 3. Probability distribution of solutions out of 20,000 simulations. (A,B) with parameters from Table 2, but σ₂ = 1 and the initial condition (13). (C,D) with parameters from Table 2 and the initial condition (14). In (A,C), the probability histogram is fit to a bimodal normal distribution at different times. (B,D) illustrate stationary joint probability histograms.

To gain better insights into the role of initial conditions, in Figure 4 we fix all parameter values, and vary initial numbers of infected cells and regulatory T cells. For the parameter combination illustrated in Figure 4A, the deterministic model exhibits a bi-stability between a stable disease-free steady state $S_1^{*}$ and a periodic oscillation around the state $S_3^{*}$, which biologically corresponds to an autoimmune regime. In the deterministic case, the black boundary provides a clear separation of the basins of attraction of these two dynamical states, in a manner similar to that investigated recently in the context of within-cell dynamics of RNA interference (Neofytou et al., 2017). For stochastic simulations, the color indicates the probability of the solution going to a disease-free state $S_1^{*}$, and it shows that even in the case where deterministically the system is in the basin of attraction of one of the states, there is a non-zero probability that it will actually end up at another state, with this probability varying smoothly across the deterministic basin boundary. This figure suggests that if the initial number of infected cells is sufficiently small, or if the number of regulatory T cells is sufficiently large, the system tends to clear the infection and approach the disease-free state. On the contrary, for higher numbers of infected cells and lower numbers of regulatory cells, autoimmune regime appears to be a more likely outcome. Qualitatively similar behavior is observed for another combination of parameters illustrated in Figure 4B, in which case the deterministic system has a bi-stability between a disease-free steady state $S_1^{*}$, and a state $S_2^{*}$ which represents the death of host cells.

FIGURE 4

Figure 4. Probability of solution entering and staying in the basin of attraction of the disease-free steady state $S_1^{*}$ in the bi-stability regime with A(0) = 18, 000 and T_in(0) = 7, 200. Black curves are the boundaries between different basins of attraction in the deterministic model. (A) With parameter values from Table 2, ${\tilde{\lambda }}_r=45$ and ${\tilde{\mu }}_a=10/9$, in the region below the black curve, the deterministic model exhibits a periodic solution around $S_3^{*}$, and above this curve is the deterministic basin of attraction of $S_1^{*}$. (B) With parameter values from Table 2, area below the black curve is the basin of attraction of $S_2^{*}$, and above it is again the basin of attraction of $S_1^{*}$.

In order to understand how biological parameters affect the size of fluctuations around steady states, in Figure 5 we explore several parameter planes by first identifying parameter regions where the deterministic system has a stable steady state $S_3^{*}$, and then for each combination of parameters inside these regions, we use the Bartels-Stewart method (Bartels and Stewart, 1972; Hammarling, 1982) to numerically solve the Lyapunov equation (12) and compute the variance in the number of regulatory T cell when the deterministic model is at the steady state $S_3^{*}$. The value of variance gives the square of the magnitude of oscillations observed in individual stochastic realizations. One should note that getting closer to the deterministic boundary of stability of $S_3^{*}$ increases the stochastic variance of fluctuations around this steady state. The reason for this is that closer parameters are to the deterministic stability boundary, the less stable is the steady state, hence the larger is the amplitude of stochastic oscillations around it. Moreover, the variance increases with the rate of production of IL-2 by autoreactive T cells and the rate at which regulatory T cells suppress autoreactive T cells; it decreases with the higher rate of production of regulatory T cells, and it appears to not depend on the rate at which autoreactive T cells destroy infected cells, or on the infection rate.

FIGURE 5

Figure 5. Variance of the number of regulatory T cells T_reg with parameter values from Table 2. Colored regions indicate areas in respective parameter planes in which the autoimmune steady state $S_3^{*}$ is deterministically stable.

4. Discussion

In this paper we have analyzed stochastic aspects of immune response against a viral infection with account for the populations of T cells with different activation thresholds, as well as cytokines mediating T cell activity. The CTMC model has provided an exact master equation, for which we applied a van Kampen's expansions to derive a linear Fokker-Planck equation that characterizes fluctuations around the deterministic solutions. We have also explored actual stochastic trajectories of the system by deriving an SDE model and solving it numerically.

One biologically important aspect we have looked at is the influence of stochasticity on the dynamics of the system in the case where deterministically it exhibits a bi-stability between either two steady states, or a steady state and a periodic solution. In such a situation, bi-stability in the deterministic version of the model translates in the stochastic case into a stationary bimodal distribution for the probability density. To obtain further insights into details of how stochasticity affects bi-stability, we have investigated how for the fixed parameter values time evolution of the system changes depending on the initial numbers of the regulatory T cells and infected cells.

Our analysis reinforces the need to distinguish mean dynamics from individuals realizations: where in the deterministic case the system can approach a stable steady state (which represents mean behavior of a very large number of simulations), individual realizations can exhibit sustained stochastic oscillations around that steady state, as we have seen in numerical simulations. Since in the clinical or laboratory setting one is usually dealing with single measurements of some specific biological quantities rather than their averaged values, the stochastic oscillations exhibited by our model may quite well explain observed variability in the measured levels of infection or T cell populations. To better understand the magnitude of stochastic fluctuations around the deterministic steady states, we have solved the Lyapunov equation, which has provided us with a quantitative information on the dependence of variance of fluctuations on system parameters.

There are several directions in which the work presented in this paper can be extended. In terms of fundamental immunology, the model can be made more realistic by including additional effects, such as the control of IL-2 secretion by regulatory T cells (Burroughs et al., 2006), or the memory T cells (Antia et al., 2005; Skapenko et al., 2005). Whilst we have used numerical simulations to compute the probability of attraction to a given steady state in the case of bi-stability, one could approach the same problem theoretically from the perspective of computing extinction probability within the framework of the CTMC model (Yuan and Allen, 2011; Mandal et al., 2014). The van Kampen's system size expansion could yield an expression for the power spectrum, which allows one to compute the peak frequency and amplification (McKane and Newman, 2005; Alonso et al., 2007; Black et al., 2009; Black and McKane, 2010). From a practical perspective, future work could focus on validating theoretical results presented in this paper using experimental measurements of the progress of autoimmune disease in animal hosts, with experimental autoimmune uveoretinitis (EAU), an autoimmune inflammation in the eyes, being one interesting possibility. In one such recent experiment, all animals were genetically identical C57BL/6 mice, but once the EAU was induced in them through inoculation, the autoimmune disease then progressed at slightly different rates (Boldison et al., 2015; Boldison and Nicholson, unpublished) and the measured variability in the numbers of infected cells and T cell responses could be compared to theoretical estimates of the variance as predicted by our model. From a clinical perspective, comparison of variance in the measured populations of different cells with the model conclusions will facilitate an efficient parameter identification and provide a set of prognostic criteria for the progress of autoimmunity, which can be used for risk stratification and assessment of patients with autoimmune disease.

Author Contributions

YK and KB designed the model; FF, SK, and AR performed the analysis and simulations, and produced the figures. All authors drafted and edited the manuscript.

Funding

FF has been sponsored by the Chancellor's Award from Sussex University. AR acknowledges funding for her PhD studies from the Engineering and Physical Sciences Research Council (EPSRC) through iCASE Award EP/N509784/1.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.00045/full#supplementary-material

References

Alexander, H. K., and Wahl, L. M. (2011). Self-tolerance and autoimmunity in a regulatory T cell model. Bull. Math. Biol. 73, 33–71. doi: 10.1007/s11538-010-9519-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, E. J. (2007). Modeling with Itô tochastic Differential Equations. Berlin: Springer.

Google Scholar

Allen, E. J., Allen, L. J. S., Arciniega, A., and Greenwood, P. E. (2008). Construction of equivalent stochastic differential equation models. Stoch. Anal. Appl. 26, 274–297. doi: 10.1080/07362990701857129

CrossRef Full Text | Google Scholar

Allen, L. J. S. (2011). An Introduction to Stochastic Processes with Applications to Biology. Boca Raton, FL: CRC Press.

Google Scholar

Alonso, D., McKane, A. J., and Pascual, M. (2007). Stochastic amplification in epidemics. J. R. Soc. Interface 4, 575–582. doi: 10.1098/rsif.2006.0192

PubMed Abstract | CrossRef Full Text | Google Scholar

Altan-Bonnet, G., and Germain, R. N. (2005). Modeling T cell antigen discrimination based on feedback control of digital ERK responses. PLoS Biol. 3:e356. doi: 10.1371/journal.pbio.0030356

PubMed Abstract | CrossRef Full Text | Google Scholar

Antia, R., Ganusov, V. V., and Ahmed, R. (2005). The role of models in understanding CD8⁺ T-cell memory. Nat. Rev. Immunol. 5, 101–111. doi: 10.1038/nri1550

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartels, R. H., and Stewart, G. W. (1972). Solution of the matrix equation AX+XB=C. Comm. ACM 15, 820–826. doi: 10.1145/361573.361582

CrossRef Full Text | Google Scholar

Ben Ezra, D., and Forrester, J. V. (1995). Fundal white dots: the spectrum of a similar pathological process. Br. J. Ophthalmol. 79, 856–860. doi: 10.1136/bjo.79.9.856

PubMed Abstract | CrossRef Full Text | Google Scholar

Bitmansour, A. D., Douek, D. C., Maino, V. C., and Picker, L. J. (2002). Direct ex vivo analysis of human CD4⁺ memory T cell activation requirements at the single clonotype level. J. Immunol. 169, 1207–1218. doi: 10.4049/jimmunol.169.3.1207

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, A. J., and McKane, A. J. (2010). Stochastic amplification in an epidemic model with seasonal forcing. J. Theor. Biol. 267, 85–94. doi: 10.1016/j.jtbi.2010.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, A. J., and McKane, A. J. (2012). Stochastic formulation of ecological models and their applications. Trends Ecol. Evol. 27, 337–345. doi: 10.1016/j.tree.2012.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, A. J., McKane, A. J., Nunes, A., and Parisi, A. (2009). Stochastic fluctuations in the susceptible-infective-recovered model with distributed infectious periods. Phys. Rev. E 80:021922. doi: 10.1103/PhysRevE.80.021922

PubMed Abstract | CrossRef Full Text | Google Scholar

Blattman, J. N., Sourdive, D.J.D., Murali-Krishna, K., Ahmed, R., and Altman, J. D. (2000). Evolution of the T cell repertoire during primary, memory, and recall response to viral infection. J. Immunol. 165, 6081–6090. doi: 10.4049/jimmunol.165.11.6081

CrossRef Full Text | Google Scholar

Blyuss, K. B., and Nicholson, L. B. (2012). The role of tunable activation thresholds in the dynamics of autoimmunity. J. Theor. Biol. 308, 45–55. doi: 10.1016/j.jtbi.2012.05.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Blyuss, K. B., and Nicholson, L. B. (2015). Understanding the roles of activation threshold and infections in the dynamics of autoimmune disease. J. Theor. Biol. 375, 13–20. doi: 10.1016/j.jtbi.2014.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Boldison, J., Khera, T. K., Copland, D. A., Stimpson, M. L., Crawford, G. L., Dick, A. D., et al. (2015). A novel pathogenic RBP-3 peptide reveals epitope spreading in persistent experimental autoimmune uveoretinitis. Immunology 146, 301–311. doi: 10.1111/imm.12503

PubMed Abstract | CrossRef Full Text | Google Scholar

Borghans, J. A. M., and De Boer, R. J. (1995). A minimal model for T-cell vaccination. Proc. R. Soc. B 259, 173–178. doi: 10.1098/rspb.1995.0025

PubMed Abstract | CrossRef Full Text | Google Scholar

Borghans, J. A. M., De Boer, R. J., Sercarz, E., and Kumar, V. (1998). T cell vaccination in experimental autoimmune encephalomyelitis: a mathematical model. J. Immunol. 161, 1087–1093.

PubMed Abstract | Google Scholar

Brauer, F., van den Driessche, P., and Wu, J. (2008). Mathematical Epidemiology. Berlin: Springer.

PubMed Abstract | Google Scholar

Braumann, C. A. (2007). Itô versus Stratonovich calculus in random population growth. Math. Biosci. 206, 81–107. doi: 10.1016/j.mbs.2004.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruna, M., Chapman, S. J., and Smith, M. J. (2014). Model reduction for slow-fast stochastic systems with metastable behaviour. J. Chem. Phys. 140:174107. doi: 10.1063/1.4871694

PubMed Abstract | CrossRef Full Text | Google Scholar

Burroughs, N. J., Ferreira, M., Oliveira, B. M. P. M., and Pinto, A. A. (2011a). A transcritical bifurcation in an immune response model. J. Diff. Eqns. Appl. 17, 1101–1106. doi: 10.1080/10236190903095291

CrossRef Full Text | Google Scholar

Burroughs, N. J., Oliveira, B. M. P. M., and Pinto, A. A. (2006). Regulatory T cell adjustment of quorum growth thresholds and the control of local immune responses. J. Theor. Biol. 241, 134–141. doi: 10.1016/j.jtbi.2005.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Burroughs, N. J., Oliveira, B. M. P. M., and Pinto, A. A. (2011b). Autoimmunity arising from bystander proliferation of T cells in an immune response model. Math. Comput. Mod. 53, 1389–1393. doi: 10.1016/j.mcm.2010.01.020

CrossRef Full Text | Google Scholar

Carneiro, J., Paixão, T., and Milutinovic, D. (2005). Immunological self-tolerance: lessons from mathematical modeling. J. Comput. Appl. Math. 184, 77–100. doi: 10.1016/j.cam.2004.10.025

CrossRef Full Text | Google Scholar

Chao, D. L., Davenport, M. P., Forrest, S., and Perelson, A. S. (2004). A stochastic model of cytotoxic T cell responses. J. Theor. Biol. 228, 227–240. doi: 10.1016/j.jtbi.2003.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Conway, J. M., and Coombs, D. (2011). A stochastic model of latently infected cell reactivation and viral blip generation in treated HIV patients. PLoS Comp. Biol. 7:e1002033. doi: 10.1371/journal.pcbi.1002033

PubMed Abstract | CrossRef Full Text | Google Scholar

Davies, T. F., Evered, D. C., Rees Smith, B., Yeo, P. P. B., Clark, F., and Hall, R. (1997). Value of thyroid-stimulating-antibody determination in predicting the short-term thyrotoxic relapse in graves' disease. Lancet 309, 1181–1182. doi: 10.1016/S0140-6736(77)92719-2

CrossRef Full Text

Deenick, E. K., Gett, A. V., and Hodgkin, P. D. (2003). Stochastic model of T cell proliferation: A calculus revealing IL-2 regulation of precursor frequencies, cell cycle time, and survival. J. Immunol. 170, 4963–4972. doi: 10.4049/jimmunol.170.10.4963

PubMed Abstract | CrossRef Full Text | Google Scholar

Detours, V., and Perelson, A. S. (2000). The paradox of alloreactivity and self MHC restriction: Quantitative analysis and statistics. Proc. Natl. Acad. Sci. U.S.A. 97, 8479–8483. doi: 10.1073/pnas.97.15.8479

PubMed Abstract | CrossRef Full Text | Google Scholar

Earnest, T. M., Roberts, E., Assaf, M. K. D., and Luthey-Schulten, Z. (2013). DNA looping increases the range of bistability in a stochastic model of the lac genetic switch. Phys. Biol. 10:26002. doi: 10.1088/1478-3975/10/2/026002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ercolini, A. M., and Miller, S. D. (2008). The role of infections in autoimmune disease. Clin. Exp. Immunol. 155, 1–15. doi: 10.1111/j.1365-2249.2008.03834.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujinami, R. S. (2011). Can virus infections trigger autoimmune disease? J. Autoimmun. 16, 229–234. doi: 10.1006/jaut.2000.0484

PubMed Abstract | CrossRef Full Text | Google Scholar

Gardiner, C. W. (2004). Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences. Berlin: Springer-Verlag.

Google Scholar

Gillespie, D. T. (2000). The chemical Langevin equation. J. Chem. Phys. 113, 297–306. doi: 10.1063/1.481811

CrossRef Full Text | Google Scholar

Gillespie, D. T. (2002). The chemical Langevin and Fokker-Planck equations for the reversible isomerization reaction. J. Phys. Chem. A 106, 5063–5071. doi: 10.1021/jp0128832

CrossRef Full Text | Google Scholar

Grossman, Z., and Paul, W. E. (1992). Adaptive cellular interactions in the immune system: the tunable activation threshold and the significance of subthreshold responses. Proc. Natl. Acad. Sci. U.S.A. 89, 10365–10369. doi: 10.1073/pnas.89.21.10365

PubMed Abstract | CrossRef Full Text | Google Scholar

Grossman, Z., and Paul, W. E. (2000). Self-tolerance: context dependent tuning of T cell antigen recognition. Sem. Immunol. 12, 197–203. doi: 10.1006/smim.2000.0232

PubMed Abstract | CrossRef Full Text | Google Scholar

Grossman, Z., and Singer, A. (1996). Tuning of activation thresholds explains flexibility in the selection and development of T cells in the thymus. Proc. Natl. Acad. Sci. U.S.A. 93, 14747–14752. doi: 10.1073/pnas.93.25.14747

PubMed Abstract | CrossRef Full Text | Google Scholar

Hammarling, S. J. (1982). Numerical solution of the stable, non-negative definite lyapunov equation. IMA J. Num. Anal. 2, 303–323. doi: 10.1093/imanum/2.3.303

CrossRef Full Text | Google Scholar

Hayot, F., and Jayaprakash, C. (2004). The linear noise approximation for molecular fluctuations within cells. Phys. Biol. 1, 205–210. doi: 10.1088/1478-3967/1/4/002

PubMed Abstract | CrossRef Full Text | Google Scholar

Hufton, P. G., Lin, Y. T., Galla, T., and McKane, A. J. (2016). Intrinsic noise in systems with switching environments. Phys. Rev. E 93:52119. doi: 10.1103/PhysRevE.93.052119

PubMed Abstract | CrossRef Full Text | Google Scholar

Iwami, S., Takeuchi, Y., Iwamoto, K., Naruo, Y., and Yasukawa, M. (2009). A mathematical design of vector vaccine against autoimmune disease. J. Theor. Biol. 256, 382–392. doi: 10.1016/j.jtbi.2008.09.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Iwami, S., Takeuchi, Y., Miura, Y., Sasaki, T., and Kajiwara, T. (2007). Dynamical properties of autoimmune disease models: tolerance, flare-up, dormancy. J. Theor. Biol. 246, 646–659. doi: 10.1016/j.jtbi.2007.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerr, E. C., Copland, D. A., Dick, A. D., and Nicholson, L. B. (2008). The dynamics of leukocyte infiltration in experimental autoimmune uveoretinitis. Eye Res. 27, 527–535. doi: 10.1016/j.preteyeres.2008.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuske, R., Gordillo, L. F., and Greenwood, P. E. (2007). Sustained oscillations via coherence resonance in SIR. J. Theo. Biol. 245, 459–469. doi: 10.1016/j.jtbi.2006.10.029

PubMed Abstract | CrossRef Full Text | Google Scholar

León, K., Faro, J., Lage, A., and Carneiro, J. (2004). Inverse correlation between the incidences of autoimmune disease and infection predicted by a model of T cell mediated tolerance. J. Autoimmun. 22, 31–42. doi: 10.1016/j.jaut.2003.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

León, K., Lage, A., and Carneiro, J. (2003). Tolerance and immunity in a mathematical model of T-cell mediated suppression. J. Theor. Biol. 225, 107–126. doi: 10.1016/S0022-5193(03)00226-1

PubMed Abstract | CrossRef Full Text | Google Scholar

León, K., Perez, R., Lage, A., and Carneiro, J. (2000). Modelling T-cell-mediated suppression dependent on interactions in multicellular conjugates. J. Theor. Biol. 207, 231–254. doi: 10.1006/jtbi.2000.2169

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandal, P. S., Allen, L. J. S., and Banerjee, M. (2014). Stochastic modeling of phytoplankton allelopathy. Appl. Math. Model. 38, 1583–1596. doi: 10.1016/j.apm.2013.08.031

CrossRef Full Text | Google Scholar

McKane, A. J., and Newman, T. J. (2005). Predator-prey cycles from resonant amplification of demographic stochasticity. Phys. Rev. Lett. 94:218102. doi: 10.1103/PhysRevLett.94.218102

PubMed Abstract | CrossRef Full Text | Google Scholar

Neofytou, G., Kyrychko, Y. N., and Blyuss, K. B. (2017). Time-delayed model of RNA interference. Ecol. Complex. 30, 11–25. doi: 10.1016/j.ecocom.2016.12.003

CrossRef Full Text | Google Scholar

Nicholson, L. B., Anderson, A. C., and Kuchroo, V. K. (2000). Tuning T cell activation threshold and effector function with cross-reactive peptide ligands. Int. Immunol. 12, 205–213. doi: 10.1093/intimm/12.2.205

PubMed Abstract | CrossRef Full Text | Google Scholar

Noest, A. J. (2000). Designing lymphocyte functional structure for optimal signal detection: voilá. J. Theor. Biol. 207, 195–216. doi: 10.1006/jtbi.2000.2164

PubMed Abstract | CrossRef Full Text | Google Scholar

Nylander, A., and Hafler, D. A. (2012). Multiple sclerosis. J. Clin. Invest. 122, 1180–1188. doi: 10.1172/JCI58649

PubMed Abstract | CrossRef Full Text

Øksendal, B. (2000). Stochastic Differential Equations: An Introduction with Applications. Berlin: Springer.

Pahle, J., Challenger, J. D., Mendes, P., and McKane, A. J. (2012). Biochemical fluctuations, optimisation and the linear noise approximation. BMC Syst. Biol. 6:86. doi: 10.1186/1752-0509-6-86

PubMed Abstract | CrossRef Full Text | Google Scholar

Perelson, A. S., and Weisbuch, G. (1997). Immunology for physicists. Rev. Mod. Phys. 69, 1219–1267. doi: 10.1103/RevModPhys.69.1219

CrossRef Full Text | Google Scholar

Prat, E., and Martin, R. (2002). The immunopathogenesis of multiple sclerosis. J. Rehabil. Res. Dev. 39, 187–199.

PubMed Abstract | Google Scholar

Reynolds, J., Coles, M., Lythe, G., and Molina-París, C. (2012). Deterministic and stochastic naive T cell population dynamics: symmetric and asymmetric cell division. Dyn. Syst. 27, 75–103. doi: 10.1080/14689367.2011.645447

CrossRef Full Text | Google Scholar

Römer, P. S., Berr, S., Avota, E., Na, S.-Y., Battaglia, M., ten Berge, I., et al. (2011). Preculture of PBMC at high cell density increases sensitivity of T-cell responses, revealing cytokine release by CD28 superagonist TGN1412. Blood 118, 6772–6782. doi: 10.1182/blood-2010-12-319780

CrossRef Full Text | Google Scholar

Root-Bernstein, R. (2015). Theories and modeling of autoimmunity. J. Theor. Biol. 375, 1–124. doi: 10.1016/j.jtbi.2015.04.003

CrossRef Full Text

Santamaria, P. (2010). The long and winding road to understanding and conquering type 1 diabetes. Immunity 32, 437–445. doi: 10.1016/j.immuni.2010.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Scherer, A., Noest, A., and deBoer, R. J. (2004). Activation-threshold tuning in an affinity model for the T-cell repertoire. Proc. R. Soc. B 271, 609–616. doi: 10.1098/rspb.2003.2653

PubMed Abstract | CrossRef Full Text | Google Scholar

Segel, L. A., Jäger, E., Elias, D., and Cohen, I. R. (1995). A quantitative model of autoimmune disease and T-cell vaccination: does more mean less? Immunol. Today 16, 80–84. doi: 10.1016/0167-5699(95)80093-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Skapenko, A., Leipe, J., Lipsky, P. E., and Schulze-Koops, H. (2005). The role of the T cell in autoimmune inflammation. Arthr. Res. Ther. 7(Suppl. 2), S4–S14. doi: 10.1186/ar1703

PubMed Abstract | CrossRef Full Text | Google Scholar

Stefanova, I., Dorfman, J. R., and Germain, R. N. (2002). Self-recognition promotes the foreign antigen sensitivity of naive T lymphocytes. Nature 420, 429–434. doi: 10.1038/nature01146

PubMed Abstract | CrossRef Full Text | Google Scholar

Stirk, E. R., Lythe, G., van den Berg, H. A., Hurst, G. A. D., and Molina-París, C. (2010a). The limiting conditional probability distribution in a stochastic model of T cell repertoire maintenance. Math. Biosci. 224, 74–86. doi: 10.1016/j.mbs.2009.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Stirk, E. R., Lythe, G., van den Berg, H. A., and Molina-París, C. (2010b). Stochastic competitive exclusion in the maintenance of the naïve T cell repertoire. J. Theor. Biol. 265, 396–410. doi: 10.1016/j.jtbi.2010.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Berg, H. A., and Rand, D. A. (2004). Dynamics of T cell activation threshold tuning. J. Theor. Biol. 228, 397–416. doi: 10.1016/j.jtbi.2004.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

van Kampen, N. G. (1981). Stochastic Processes in Physics and Chemistry. Amsterdam: Elsevier.

Google Scholar

von Herrath, M. G., and Oldstone, M. B. A. (1996). Virus-induced autoimmune disease. Curr. Opin. Immunol. 8, 878–885. doi: 10.1016/S0952-7915(96)80019-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wallace, E. W. J., Gillespie, D. T., Sanft, K. R., and Petzold, L. R. (2012). Linear noise approximation is valid over limited times for any chemical system that is sufficiently large. IET Syst. Biol. 6, 102–115. doi: 10.1049/iet-syb.2011.0038

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, S. D., Dittel, B. N., Hardardottir, F., and Janeway Jr, C. A. (1996). Experimental autoimmune encephalomyelitis induction in genetically B cell-deficient mice. J. Exp. Med. 184, 2271–2278. doi: 10.1084/jem.184.6.2271

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H.-J., Ivanov, I. I., Darceet, J., Hattori, K., Shima, T., Umesaki, Y., et al. (2010). Gut-residing segmented filamentous bacteria drive autoimmune arthritis via T helper 17 cells. Immunity 32, 815–827. doi: 10.1016/j.immuni.2010.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, Y., and Allen, L. J. S. (2011). Stochastic models for virus and immune system dynamics. Math. Biosci. 234, 84–94. doi: 10.1016/j.mbs.2011.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar