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.


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. (2000León et al. ( , 2003León et al. ( , 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. (2007Iwami et al. ( , 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 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. 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.

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

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 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 1 of them will develop into regulatory T cells, a proportion p 2 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 σ 1 and σ 2 , and IL-2 in turn facilitates proliferation of regular, normal and autoreactive T cells at rates ρ 1 , ρ 2 , and ρ 3 , 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. 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 1 (t), . . . , X 7 (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 X 0 = (X 1 (0), . . . , X 7 (0)) = (n 10 , n 20 , n 30 , n 40 , n 50 , n 60 , n 70 ).
Here, b 1 n 1 + b 2 n 2 1 and d 1 n 1 + d 2 n 2 1 are natural birth and death rates for uninfected cells with b 1 − d 1 = r and d 2 − b 2 = r/N (Allen, 2011).
The probabilities P(n, t) satisfy the following master equation where the operators ε ± i are defined as follows, ε ± i f (n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , t) = f (n 1 , ..., n i ± 1, ..., n 7 , t), 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 highdimensional 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(Gillespie, , 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(Allen, , 2011Braumann, 2007).

Stochastic Differential Equation Model
To derive Itô SDE model, let Y(t) = (Y 1 (t), Y 2 (t), Y 3 (t), Y 4 (t), Y 5 (t), Y 6 (t), Y 7 (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 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 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., To derive Itô SDE model, we need to find a diffusion matrix H defined according to HH T = .

Now if we consider
and W(t) = [W 1 (t), W 2 (t), ..., W 11 (t)] 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 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).

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 1 /2 as follows, where x i (t) and ζ i (t) are two continuous variables, and x i (t) = E[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 + 1 /2 ζ , t , which implies To expand the master equation (3) in a power series in − 1 /2 , we use the following expansion for the step operators Substituting expressions (6, 7) into the master equation (see Supplementary Material for details) and collecting terms of order 1 /2 yields the following deterministic model for macroscopic behavior 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 and it is stable if d F > β. The second and third steady states can be found as where x * 4 satisfies the following quadratic equation These steady states are stable, provided where K = 1 for S * 2 , and K = β − d F / 1 + β 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 where A is the Jacobian matrix of system (8) Frontiers in Physiology | www.frontiersin.org and B is a 7 × 7 symmetric matrix given by if (i, j) = (3, 6) or (6, 3), 0, otherwise.
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 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 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) AC + CA T + B = 0.
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.

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 (x 1 (0), x 2 (0), x 3 (0), x 4 (0), x 5 (0), x 6 (0), x 7 (0)) = (18, 2, 7.2, 6.3, 0, 0, 0), which corresponds to a small number of host cells being initially infected. Figure 2 shows the results of 20,000 simulations with the initial condition (13) and σ 2 = 1. In the deterministic model (8), for σ 2 = 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). 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  Table 2, σ 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. 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 (x 1 (0), x 2 (0), x 3 (0), x 4 (0), x 5 (0), x 6 (0), x 7 (0)) = (18, 9, 7.2, 6.3, 0, 0, 0) 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).
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.
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.

DISCUSSION
In this paper we have analyzed stochastic aspects of immune response against a viral infection with account for the  Table 2, but σ 2 = 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. 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  Table 2, λ r = 45 and µ 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 .
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. 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.