Stochastic effects in autoimmune dynamics
F. Fatehi, S.N. Kyrychko, A. Ross, Y.N. Kyrychko, K.B. Blyuss
SStochastic effects in autoimmune dynamics
F. Fatehi, S.N. Kyrychko, A. Ross, Y.N. Kyrychko, K.B. Blyuss ∗ Department of Mathematics, University of Sussex, Falmer,Brighton, BN1 9QH, United KingdomNovember 9, 2018
Abstract
Among various possible causes of autoimmune disease, an important role is played by infectionsthat can result in a breakdown of immune tolerance, primarily through the mechanism of “molecularmimicry”. In this paper we propose and analyse a stochastic model of immune response to a viralinfection and subsequent autoimmunity, with account for the populations of T cells with differentactivation thresholds, regulatory T cells, and cytokines. We show analytically and numericallyhow 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 providea possible explanation for experimentally observed variations in the progression of autoimmunedisease. Computations of the variance of stochastic fluctuations provide practically importantinsights into how the size of these fluctuations depends on various biological parameters, and thisalso gives a headway for comparison with experimental data on variation in the observed numbersof T cells and organ cells affected by infection.
Breakdown of immune tolerance and the resulting autoimmune disease occur when the immune systemfails to distinguish the host’s own healthy cells from the cells affected by the infection, thus triggeringan immune response that also targets healthy cells. Autoimmune disease is usually focused in a specificorgan or part of the body, such as retina in the case of uveitis, central nervous system in multiplesclerosis, or pancreatic β -cells in insulin-dependent diabetes mellitus type-1 [1, 2, 3]. Whilst it is closeto impossible to pinpoint precise causes of autoimmunity in each individual case, it can usually beattributed to a number of factors, which can include the genetic predisposition, age, previous immunechallenges, exposure to pathogens etc. A number of distinct mechanisms have been identified for howan infection of the host with a pathogen can result in the subsequent onset of autoimmune disease, andthese include bystander activation [4] and molecular mimicry [5, 6], which is particularly importantin the context of autoimmunity caused by viral infections.Over the years, a number of mathematical models have investigated various origins and aspects ofimmune response, with an emphasis on the onset and the development of autoimmune disease. Someof the earlier models studied interactions between regulatory and effector T cells without looking atcauses of autoimmunity but instead focusing on T cell vaccination [7]. Borghans et al. [8, 9] looked intothis process in more detail and showed the onset of autoimmunity, which was defined as oscillationsin the number of autoreactive cells that exceeded a certain threshold. Le´on et al. [10, 11, 12] and ∗ Corresponding author. Email: [email protected] a r X i v : . [ q - b i o . T O ] M a y arneiro et al. [13] have analysed interactions between different T cells and their effect on regulationof immune response and control of autoimmunity. More recently, Iwami et al. [14, 15] considereda model of immune response to a viral infection, in which they explicitly included the dynamics ofa virus population. Although this model is able to demonstrate an emergence of autoimmunity, itfails to produce a regime of normal viral clearance. Alexander and Wahl [16] have focused on therole of professional antigen-presenting cells (APCs) and their interactions with regulatory and effectorcells for the purposes of controlling autoimmune response. Burroughs et al. [17, 18] have analysedthe emergence of autoimmunity through the mechanism of cytokine-mediated bystander activation.A special issue on “Theories and modelling of autoimmunity” provides an excellent overview of re-cent research in the area of mathematical modelling of various aspects of onset and development ofautoimmune disease [19].These are several different frameworks for modelling the role of T cells in controlling autoimmuneresponse. Alexander and Wahl [16] and Burroughs et al. [17, 18] have explicitly included a separatecompartment for regulatory T cells that are activated by autoantigens and suppress the activity ofautoreactive T cells. Another modelling approach is to consider the possibility of the same T cellsperforming different immune functions through having different or tunable activation thresholds, whichallows T cells to adjust their response to T cell antigen receptor stimulation by autoantigens. Thismethodology was originally proposed theoretically to study peripheral and central T cell activation [20,21, 22], and has been subsequently used to analyse differences in activation/response thresholds thatare dependent on the activation state of the T cell [23]. van den Berg and Rand [24] and Scherer etal. [25] have studied stochastic tuning of activation thresholds. Interestingly, the need for T cells tohave tunable activation can be shown to emerge from the fundamental principles of of signal detectiontheory [26]. A number of murine and human experiments have confirmed that activation of T cellscan indeed dynamically change during their circulation [27, 28, 29, 30], thus supporting the theorydeveloped in earlier papers.Since immune response is known to be a complex multi-factor process [31], a number of studies havelooked into various stochastic aspects of immune dynamics, such as T cell selection and proliferation.Deenick et al. [32] have analysed stochastic effects of interleukin-2 (IL-2) on T cell proliferationfrom precursors. Blattman et al. [33] have shown that repertoires of the CTL (cytotoxic T celllymphocyte) populations during primary response to a viral infection and in the memory pool aresimilar, thus providing further support to the theory of stochastic selection for the memory pool.Detours and Perelson [34] have explored the distribution of possible outcomes during T cell selectionwith account for variable affinity between T cell receptors and MHC-peptide complexes. Chao etal. [35] analysed a detailed stochastic model of T cell recruitment during immune response to aviral infection. Stirk et al. [36, 37] have developed a stochastic model for T cell repertoire andinvestigated the role of competitive exclusion between different clonotypes. Using the methodologyof continuous-time Markov processes, the authors computed extinction times, a limited multivariateprobability distribution, as well as the size of fluctuations around the deterministic steady states.Reynolds [38] have used a similar methodology to investigate an important question of asymmetriccell division and its impact on the extinction of different T cell populations and the expected lifetimesof na¨ıve T cell clones. With regards to modelling autoimmune dynamics, Alexander and Wahl [16]have studied the stochastic model of immune response with an emphasis on professional APCs to showthat the probability of developing a chronic autoimmune response increases with the initial exposureto self-antigen or autoreactive effector T cells. An important aspect of stochastic dynamics that hasto be accounted for in the models is the so-called stochastic amplification [39, 40], which denotes asituation where periodic solutions with decaying amplitudes in the deterministic model can result insustained stochastic periodic oscillations in individual realisations of the same model. This suggeststhat whilst on average the behaviour may show decaying-amplitude oscillations, individual realisations2epresented by stochastic oscillations can explain relapses/remissions in clinical manifestations of thedisease as caused by endogenous stochasticity of the immune processes.Blyuss and Nicholson [41, 42] have proposed and analysed a mathematical model of immuneresponse to a viral infection that explicitly takes into account the populations of two types of T cellswith different activation thresholds and also allows for infection and autoimmune response to occurin 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 withepisodes of high viral production followed by long periods of quiescence. Such behaviour, that in theclinical observation could be associated with relapses and remissions, has been observed in a numberof autoimmune diseases, such as MS, autoimmune thyroid disease and uveitis [43, 44, 45]. Despite itssuccesses, this model has a limitation that the periodic oscillations are only possible when the amountof free virus and the number of infected cells are also exhibiting oscillations, while in laboratory andclinical 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 etal. (2017) have recently extended the model of Blyuss and Nicholson to also include the populationof 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 modelwith 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 themethodology of continuous-time Markov chains (CTMC) to derive a Kolmogorov, or chemical masterequation, describing the dynamics of a probability distribution of finding the system in a particularstate. To make further analytical and numerical progress, we derive an Itˆo stochastic differentialequation, whose solutions provide similar stochastic paths to those of the CTMC models. This thenallows us to numerically study the stationary multivariate probability distributions for the states inthe model, explore stochastic amplification, determine how the magnitude of stochastic fluctuationsaround deterministic steady states depends on various parameters, and investigate the effects of initialconditions on the outcome in the case of bi-stability between different dynamical states. These resultssuggest that the experimentally observed variation in the progression of autoimmune disease canbe attributed to stochastic amplification, and they also provide insights into how the variance offluctuations depends on parameters, which can guide new laboratory experiments. In a recent paper we introduced and analysed a deterministic model for autoimmune dynamics withaccount for the populations of T cells with different activation thresholds and cytokines (Fatehi etal. 2017). The analysis showed that depending on parameters and initial conditions, the model cansupport the regimes of normal disease clearance , where an initial infection is cleared without furtherconsequences for immune dynamics, chronic infection characterised by a persistent presence of infectedcells in the body, and the state of autoimmune behaviour where after clearance of initial infection, theimmune system supports stable endogenous oscillations in the number of autoreactive T cells, whichcan be interpreted in the clinical practice of autoimmune disease as periods of relapses and remissions.This work extended earlier results on modelling the effects of tunable activation thresholds [41, 42] byincluding regulatory T cells, as well as the cytokine mediating proliferation and activity of differenttypes of T cells.A deterministic model for immune response to a viral infection, as illustrated in a diagram shownin Figure 1, has the form 3igure 1: A schematic diagram of immune response to infection. Blue indicates host cells (susceptibleand infected), red denotes T cells (na¨ıve, regulatory, normal activated, and autoreactive), yellow showscytokines (interleukin-2). 4
Sdt = rS (cid:18) − SN (cid:19) − βSF − µ a T aut S,dFdt = βSF − d F F − µ F T nor F − µ a T aut F,dT in dt = λ in − d in T in − αT in F,dT reg dt = λ r − d r T reg + p αT in F + ρ IT reg ,dT nor dt = p αT in F − d n T nor + ρ IT nor ,dT aut dt = (1 − p − p ) αT in F − d a T aut − δT reg T aut + ρ IT aut ,dIdt = σ T nor + σ T aut − d i I, (1)where S ( t ) is the number of susceptible organ cells, F ( t ) is the number of infected cells, T in ( t ) isthe number of na¨ıve T cells, T reg ( t ) is the number of regulatory T cells, T nor ( t ) is the number ofactivated T cells which recognise foreign antigen and destroy infected cells, T aut ( t ) is the numberof autoreactive T cells which destroy cells presenting both foreign and self-antigen, and I ( t ) is theamount 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 , andthey can become infected at rate β by already infected cells that are producing new virus particles.Unlike earlier models [41, 42] and (Fatehi et al. 2017), we consider the situation where the processof producing virions by infected cells is quite fast, hence, we do not explicitly incorporate a separatecompartment for free virus. Regarding immune response, we assume that na¨ıve T cells remain inhomeostasis, and upon activation at rate α by a signal from infected cells, a proportion p of themwill develop into regulatory T cells, a proportion p will become normal activated T cells able todestroy infected cells at rate µ F , and the remaining T cells will become autoreactive, in which casetheir threshold for activation by susceptible cells is reduced, and hence, they will be destroying bothinfected and susceptible host cells at rate µ a . The effect of regulatory T cells is in reducing the numberof autoreactive T cells at rate δ , and regulatory T cells are themselves assumed to be in a state ofhomeostasis. Finally, normal and autoreactive T cells produce IL-2 at rates σ and σ , and IL-2 in turnfacilitates 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 inthe absence of B cells [46] and the fact that the development of antibodies can itself depend on priorT cell activation with bacteria [47], 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 CTMCmodel based on the ODE model (1) using the methodology developed earlier in the context of modellingstochastic effects in epidemic and immunological models [48, 49, 36]. To this end, we introducevariables X ( t ) , . . . , X ( t ) ∈ { , , , . . . } as discrete random variables representing the number ofuninfected cells, infected cells, na¨ıve T cells, regulatory T cells, normal activated T cells, autoreactiveT cells, and interleukin-2 at time t , respectively. Let the initial condition be fixed as X = ( X (0) , . . . , X (0)) = ( n , n , n , n , n , n , n ) . n = ( n , n , n , n , n , n , n ) with n i ∈ { , , , ... } at time t can be defined as P ( n , t ) = Prob { X ( t ) = n | X (0) = X } . Let ∆ t be sufficiently small such that ∆ X i ( t ) = X i ( t + ∆ t ) − X i ( t ) ∈ {− , , } for 1 ≤ i ≤ X = i | X = n ) = q ∆ t + o (∆ t ) , i = (1 , , , , , , ,q ∆ t + o (∆ t ) , i = ( − , , , , , , ,q ∆ t + o (∆ t ) , i = ( − , , , , , , ,q ∆ t + o (∆ t ) , i = (0 , , , , , , ,q ∆ t + o (∆ t ) , i = (0 , , − , , , , ,q ∆ t + o (∆ t ) , i = (0 , , − , , , , ,q ∆ t + o (∆ t ) , i = (0 , , − , , , , ,q ∆ t + o (∆ t ) , i = (0 , , − , , , , ,q ∆ t + o (∆ t ) , i = (0 , − , , , , , ,q ∆ t + o (∆ t ) , i = (0 , , , , , , ,q ∆ t + o (∆ t ) , i = (0 , , , − , , , ,q ∆ t + o (∆ t ) , i = (0 , , , , , , ,q ∆ t + o (∆ t ) , i = (0 , , , , − , , ,q ∆ t + o (∆ t ) , i = (0 , , , , , , ,q ∆ t + o (∆ t ) , i = (0 , , , , , − , ,q ∆ t + o (∆ t ) , i = (0 , , , , , , ,q ∆ t + o (∆ t ) , i = (0 , , , , , , − , − (cid:80) i =1 q i ∆ t + o (∆ t ) , i = (0 , , , , , , ,o (∆ t ) , otherwise, (2)where q = b n + b n , q = d n + d n + µ a n n , q = βn n , q = λ in ,q = d in n , q = p αn n , q = p αn n , q = (1 − p − p ) αn n ,q = ( d F + µ F n + µ a n ) n , q = λ r + ρ n n , q = d r n , q = ρ n n ,q = d n n , q = ρ n n , q = ( d a + δn ) n , q = σ n + σ n , q = d i n . Here, b n + b n and d n + d n are natural birth and death rates for uninfected cells with b − d = r and d − b = r/N [48].The probabilities P ( n , t ) satisfy the following master equation (forward Kolmogorov equation) dP ( n , t ) dt = { ( ε − − q + ( ε +1 − q + ( ε +1 ε − − q + ( ε − − q + ( ε +3 − q + ( ε +3 ε − − q + ( ε +3 ε − − q + ( ε +3 ε − − q + ( ε +2 − q + ( ε − − q + ( ε +4 − q + ( ε − − q + ( ε +5 − q + ( ε − − q + ( ε +6 − q + ( ε − − q + ( ε +7 − q } P ( n , t ) . (3)6able 1: Possible state changes ∆ Y during a small time interval ∆ t i (∆ Y ) Ti Probability P i ∆ t , , , , , ,
0) ( b Y + b Y )∆ t − , , , , , ,
0) ( d Y + d Y + µ a Y Y )∆ t − , , , , , , βY Y ∆ t , , , , , , λ in ∆ t , , − , , , , d in Y ∆ t , , − , , , , p αY Y ∆ t , , − , , , , p αY Y ∆ t , , − , , , ,
0) (1 − p − p ) αY Y ∆ t , − , , , , ,
0) ( d F + µ F Y + µ a Y ) Y ∆ t
10 (0 , , , , , ,
0) ( λ r + ρ Y Y )∆ t
11 (0 , , , − , , , d r Y ∆ t
12 (0 , , , , , , ρ Y Y ∆ t
13 (0 , , , , − , , d n Y ∆ t
14 (0 , , , , , , ρ Y Y ∆ t
15 (0 , , , , , − ,
0) ( d a + δY ) Y ∆ t
16 (0 , , , , , ,
1) ( σ Y + σ Y )∆ t
17 (0 , , , , , , − d i Y ∆ t
18 (0 , , , , , ,
0) 1 − (cid:80) i =1 P i ∆ t where the operators ε ± i are defined as follows, ε ± i f ( n , n , n , n , n , n , n , t ) = f ( n , ..., n i ± , ..., n , t ) , for each 1 ≤ i ≤ , and if n i < ≤ 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 chal-lenging task. Normally, the number of events occurring in a small time step in the CTMC model isextremely large, hence using the CTMC model for plotting stochastic trajectories is very computa-tionally intensive [50]. A much more computationally efficient approach is to use chemical Langevinequations [51, 52], also known as Itˆo stochastic differential equation (SDE) models, which provide verysimilar sample paths to those of the CTMC models [50]. While both Itˆo and Stratonovich interpre-tations of stochastic calculus can be applied [53], in biological applications Itˆo formulation is morefrequently used due to its non-anticipatory nature and a closer connection to numerical implementa-tion [54, 55, 48]. To derive Itˆo SDE model, let Y ( t ) = ( Y ( t ) , Y ( t ) , Y ( t ) , Y ( t ) , Y ( t ) , Y ( t ) , Y ( t )) be a continuousrandom vector for the sizes of various cell compartments at time t . Similar to the CTMC model, weassume that ∆ t is small enough so that during this time interval at most one change can occur in statevariables. These changes together with their probabilities are listed in Table 1, which is again basedon Figure 1 and transitions in the CTMC model (2). Using this table of possible state changes, onecan compute the expectation vector and covariance matrix of ∆ Y for sufficiently small ∆ t [50, 56].7he expectation vector to order ∆ t is given by E (∆ Y ) ≈ (cid:88) i =1 P i (∆ Y ) i ∆ t = µ ∆ t, where µ = P − P − P P − P P − P − P − P − P P + P − P P + P − P P + P − P P − P is the drift vector, which can be easily seen to be identical to the right-hand side of the deterministicmodel (1). The covariance matrix is obtained by keeping terms of order ∆ t only, i.e.cov(∆ Y ) = E (cid:2) (∆ Y )(∆ Y ) T (cid:3) − E [∆ Y ] ( E [∆ Y ]) T ≈ E (cid:2) (∆ Y )(∆ Y ) T (cid:3) = (cid:80) i =1 P i (∆ Y ) i (∆ Y i ) T ∆ t = Σ∆ t, where Σ = P + P + P − P − P P + P P + P + P + P + P − P − P − P
00 0 − P P + P + P − P P + P + P − P P + P + P
00 0 0 0 0 0 P + P is a 7 × H definedaccording to HH T = Σ. Although this matrix is not unique, different forms of this matrix giveequivalent systems [55, 56].If one rewrites the covariance matrix Σ in the formΣ = U W
00 0 Z , with U = (cid:18) P + P + P − P − P P + P (cid:19) , Z = P + P , and W = P + P + P + P + P − P − P − P − P P + P + P − P P + P + P − P P + P + P , we can define three matrices H , H and H as follows, H = (cid:18) √ P + P −√ P √ P √ P (cid:19) , H = (cid:112) P + P , = √ P + P −√ P −√ P −√ P √ P √ P + P √ P √ P + P
00 0 0 √ P √ P + P . Now if we consider H = H H
00 0 H , then HH T = Σ, where H is a 7 ×
11 matrix. The Itˆo SDE model now has the form (cid:40) d Y ( t ) = µ dt + Hd W ( t ) , Y (0) = ( A (0) , F (0) , T in (0) , T reg (0) , T nor (0) , T aut (0) , I (0)) T , (4)and W ( t ) = [ W ( t ) , W ( t ) , ..., W ( t )] T is a vector of eleven independent Wiener processes [55].In order to make further analytical progress, we find an approximate probability density functionfor the model (4) as given by an approximate solution of the master equation [55, 57]. Let P ( Y , t ) bethe probability density function of the model (4). Then P ( Y , t ) satisfies the following Fokker-Planckequation [55, 58] which is an approximation of the master equation ∂P ( Y , t ) ∂t = − (cid:88) i =1 ∂∂y i [ µ i P ( Y , t )] + 12 (cid:88) i =1 7 (cid:88) j =1 ∂ ∂y i ∂y j [Σ ij P ( Y , t )] ,P ( Y ,
0) = δ ( Y − Y ) . By solving this PDE, one can find the probability density function of our model, but since thisequation is high-dimensional and nonlinear, solving it analytically is impossible. Hence, we use anotherapproach, a so-called system size expansion or van Kampen’s Ω-expansion [57], which is a method forconstructing a continuous approximation to a discrete stochastic model [36, 37], which allows one tostudy stochastic fluctuations around deterministic attractors [59].
In order to apply the van Kampen’s approach, we consider fluctuations within a systematic expansionof the master equation for a large system size Ω. Specifically, we write each n i ( t ) as a deterministicpart of order Ω plus a fluctuation of order Ω / as follows, n i ( t ) = Ω x i ( t ) + Ω / ζ i ( t ) , i = 1 , . . . , , (5)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 (cid:0) Ω x + Ω / ζ , t (cid:1) , which implies dP ( n , t ) dt = ∂ Π ∂t − (cid:88) i =1 Ω / dx i dt ∂ Π ∂ζ i . (6)To expand the master equation (3) in a power series in Ω − / , we use the following expansion for thestep operators ε ± i = 1 ± Ω − / ∂∂ζ i + 12 Ω − ∂ ∂ζ i ± · · · . (7)9ubstituting expressions (6) and (7) into the master equation (see Supplementary Material fordetails) and collecting terms of order Ω / yields the following deterministic model for macroscopicbehaviour dx dt = b x + (cid:101) b x − d x − (cid:101) d x − (cid:101) βx x − (cid:101) µ a x x ,dx dt = (cid:101) βx x − d F x − (cid:101) µ F x x − (cid:101) µ a x x ,dx dt = (cid:101) λ in − d in x − (cid:101) αx x ,dx dt = (cid:101) λ r − d r x + p (cid:101) αx x + (cid:101) ρ x x ,dx dt = p (cid:101) αx x − d n x + (cid:101) ρ x x ,dx dt = (1 − p − p ) (cid:101) αx x − d a x − (cid:101) δx x + (cid:101) ρ x x ,dx dt = σ x + σ x − d i x , (8)where b = (cid:101) b Ω , d = (cid:101) d Ω , β = (cid:101) β Ω , µ a = (cid:101) µ a Ω , µ F = (cid:101) µ F Ω , α = (cid:101) α Ω , δ = (cid:101) δ Ω ,ρ i = (cid:101) ρ i Ω , i = 1 , , , λ in = (cid:101) λ in Ω , λ r = (cid:101) λ r Ω . Model (8) has been analysed in (Fatehi et al. 2017), and it can have at most four biologicallyfeasible steady states. The first one, a disease-free steady state, is given by S ∗ = (cid:32) b − d (cid:101) d − (cid:101) b , , (cid:101) λ in d in , (cid:101) λ r d r , , , (cid:33) , and it is stable if d F > (cid:101) β . The second and third steady states can be found as S ∗ = , , (cid:101) λ in d in , x ∗ , , d i (cid:16) d a + (cid:101) δx ∗ (cid:17)(cid:101) ρ σ , d a + (cid:101) δx ∗ (cid:101) ρ , and S ∗ = (cid:101) ρ σ ( b − d ) − (cid:101) µ a d i (cid:16) d a + (cid:101) δx ∗ (cid:17)(cid:101) ρ σ (cid:16) (cid:101) d − (cid:101) b (cid:17) , , (cid:101) λ in d in , x ∗ , , d i (cid:16) d a + (cid:101) δx ∗ (cid:17)(cid:101) ρ σ , d a + (cid:101) δx ∗ (cid:101) ρ , where x ∗ satisfies the following quadratic equation (cid:101) ρ (cid:101) δ ( x ∗ ) + ( (cid:101) ρ d a − (cid:101) ρ d r ) x ∗ + (cid:101) ρ (cid:101) λ r = 0 . (9)These steady states are stable, provided σ (cid:101) µ a d i K < d a + (cid:101) δx ∗ (cid:101) ρ < d n (cid:101) ρ , (cid:101) δ (cid:101) ρ ( x ∗ ) > (cid:101) λ r (cid:101) ρ , (cid:101) ρ (cid:101) λ r + (cid:101) ρ d i (cid:101) λ r x ∗ − (cid:101) ρ d i d a ( x ∗ ) − (cid:101) δ ( (cid:101) ρ d a + (cid:101) ρ d i )( x ∗ ) − (cid:101) ρ (cid:101) δ ( x ∗ ) > , K = 1 for S ∗ , and K = (cid:16) (cid:101) β − d F (cid:17) / (cid:16) (cid:101) β (cid:17) for S ∗ . Biologically, the steady state S ∗ representsthe death of organ cells, while S ∗ corresponds to an autoimmune regime.The last steady state S ∗ has all of its components positive and corresponds to the state of chronicinfection.At the next order, stochastic fluctuations are determined by linear stochastic processes, hence, thisis known as a linear noise approximation [57, 60]. The dynamics of these fluctuations is described bythe following linear Fokker-Planck equation ∂ Π( ζ , t ) ∂t = − (cid:88) i,j A ij ∂∂ζ i ( ζ j Π) + 12 (cid:88) i,j B ij ∂ Π ∂ζ i ∂ζ j , (10)where A is the Jacobian matrix of system (8) A = b + 2 (cid:101) b x − d − (cid:101) d x − (cid:101) µ a x − (cid:101) βx − (cid:101) βx − (cid:101) µ a x (cid:101) βx (cid:101) βx − d F − (cid:101) µ F x − (cid:101) µ a x − (cid:101) µ F x − (cid:101) µ a x − (cid:101) αx − d in − (cid:101) αx p (cid:101) αx p (cid:101) αx (cid:101) ρ x − d r (cid:101) ρ x p (cid:101) αx p (cid:101) αx (cid:101) ρ x − d n (cid:101) ρ x − p − p ) (cid:101) αx (1 − p − p ) (cid:101) αx − (cid:101) δx (cid:101) ρ x − d a − (cid:101) δx (cid:101) ρ x σ σ − d i , and B is a 7 × B ij = b x + (cid:101) b x + d x + (cid:101) d x + (cid:101) βx x + (cid:101) µ a x x , if ( i, j ) = (1 , , (cid:101) βx x + d F x + (cid:101) µ F x x + (cid:101) µ a x x , if ( i, j ) = (2 , , (cid:101) λ in + d in x + (cid:101) αx x , if ( i, j ) = (3 , , (cid:101) λ r + d r x + p (cid:101) αx x + (cid:101) ρ x x , if ( i, j ) = (4 , ,p (cid:101) αx x + d n x + (cid:101) ρ x x , if ( i, j ) = (5 , , (1 − p − p ) (cid:101) αx x + d a x + (cid:101) δx x + (cid:101) ρ x x , if ( i, j ) = (6 , ,σ x + σ x + d i x , if ( i, j ) = (7 , , − (cid:101) βx x , if ( i, j ) = (1 ,
2) or (2 , , − p (cid:101) αx x , if ( i, j ) = (3 ,
4) or (4 , , − p (cid:101) αx x , if ( i, j ) = (3 ,
5) or (5 , , − (1 − p − p ) (cid:101) αx x , if ( i, j ) = (3 ,
6) or (6 , , , 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 characterise it [61, 62]. Due to the way the system sizeexpansion was introduced in (5), the mean values of fluctuations for all variables are zero, i.e. (cid:104) ζ i ( t ) (cid:105) =0 for all 1 ≤ i ≤
7, while the covariance matrix Ξ with Ξ ij = (cid:104) ζ i ( t ) ζ j ( t ) (cid:105) − (cid:104) ζ i ( t ) (cid:105)(cid:104) ζ j ( t ) (cid:105) = (cid:104) ζ i ( t ) ζ j ( t ) (cid:105) satisfies the following equation [57, 62] ∂ t Ξ = A Ξ + Ξ A T + B, (11)where A T is the transpose of A .We are mainly interested in the dynamics of fluctuations when the oscillations of the deterministicmodel have died out, and the system is in a stationary state, i.e. the fluctuations take place aroundthe steady states [59]. 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 state11able 2: Table of parametersParameter value Parameter value b d r (cid:101) b p d (cid:101) ρ / (cid:101) d p (cid:101) β d n (cid:101) µ a / (cid:101) ρ / d F d a (cid:101) µ F / (cid:101) δ / (cid:101) λ in (cid:101) ρ / d in σ (cid:101) α σ (cid:101) λ r d i A Ξ + Ξ A T + B = 0 . In order to be able to relate the results of this analysis to simulations, it is convenient to express thecovariance matrix in terms of actual numbers of cells in each compartment, rather than deviations fromstationary values. To this end, we instead use the covariance matrix C defined as C ij = (cid:104) ( n i −(cid:104) n i (cid:105) )( n j −(cid:104) n j (cid:105) ) (cid:105) , which, in light of the relation C ij = ΩΞ ij , satisfies the following Lyapunov equation [62] AC + CA T + Ω B = 0 . (12)This equation can be solved numerically for each of the stable steady states to determine the varianceof fluctuations around that steady state depending on system parameters. To simulate the dynamics of the model, we solve the system (4) numerically using the Euler-Maruyamamethod with parameter values given in Table 2, and Ω = 1000. The initial condition is chosen to beof the form ( x (0) , x (0) , x (0) , x (0) , x (0) , x (0) , x (0)) = (18 , , . , . , , , , (13)which corresponds to a small number of host cells being initially infected.Figure 2 shows the results of 20000 simulations with the initial condition (13) and σ = 1. In thedeterministic model (8), for σ = 1 both steady states S ∗ (disease-free) and S ∗ (autoimmune state) arestable, but with the initial condition (13) the system is in the basin of attraction of S ∗ . In the stochasticmodel, the majority of trajectories also enter the attraction region of S ∗ , but a small proportion ofthem went into the basin of attraction of S ∗ . This figure illustrates a single stochastic path around S ∗ ,and a single stochastic path around S ∗ , together with the deterministic trajectory. These individualsolutions indicate that whilst deterministically, the system exhibits decaying oscillations around S ∗ ,the same behaviour is observed in the stochastic simulations only upon taking an average of a very12igure 2: Numerical simulation of the model (4) with parameter values from Table 2, σ = 1, andthe initial condition (13). Red curves are two sample paths that have entered the basins of attractionof S ∗ or S ∗ , black curve is the deterministic trajectory from (1), and the shaded areas indicate theregions of one standard deviation from the mean.large number of simulations. At the same time, individual realisations exhibit sustained stochasticoscillations in a manner similar to that observed in models of stochastic amplification in epidemics [39,40]. Figure 2 also illustrates the size of areas of one standard deviation from the mean for trajectoriesin the basins of attraction S ∗ and S ∗ , in which individual stochastic trajectories may exhibit stochasticoscillations [63, 38].Figures 3A and B show temporal evolution of the probability distribution in the case of bi-stabilitybetween the steady states S ∗ and S ∗ , as illustrated in Figure 2. They indicate that after some initialtransient, the system reaches a stationary bimodal normal distribution. The width of the probabilitydistribution around each stable steady state, as described by its variance or standard deviation, givesthe size of fluctuations around this steady state observed in individual stochastic realisations, as isshown in Fig. 2. Similar behaviour has been observed in stochastic realisations of other deterministicmodels with bi-stability [64, 65, 66]. For the parameter values given in Table 2, the deterministicsystem exhibits a bi-stability between S ∗ and S ∗ , and with the initial condition( x (0) , x (0) , x (0) , x (0) , x (0) , x (0) , x (0)) = (18 , , . , . , , , , (14)it is in the basin of attraction of S ∗ . Due to stochasticity, the stationary probability distribution inthis case is also bimodal, with the majority of solutions being distributed around S ∗ , and a very smallnumber being centred around S ∗ , as can be seen in Figures 3C and D. Increasing the system size Ω isknown to result in the bimodal distribution becoming unimodal due to the size of fluctuations scalingas Ω − / , which results in a reduced variability in trajectories [67, 66], and the same conclusion holdsfor the system (4).To gain better insights into the role of initial conditions, in Figure 4 we fix all parameter values,13igure 3: Probability distribution of solutions out of 20000 simulations. (A) and (B) with parametersfrom Table 2, but σ = 1 and the initial condition (13). (C) and (D) with parameters from Table 2and the initial condition (14). In (A) and (C), the probability histogram is fit to a bimodal normaldistribution at different times. (B) and (D) illustrate stationary joint probability histograms.and vary initial numbers of infected cells and regulatory T cells. For the parameter combinationillustrated in Figure 4A, the deterministic model exhibits a bi-stability between a stable disease-freesteady state S ∗ and a periodic oscillation around the state S ∗ , which biologically corresponds to anautoimmune regime. In the deterministic case, the black boundary provides a clear separation of thebasins of attraction of these two dynamical states, in a manner similar to that investigated recentlyin the context of within-cell dynamics of RNA interference [68]. For stochastic simulations, the colourindicates the probability of the solution going to a disease-free state S ∗ , and it shows that even inthe case where deterministically the system is in the basin of attraction of one of the states, thereis a non-zero probability that it will actually end up at another state, with this probability varyingsmoothly across the deterministic basin boundary. This figure suggests that if the initial number ofinfected cells is sufficiently small, or if the number of regulatory T cells is sufficiently large, the systemtends to clear the infection and approach the disease-free state. On the contrary, for higher numbersof infected cells and lower numbers of regulatory cells, autoimmune regime appears to be a more likelyoutcome. Qualitatively similar behaviour is observed for another combination of parameters illustrated14igure 4: Probability of solution entering and staying in the basin of attraction of the disease-freesteady state S ∗ in the bi-stability regime with A (0) = 18000 and T in (0) = 7200. Black curves arethe boundaries between different basins of attraction in the deterministic model. (A) With parametervalues from Table 2, (cid:101) λ r = 45 and (cid:101) µ a = 10 /
9, in the region below the black curve, the deterministicmodel exhibits a periodic solution around S ∗ , and above this curve is the deterministic basin ofattraction of S ∗ . (B) With parameter values from Table 2, area below the black curve is the basin ofattraction of S ∗ , and above it is again the basin of attraction of S ∗ .in Figure 4B, in which case the deterministic system has a bi-stability between a disease-free steadystate S ∗ , and a state S ∗ which represents the death of host cells.In order to understand how biological parameters affect the size of fluctuations around steadystates, in Figure 5 we explore several parameter planes by first identifying parameter regions wherethe deterministic system has a stable steady state S ∗ , and then for each combination of parametersinside these regions, we use the Bartels-Stewart method [69, 70] to numerically solve the Lyapunovequation (12) and compute the variance in the number of regulatory T cell when the deterministicmodel is at the steady state S ∗ . The value of variance gives the square of the magnitude of oscillationsobserved in individual stochastic realisations. One should note that getting closer to the deterministicboundary of stability of S ∗ 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 stableis 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 whichregulatory T cells suppress autoreactive T cells; it decreases with the higher rate of production ofregulatory T cells, and it appears to not depend on the rate at which autoreactive T cells destroyinfected cells, or on the infection rate. In this paper we have analysed stochastic aspects of immune response against a viral infection withaccount for the populations of T cells with different activation thresholds, as well as cytokines medi-ating T cell activity. The CTMC model has provided an exact master equation, for which we applieda van Kampen’s expansions to derive a linear Fokker-Planck equation that characterises fluctuationsaround the deterministic solutions. We have also explored actual stochastic trajectories of the system15igure 5: Variance of the number of regulatory T cells T reg with parameter values from Table 2.Coloured regions indicate areas in respective parameter planes in which the autoimmune steady state S ∗ is deterministically stable.by deriving an SDE model and solving it numerically.One biologically important aspect we have looked at is the influence of stochasticity on the dynam-ics of the system in the case where deterministically it exhibits a bi-stability between either two steadystates, or a steady state and a periodic solution. In such a situation, bi-stability in the deterministicversion of the model translates in the stochastic case into a stationary bimodal distribution for theprobability density. To obtain further insights into details of how stochasticity affects bi-stability, wehave investigated how for the fixed parameter values time evolution of the system changes dependingon the initial numbers of the regulatory T cells and infected cells.Our analysis reinforces the need to distinguish mean dynamics from individuals realisations: wherein the deterministic case the system can approach a stable steady state (which represents mean be-haviour of a very large number of simulations), individual realisations can exhibit sustained stochasticoscillations around that steady state, as we have seen in numerical simulations. Since in the clini-cal or laboratory setting one is usually dealing with single measurements of some specific biologicalquantities rather than their averaged values, the stochastic oscillations exhibited by our model mayquite well explain observed variability in the measured levels of infection or T cell populations. Tobetter understand the magnitude of stochastic fluctuations around the deterministic steady states, we16ave solved the Lyapunov equation, which has provided us with a quantitative information on thedependence of variance of fluctuations on system parameters.There are several directions in which the work presented in this paper can be extended. In termsof 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 [71], or the memory T cells [72, 73]. Whilstwe have used numerical simulations to compute the probability of attraction to a given steady statein the case of bi-stability, one could approach the same problem theoretically from the perspectiveof computing extinction probability within the framework of the CTMC model [50, 74]. The vanKampen’s system size expansion could yield an expression for the power spectrum, which allowsone to compute the peak frequency and amplification [59, 75, 76, 39]. From a practical perspective,future work could focus on validating theoretical results presented in this paper using experimentalmeasurements of the progress of autoimmune disease in animal hosts, with experimental autoimmuneuveoretinitis (EAU), an autoimmune inflammation in the eyes, being one interesting possibility. Inone such recent experiment, all animals were genetically identical C57BL/6 mice, but once the EAUwas induced in them through inoculation, the autoimmune disease then progressed at slightly differentrates [77], (Boldison and Nicholson, 2012) and the measured variability in the numbers of infected cellsand T cell responses could be compared to theoretical estimates of the variance as predicted by ourmodel. From a clinical perspective, comparison of variance in the measured populations of differentcells with the model conclusions will facilitate an efficient parameter identification and provide a setof prognostic criteria for the progress of autoimmunity, which can be used for risk stratification andassessment of patients with autoimmune disease. Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financialrelationships that could be construed as a potential conflict of interest.
Author Contributions
YK and KB designed the model; FF, SK and AR performed the analysis and simulations, and producedthe figures. All authors drafted and edited the manuscript.
Funding
FF has been sponsored by the Chancellor’s Award from Sussex University. AR acknowledges fundingfor her PhD studies from the Engineering and Physical Sciences Research Council (EPSRC) throughiCASE Award EP/N509784/1.
Supplementary Material
The Supplementary Material for this article is attached.
References [1] E. C. Kerr, D. A. Copland, A. D. Dick, and L. B. Nicholson, “The dynamics of leukocyte infil-tration in experimental autoimmune uveoretinitis,”
Eye Res. , vol. 27, pp. 527–535., 2008.172] E. Prat and R. Martin, “The immunopathogenesis of multiple sclerosis,”
J. Rehabil. Res. Dev. ,vol. 39, pp. 187–199., 2002.[3] P. Santamaria, “The long and winding road to understanding and conquering type 1 diabetes,”
Immunity , vol. 32, pp. 437–445., 2010.[4] R. S. Fujinami, “Can virus infections trigger autoimmune disease?,”
J. Autoimmun. , vol. 16,pp. 229–234., 2001.[5] M. G. von Herrath and M. B. A. Oldstone, “Virus-induced autoimmune disease,”
Curr. Opin.Immunol. , vol. 8, pp. 878–885., 1996.[6] A. M. Ercolini and S. D. Miller, “The role of infections in autoimmune disease,”
Clin. Exp.Immunol. , vol. 155, pp. 1–15., 2009.[7] L. A. Segel, E. J¨ager, D. Elias, and I. R. Cohen, “A quantitative model of autoimmune diseaseand T-cell vaccination: does more mean less?,”
Immunol. Today , vol. 16, no. 2, pp. 80–84., 1995.[8] J. A. M. Borghans and R. J. De Boer, “A minimal model for T-cell vaccination,”
Proc. R. Soc.B , vol. 259, no. 1355, pp. 173–178., 1995.[9] J. A. M. Borghans, R. J. De Boer, E. Sercarz, and V. Kumar, “T cell vaccination in experimentalautoimmune encephalomyelitis: a mathematical model,”
J. Immunol. , vol. 161, no. 3, pp. 1087–1093., 1998.[10] K. Le´on, R. Perez, A. Lage, and J. Carneiro, “Modelling T-cell-mediated suppression dependenton interactions in multicellular conjugates,”
J. Theor. Biol. , vol. 207, no. 2, pp. 231–254., 2000.[11] K. Le´on, A. Lage, and J. Carneiro, “Tolerance and immunity in a mathematical model of T-cellmediated suppression,”
J. Theor. Biol. , vol. 225, no. 1, pp. 107–126., 2003.[12] K. Le´on, J. Faro, A. Lage, and J. Carneiro, “Inverse correlation between the incidences of autoim-mune disease and infection predicted by a model of T cell mediated tolerance,”
J. Autoimmun. ,vol. 22, no. 1, pp. 31–42., 2004.[13] J. Carneiro, T. Paix˜ao, and D. Milutinovic, “Immunological self-tolerance: lessons from mathe-matical modeling,”
J. Comput. Appl. Math. , vol. 184, pp. 77–100., 2005.[14] S. Iwami, Y. Takeuchi, Y. Miura, T. Sasaki, and T. Kajiwara, “Dynamical properties of autoim-mune disease models: tolerance, flare-up, dormancy,”
J. Theor. Biol. , vol. 246, no. 4, pp. 646–659.,2007.[15] S. Iwami, Y. Takeuchi, K. Iwamoto, Y. Naruo, and M. Yasukawa, “A mathematical design ofvector vaccine against autoimmune disease,”
J. Theor. Biol. , vol. 256, no. 3, pp. 382–392., 2009.[16] H. K. Alexander and L. M. Wahl, “Self-tolerance and autoimmunity in a regulatory T cell model,”
Bull. Math. Biol. , vol. 73, no. 1, pp. 33–71., 2011.[17] N. J. Burroughs, B. M. P. M. Oliveira, and A. A. Pinto, “Autoimmunity arising from bystanderproliferation of T cells in an immune response model,”
Math. Comput. Mod. , vol. 53, pp. 1389–1393., 2011.[18] N. J. Burroughs, M. Ferreira, B. M. P. M. Oliveira, and A. A. Pinto, “A transcritical bifurcationin an immune response model,”
J. Diff. Eqns. Appl. , vol. 17, pp. 1101–1106., 2011.1819] R. Root-Bernstein, “Theories and modeling of autoimmunity,”
J. Theor. Biol. , vol. 375, pp. 1–124., 2015.[20] Z. Grossman and W. E. Paul, “Adaptive cellular interactions in the immune system: the tun-able activation threshold and the significance of subthreshold responses,”
Proc. Natl. Acad. Sci.U.S.A. , vol. 89, pp. 10365–10369., 1992.[21] Z. Grossman and W. E. Paul, “Self-tolerance: context dependent tuning of T cell antigen recog-nition,”
Sem. Immunol. , vol. 12, pp. 197–203., 2000.[22] Z. Grossman and A. Singer, “Tuning of activation thresholds explains flexibility in the selectionand development of T cells in the thymus,”
Proc. Natl. Acad. Sci. U.S.A. , vol. 93, pp. 14747–14752., 1996.[23] G. Altan-Bonnet and R. N. Germain, “Modeling T cell antigen discrimination based on feedbackcontrol of digital ERK responses,”
PLoS Biol. , vol. 3, p. e356., 2005.[24] H. A. van den Berg and D. A. Rand, “Dynamics of T cell activation threshold tuning,”
J. Theor.Biol. , vol. 228, pp. 397–416., 2004.[25] A. Scherer, A. Noest, and R. J. deBoer, “Activation-threshold tuning in an affinity model for theT-cell repertoire,”
Proc. R. Soc. B , vol. 271, pp. 609–616., 2004.[26] A. J. Noest, “Designing lymphocyte functional structure for optimal signal detection: voil´a,”
J.Theor. Biol. , vol. 207, pp. 195–216., 2000.[27] A. D. Bitmansour, D. C. Douek, V. C. Maino, and L. J. Picker, “Direct ex vivo analysis of humanCD4 + memory T cell activation requirements at the single clonotype level,” J. Immunol. , vol. 169,pp. 1207–1218., 2002.[28] L. B. Nicholson, A. C. Anderson, and V. K. Kuchroo, “Tuning T cell activation threshold andeffector function with cross-reactive peptide ligands,”
Int. Immunol. , vol. 12, pp. 205–213., 2000.[29] P. S. R¨omer, S. Berr, E. Avota, S.-Y. Na, M. Battaglia, I. ten Berge, and et al., “Preculture ofPBMC at high cell density increases sensitivity of T-cell responses, revealing cytokine release byCD28 superagonist TGN1412,”
Blood , vol. 118, pp. 6772–6782., 2011.[30] I. Stefanova, J. R. Dorfman, and R. N. Germain, “Self-recognition promotes the foreign antigensensitivity of naive T lymphocytes,”
Nature , vol. 420, pp. 429–434., 2002.[31] A. S. Perelson and G. Weisbuch, “Immunology for physicists,”
Rev. Mod. Phys. , vol. 69, no. 4,pp. 1219–1267., 1997.[32] E. K. Deenick, A. V. Gett, and P. D. Hodgkin, “Stochastic model of T cell proliferation: A calculusrevealing IL-2 regulation of precursor frequencies, cell cycle time, and survival,”
J. Immunol. ,vol. 170, pp. 4963–4972., 2003.[33] J. N. Blattman, K. Sourdive D J D, Murali-Krishna, R. Ahmed, and J. D. Altman, “Evolution ofthe T cell repertoire during primary, memory, and recall response to viral infection,”
J. Immunol. ,vol. 165, pp. 6081–6090., 2000.[34] V. Detours and A. S. Perelson, “The paradox of alloreactivity and self MHC restriction: Quan-titative analysis and statistics,”
Proc. Natl. Acad. Sci. U.S.A. , vol. 97, pp. 8479–8483., 2000.1935] D. L. Chao, M. P. Davenport, S. Forrest, and A. S. Perelson, “A stochastic model of cytotoxic Tcell responses,”
J. Theor. Biol. , vol. 228, pp. 227–240, 2004.[36] E. R. Stirk, G. Lythe, H. A. van den Berg, G. A. D. Hurst, and C. Molina-Par´ıs, “The limitingconditional probability distribution in a stochastic model of T cell repertoire maintenance,”
Math.Biosci. , vol. 224, pp. 74–86., 2010.[37] E. R. Stirk, G. Lythe, H. A. van den Berg, and C. Molina-Par´ıs, “Stochastic competitive exclusionin the maintenance of the na¨ıve T cell repertoire,”
J. Theor. Biol. , vol. 265, no. 3, pp. 396–410.,2010.[38] J. Reynolds, M. Coles, G. Lythe, and C. Molina-Par´ıs, “Deterministic and stochastic naive T cellpopulation dynamics: symmetric and asymmetric cell division,”
Dyn. Syst. , vol. 27, pp. 75–103.,2012.[39] D. Alonso, A. J. McKane, and M. Pascual, “Stochastic amplification in epidemics,”
J. R. Soc.Interface , vol. 4, no. 14, pp. 575–582., 2007.[40] R. Kuske, L. F. Gordillo, and P. E. Greenwood, “Sustained oscillations via coherence resonancein SIR,”
J. Theo. Biol. , vol. 245, no. 3, pp. 459–469., 2007.[41] K. B. Blyuss and L. B. Nicholson, “The role of tunable activation thresholds in the dynamics ofautoimmunity,”
J. Theor. Biol. , vol. 308, pp. 45–55., 2012.[42] K. B. Blyuss and L. B. Nicholson, “Understanding the roles of activation threshold and infectionsin the dynamics of autoimmune disease,”
J. Theor. Biol. , vol. 375, pp. 13–20., 2015.[43] D. Ben Ezra and J. V. Forrester, “Fundal white dots: the spectrum of a similar pathologicalprocess,”
Brit. J. Ophthalmol. , vol. 79, pp. 856–860., 1995.[44] T. F. Davies, D. C. Evered, B. Rees Smith, P. P. B. Yeo, F. Clark, and R. Hall, “Value ofthyroid-stimulating-antibody determination in predicting the short-term thyrotoxic relapse ingraves’ disease,”
Lancet , vol. 309, pp. 1181–1182., 1997.[45] A. Nylander and D. A. Hafler, “Multiple sclerosis,”
J. Clin. Invest. , vol. 122, pp. 1180–1188.,2012.[46] S. D. Wolf, B. N. Dittel, F. Hardardottir, and C. A. Janeway Jr, “Experimental autoimmuneencephalomyelitis induction in genetically B cell-deficient mice,”
J. Exp. Med. , vol. 184, pp. 2271–2278., 1996.[47] H.-J. Wu, I. I. Ivanov, J. Darceet, K. Hattori, T. Shima, Y. Umesaki, and et al, “Gut-residingsegmented filamentous bacteria drive autoimmune arthritis via T helper 17 cells,”
Immunity ,vol. 32, pp. 815–827., 2010.[48] L. J. S. Allen,
An introduction to stochastic processes with applications to biology.
Mathematical Epidemiology.
Appl. Math. Model. , vol. 38, no. 5, pp. 1583–1596., 2014.2051] D. T. Gillespie, “The chemical Langevin equation,”
J. Chem. Phys. , vol. 113, no. 1, pp. 297–306.,2000.[52] D. T. Gillespie, “The chemical Langevin and Fokker-Planck equations for the reversible isomer-ization reaction,”
J. Phys. Chem. A , vol. 106, no. 20, pp. 5063–5071., 2002.[53] B. Øksendal,
Stochastic Differential Equations: An Introduction with Applications.
Math. Biosci. ,vol. 206, pp. 81–107., 2007.[55] E. J. Allen,
Modeling with Itˆo stochastic differential equations . Springer, 2007. Berlin.[56] E. J. Allen, L. J. S. Allen, A. Arciniega, and P. E. Greenwood, “Construction of equivalentstochastic differential equation models,”
Stoch. Anal. Appl. , vol. 26, no. 2, pp. 274–297., 2008.[57] N. G. van Kampen,
Stochastic processes in physics and chemistry.
Handbook of stochastic methods for physics, chemistry and the natural sciences.
Berlin: Springer-Verlag., 2004.[59] A. J. Black, A. J. McKane, A. Nunes, and A. Parisi, “Stochastic fluctuations in the susceptible-infective-recovered model with distributed infectious periods,”
Phys. Rev. E , vol. 80, no. 2,p. 021922., 2009.[60] E. W. J. Wallace, D. T. Gillespie, K. R. Sanft, and L. R. Petzold, “Linear noise approximation isvalid over limited times for any chemical system that is sufficiently large,”
IET Syst. Biol. , vol. 6,no. 4, pp. 102–115., 2012.[61] F. Hayot and C. Jayaprakash, “The linear noise approximation for molecular fluctuations withincells,”
Phys. Biol. , vol. 1, pp. 205–210., 2004.[62] J. Pahle, J. D. Challenger, P. Mendes, and A. J. McKane, “Biochemical fluctuations, optimisationand the linear noise approximation,”
BMC Syst. Biol. , vol. 6, p. 86., 2012.[63] J. M. Conway and D. Coombs, “A stochastic model of latently infected cell reactivation and viralblip generation in treated HIV patients,”
PLoS Comp. Biol. , vol. 7, no. 4, p. e1002033., 2011.[64] M. Bruna, S. J. Chapman, and M. J. Smith, “Model reduction for slow-fast stochastic systemswith metastable behaviour,”
J. Chem. Phys. , vol. 140, p. 174107., 2014.[65] T. M. Earnest, E. Roberts, M. Assaf, D. K., and Z. Luthey-Schulten, “DNA looping increasesthe range of bistability in a stochastic model of the lac genetic switch,”
Phys. Biol. , vol. 10,p. 026002., 2013.[66] P. G. Hufton, Y. T. Lin, T. Galla, and A. J. McKane, “Intrinsic noise in systems with switchingenvironments,”
Phys. Rev. E , vol. 93, p. 052119., 2016.[67] A. J. Black and A. J. McKane, “Stochastic formulation of ecological models and their applica-tions,”
Trends Ecol. Evol. , vol. 27, no. 6, pp. 337–345., 2012.[68] G. Neofytou, Y. N. Kyrychko, and K. B. Blyuss, “Time-delayed model of RNA interference,”
Ecol. Complex. , vol. 30, pp. 11–25., 2017. 2169] R. H. Bartels and G. W. Stewart, “Solution of the matrix equation AX+XB=C,”
Comm. ACM ,vol. 15, pp. 820–826., 1972.[70] S. J. Hammarling, “Numerical solution of the stable, non-negative definite lyapunov equation,”
IMA J. Num. Anal. , vol. 2, pp. 303–323., 1982.[71] N. J. Burroughs, B. M. P. M. Oliveira, and A. A. Pinto, “Regulatory T cell adjustment of quorumgrowth thresholds and the control of local immune responses,”
J. Theor. Biol. , vol. 241, pp. 134–141., 2006.[72] A. Skapenko, J. Leipe, P. E. Lipsky, and H. Schulze-Koops, “The role of the T cell in autoimmuneinflammation,”
Arthr. Res. Ther. , vol. 7(Suppl 2), pp. S4–S14., 2005.[73] R. Antia, V. V. Ganusov, and R. Ahmed, “The role of models in understanding CD8 + T-cellmemory,”
Nature Rev. Immunol. , vol. 5, pp. 101–111., 2005.[74] Y. Yuan and L. J. S. Allen, “Stochastic models for virus and immune system dynamics,”
Math.Biosci. , vol. 234, no. 2, pp. 84–94., 2011.[75] A. J. Black and A. J. McKane, “Stochastic amplification in an epidemic model with seasonalforcing,”
J. Theor. Biol. , vol. 267, no. 1, pp. 85–94., 2010.[76] A. J. McKane and T. J. Newman, “Predator-prey cycles from resonant amplification of demo-graphic stochasticity,”
Phys. Rev. Lett. , vol. 94, no. 21, p. 218102., 2005.[77] J. Boldison, T. K. Khera, D. A. Copland, M. L. Stimpson, G. L. Crawford, A. D. Dick, andL. B. Nicholson, “A novel pathogenic RBP-3 peptide reveals epitope spreading in persistentexperimental autoimmune uveoretinitis,”
Immunol. , vol. 146, pp. 301–311., 2015.22
Supplementary Material
Derivation of the van Kampen’s system size expansion
As described in Section 2.2 of the main text, the CTMC model based on the transitions probabilitiesyields the following master equation: dP ( n , t ) dt = { ( ε − − q + ( ε +1 − q + ( ε +1 ε − − q + ( ε − − q + ( ε +3 − q + ( ε +3 ε − − q + ( ε +3 ε − − q + ( ε +3 ε − − q + ( ε +2 − q + ( ε − − q + ( ε +4 − q + ( ε − − q + ( ε +5 − q + ( ε − − q + ( ε +6 − q + ( ε − − q + ( ε +7 − q } P ( n , t ) , (15)where n = ( n , n , n , n , n , n , n ) is the current state of the system, the coefficients q i are given by q = b n + b n , q = d n + d n + µ a n n , q = βn n , q = λ in ,q = d in n , q = p αn n , q = p αn n , q = (1 − p − p ) αn n ,q = ( d F + µ F n + µ a n ) n , q = λ r + ρ n n , q = d r n , q = ρ n n ,q = d n n , q = ρ n n , q = ( d a + δn ) n , q = σ n + σ n , q = d i n , and the operators ε ± i are defined as follows, ε ± i f ( n , n , n , n , n , n , n , t ) = f ( n , ..., n i ± , ..., n , t ) , ≤ i ≤ . If n i < ≤ i ≤
7, then P ( n , t ) = 0.To derive the van Kampen’s system size expansion of the master equation (3), we rewrite each n i ( t ) in the form n i ( t ) = Ω x i ( t ) + Ω / ζ i ( t ) , ≤ i ≤ , where Ω x i ( t ) = E [ n i ( t )], so ζ i ( t ) represent the fluctuations. Replacing the probability density P ( n , t )by the equivalent probability density Π( ζ , t ), i.e. with Π( ζ , t ) = P ( n , t ) = P (cid:0) Ω x + Ω / ζ, t (cid:1) , theleft-hand side of the master equation (3) transforms into dP ( n , t ) dt = ∂ Π ∂t − (cid:88) i =1 Ω / dx i dt ∂ Π ∂ζ i . (16)The operators ε ± i and their product now satisfy the following expansions ε ± i = 1 ± Ω − / ∂∂ζ i + 12 Ω − ∂ ∂ζ i ± · · · ,ε + i ε − j = (cid:18) − / ∂∂ζ i + 12 Ω − ∂ ∂ζ i + · · · (cid:19) (cid:32) − Ω − / ∂∂ζ j + 12 Ω − ∂ ∂ζ j − · · · (cid:33) = 1 + Ω − / (cid:18) ∂∂ζ i − ∂∂ζ j (cid:19) + Ω − (cid:32) ∂ ∂ζ i − ∂ ∂ζ i ∂ζ j + 12 ∂ ∂ζ j (cid:33) + · · · , ≤ i, j ≤ . (17)23ne can easily obtain/introduce Ω-expansions for all parameters q i : q = b n + b n = b (cid:16) Ω x + Ω / ζ (cid:17) + b (cid:16) Ω x + Ω ζ + 2Ω / x ζ (cid:17) = b (cid:16) Ω x + Ω / ζ (cid:17) + b Ω (cid:124)(cid:123)(cid:122)(cid:125) (cid:101) b (cid:16) Ω x + ζ + 2Ω / x ζ (cid:17) = (cid:101) b ζ + (cid:16) b ζ + 2 (cid:101) b x ζ (cid:17) Ω / + (cid:16) b x + (cid:101) b x (cid:17) Ω , q = d n + d n + µ a n n = d (cid:16) Ω x + Ω / ζ (cid:17) + d (cid:16) Ω x + Ω ζ + 2Ω / x ζ (cid:17) + µ a (cid:16) Ω x + Ω / ζ (cid:17) (cid:16) Ω x + Ω / ζ (cid:17) = d (cid:16) Ω x + Ω / ζ (cid:17) + d Ω (cid:124)(cid:123)(cid:122)(cid:125) (cid:101) d (cid:16) Ω x + ζ + 2Ω / x ζ (cid:17) + µ a Ω (cid:124)(cid:123)(cid:122)(cid:125) (cid:101) µ a (cid:16) Ω / x + ζ (cid:17) (cid:16) Ω / x + ζ (cid:17) = d (cid:16) Ω x + Ω / ζ (cid:17) + (cid:101) d (cid:16) Ω x + ζ + 2Ω / x ζ (cid:17) + (cid:101) µ a (cid:16) Ω x x + Ω / x ζ + Ω / x ζ + ζ ζ (cid:17) = (cid:101) µ a ζ ζ + (cid:16) d ζ + 2 (cid:101) d x ζ + (cid:101) µ a x ζ + (cid:101) µ a x ζ (cid:17) Ω / + (cid:16) d x + (cid:101) d x + (cid:101) µ a x x (cid:17) Ω , q = βn n = β (cid:16) Ω x + Ω / ζ (cid:17) (cid:16) Ω x + Ω / ζ (cid:17) = β Ω (cid:124)(cid:123)(cid:122)(cid:125) (cid:101) β (cid:16) Ω / x + ζ (cid:17) (cid:16) Ω / x + ζ (cid:17) = (cid:101) β (cid:16) Ω x x + Ω / x ζ + Ω / x ζ + ζ ζ (cid:17) = (cid:101) βζ ζ + (cid:16) (cid:101) βx ζ + (cid:101) βx ζ (cid:17) Ω / + (cid:101) βx x Ω ,q = λ in = λ in Ω (cid:124)(cid:123)(cid:122)(cid:125) (cid:101) λ in Ω = (cid:101) λ in Ω ,q = d in n = d in (cid:16) Ω x + Ω / ζ (cid:17) = d in ζ Ω / + d in x Ω , In a similar way we can easily show q = p (cid:101) αζ ζ + ( p (cid:101) αx ζ + p (cid:101) αx ζ ) Ω / + p (cid:101) αx x Ω ,q = p (cid:101) αζ ζ + ( p (cid:101) αx ζ + p (cid:101) αx ζ ) Ω / + p (cid:101) αx x Ω ,q = (1 − p − p ) (cid:104)(cid:101) αζ ζ + (cid:101) αx ζ + (cid:101) αx ζ Ω / + (cid:101) αx x Ω (cid:105) ,q = ( (cid:101) µ F ζ ζ + (cid:101) µ a ζ ζ ) + ( d F + (cid:101) µ F x ζ + (cid:101) µ F x ζ + (cid:101) µ a x ζ + (cid:101) µ a x ζ ) Ω / + ( d F x + (cid:101) µ F x x + (cid:101) µ a x x ) Ω ,q = (cid:101) ρ ζ ζ + ( (cid:101) ρ x ζ + (cid:101) ρ x ζ ) Ω / + (cid:16)(cid:101) λ r + (cid:101) ρ x x (cid:17) Ω , = d r ζ Ω / + d r x Ω , q = (cid:101) ρ ζ ζ + ( (cid:101) ρ x ζ + (cid:101) ρ x ζ ) Ω / + (cid:101) ρ x x Ω ,q = d n ζ Ω / + d n x Ω , q = (cid:101) ρ ζ ζ + ( (cid:101) ρ x ζ + (cid:101) ρ x ζ ) Ω / + (cid:101) ρ x x Ω ,q = (cid:101) δζ ζ + (cid:16) d a ζ + (cid:101) δx ζ + (cid:101) δx ζ (cid:17) Ω / + (cid:16) d a x + (cid:101) δx x (cid:17) Ω ,q = ( σ ζ + σ ζ ) Ω / + ( σ x + σ x )Ω , q = d i ζ Ω / + d i x Ω , where λ r = (cid:101) λ r Ω , µ F = (cid:101) µ F Ω , α = (cid:101) α Ω , δ = (cid:101) δ Ω , ρ i = (cid:101) ρ i Ω , i = 1 , , . Substituting expressions (16), (17) and q i ’s into the master equation (3) shows that the left-hand sideof the equation only contains terms of the order Ω / and Ω , while the right-hand side has terms ofthe order Ω / , Ω , and Ω − n/ , for n ∈ N . To derive a linear Fokker-Planck equation, we ignore theterms of order Ω − n/ , for n ∈ N .Considering the terms of order Ω / , i.e. only the terms that are proportional to ∂ Π /∂ζ i , yields − Ω / dx dt ∂ Π ∂ζ = (cid:18) − Ω − / ∂∂ζ (cid:19) (cid:104)(cid:16) b x + (cid:101) b x (cid:17) Ω (cid:105) Π + (cid:18) Ω − / ∂∂ζ (cid:19) (cid:104)(cid:16) d x + (cid:101) d x + (cid:101) µ a x x (cid:17) Ω (cid:105) Π+ (cid:18) Ω − / ∂∂ζ (cid:19) (cid:104) (cid:101) βx x Ω (cid:105) Π , − Ω / dx dt ∂ Π ∂ζ = (cid:18) − Ω − / ∂∂ζ (cid:19) (cid:104) (cid:101) βx x Ω (cid:105) Π + (cid:18) Ω − / ∂∂ζ (cid:19) [( d F x + (cid:101) µ F x x + (cid:101) µ a x x ) Ω] Π , − Ω / dx dt ∂ Π ∂ζ = (cid:18) − Ω − / ∂∂ζ (cid:19) (cid:104)(cid:101) λ in Ω (cid:105) Π + (cid:18) Ω − / ∂∂ζ (cid:19) [ d in x Ω] Π + (cid:18) Ω − / ∂∂ζ (cid:19) [ p (cid:101) αx x Ω] Π (cid:18) Ω − / ∂∂ζ (cid:19) [ p (cid:101) αx x Ω] Π + (cid:18) Ω − / ∂∂ζ (cid:19) [(1 − p − p ) (cid:101) αx x Ω] Π , − Ω / dx dt ∂ Π ∂ζ = (cid:18) − Ω − / ∂∂ζ (cid:19) [ p (cid:101) αx x Ω] Π + (cid:18) − Ω − / ∂∂ζ (cid:19) (cid:104)(cid:16)(cid:101) λ r + (cid:101) ρ x x (cid:17) Ω (cid:105) Π+ (cid:18) Ω − / ∂∂ζ (cid:19) [ d r x Ω] Π , − Ω / dx dt ∂ Π ∂ζ = (cid:18) − Ω − / ∂∂ζ (cid:19) [ p (cid:101) αx x Ω] Π + (cid:18) − Ω − / ∂∂ζ (cid:19) [ (cid:101) ρ x x Ω] Π − (cid:18) Ω − / ∂∂ζ (cid:19) [ d n x Ω] Π , − Ω / dx dt ∂ Π ∂ζ = (cid:18) − Ω − / ∂∂ζ (cid:19) [(1 − p − p ) (cid:101) αx x Ω] Π + (cid:18) − Ω − / ∂∂ζ (cid:19) [ (cid:101) ρ x x Ω] Π+ (cid:18) Ω − / ∂∂ζ (cid:19) (cid:104)(cid:16) d a x + (cid:101) δx x (cid:17) Ω (cid:105) Π , − Ω / dx dt ∂ Π ∂ζ = (cid:18) − Ω − / ∂∂ζ (cid:19) [( σ x + σ x )Ω] Π + (cid:18) Ω − / ∂∂ζ (cid:19) [ d i x Ω] Π . dx dt = b x + (cid:101) b x − d x − (cid:101) d x − (cid:101) βx x − (cid:101) µ a x x ,dx dt = (cid:101) βx x − d F x − (cid:101) µ F x x − (cid:101) µ a x x ,dx dt = (cid:101) λ in − d in x − (cid:101) αx x ,dx dt = (cid:101) λ r − d r x + p (cid:101) αx x + (cid:101) ρ x x ,dx dt = p (cid:101) αx x − d n x + (cid:101) ρ x x ,dx dt = (1 − p − p ) (cid:101) αx x − d a x − (cid:101) δx x + (cid:101) ρ x x ,dx dt = σ x + σ x − d i x . (18)Terms of order Ω give the following Fokker-Planck equation ∂ Π ∂t = − (cid:20) (cid:16) b + 2 (cid:101) b x − d − (cid:101) d x − (cid:101) µ a x − (cid:101) βx (cid:17) ∂ ( ζ Π) ∂ζ − (cid:101) βx ∂ ( ζ Π) ∂ζ − (cid:101) µ a x ∂ ( ζ Π) ∂ζ + (cid:101) βx ∂ ( ζ Π) ∂ζ + (cid:16) (cid:101) βx − d F − (cid:101) µ F x − (cid:101) µ a x (cid:17) ∂ ( ζ Π) ∂ζ − (cid:101) µ F x ∂ ( ζ Π) ∂ζ − (cid:101) µ a x ∂ ( ζ Π) ∂ζ − (cid:101) αx ∂ ( ζ Π) ∂ζ − ( d in + (cid:101) αx ) ∂ ( ζ Π) ∂ζ + p (cid:101) αx ∂ ( ζ Π) ∂ζ + p (cid:101) αx ∂ ( ζ Π) ∂ζ + ( (cid:101) ρ x − d r ) ∂ ( ζ Π) ∂ζ + (cid:101) ρ x ∂ ( ζ Π) ∂ζ + p (cid:101) αx ∂ ( ζ Π) ∂ζ + p (cid:101) αx ∂ ( ζ Π) ∂ζ + ( (cid:101) ρ x − d n ) ∂ ( ζ Π) ∂ζ + (cid:101) ρ x ∂ ( ζ Π) ∂ζ + (1 − p − p ) (cid:101) αx ∂ ( ζ Π) ∂ζ + (1 − p − p ) (cid:101) αx ∂ ( ζ Π) ∂ζ − (cid:101) δx ∂ ( ζ Π) ∂ζ + (cid:16)(cid:101) ρ x − d a − (cid:101) δ (cid:17) x ∂ ( ζ Π) ∂ζ + (cid:101) ρ x ∂ ( ζ Π) ∂ζ + σ ∂ ( ζ Π) ∂ζ + σ ∂ ( ζ Π) ∂ζ − d i ∂ ( ζ Π) ∂ζ (cid:21) + 12 (cid:26) (cid:16) b x + (cid:101) b x + d x + (cid:101) d x + (cid:101) βx x + (cid:101) µ a x x (cid:17) ∂ Π ∂ζ − (cid:101) βx x ∂ Π ∂ζ ∂ζ + (cid:16) (cid:101) βx x + d F x + (cid:101) µ F x x + (cid:101) µ a x x (cid:17) ∂ Π ∂ζ + (cid:16)(cid:101) λ in + d in x + (cid:101) αx x (cid:17) ∂ Π ∂ζ − p (cid:101) αx x ∂ Π ∂ζ ∂ζ − p (cid:101) αx x ∂ Π ∂ζ ∂ζ − − p − p ) (cid:101) αx x ∂ Π ∂ζ ∂ζ + (cid:16)(cid:101) λ r + d r x + p (cid:101) αx x + (cid:101) ρ x x (cid:17) ∂ Π ∂ζ + ( p (cid:101) αx x + d n x + (cid:101) ρ x x ) ∂ Π ∂ζ + (cid:104) (1 − p − p ) (cid:101) αx x + d a x + (cid:101) δx x + (cid:101) ρ x x (cid:105) ∂ Π ∂ζ + ( σ x + σ x + d i x ) ∂ Π ∂ζ (cid:27) . This equation can be equivalently rewritten in the form ∂ Π( ζ , t ) ∂t = − (cid:88) i,j A ij ∂∂ζ i ( ζ j Π) + 12 (cid:88) i,j B ij ∂ Π ∂ζ i ∂ζ j , A is the Jacobian matrix of system (18) A = b + 2 (cid:101) b x − d − (cid:101) d x − (cid:101) µ a x − (cid:101) βx − (cid:101) βx − (cid:101) µ a x (cid:101) βx (cid:101) βx − d F − (cid:101) µ F x − (cid:101) µ a x − (cid:101) µ F x − (cid:101) µ a x − (cid:101) αx − d in − (cid:101) αx p (cid:101) αx p (cid:101) αx (cid:101) ρ x − d r (cid:101) ρ x p (cid:101) αx p (cid:101) αx (cid:101) ρ x − d n (cid:101) ρ x − p − p ) (cid:101) αx (1 − p − p ) (cid:101) αx − (cid:101) δx (cid:101) ρ x − d a − (cid:101) δx (cid:101) ρ x σ σ − d i B is a 7 × B ij = b x + (cid:101) b x + d x + (cid:101) d x + (cid:101) βx x + (cid:101) µ a x x , if ( i, j ) = (1 , , (cid:101) βx x + d F x + (cid:101) µ F x x + (cid:101) µ a x x , if ( i, j ) = (2 , , (cid:101) λ in + d in x + (cid:101) αx x , if ( i, j ) = (3 , , (cid:101) λ r + d r x + p (cid:101) αx x + (cid:101) ρ x x , if ( i, j ) = (4 , ,p (cid:101) αx x + d n x + (cid:101) ρ x x , if ( i, j ) = (5 , , (1 − p − p ) (cid:101) αx x + d a x + (cid:101) δx x + (cid:101) ρ x x , if ( i, j ) = (6 , ,σ x + σ x + d i x , if ( i, j ) = (7 , , − (cid:101) βx x , if ( i, j ) = (1 ,
2) or (2 , , − p (cid:101) αx x , if ( i, j ) = (3 ,
4) or (4 , , − p (cid:101) αx x , if ( i, j ) = (3 ,
5) or (5 , , − (1 − p − p ) (cid:101) αx x , if ( i, j ) = (3 ,
6) or (6 , ,,
6) or (6 , ,, ,,