Effects of viral and cytokine delays on dynamics of autoimmunity

A major contribution to the onset and development of autoimmune disease is known to come from infections. An important practical problem is identifying the precise mechanism by which the breakdown of immune tolerance as a result of immune response to infection leads to autoimmunity. In this paper, we develop a mathematical model of immune response to a viral infection, which includes T cells with different activation thresholds, regulatory T cells (Tregs), and~a cytokine mediating immune dynamics. Particular emphasis is made on the role of time delays associated with the processes of infection and mounting the immune response. Stability analysis of various steady states of the model allows us to identify parameter regions associated with different types of immune behaviour, such as, normal clearance of infection, chronic infection, and autoimmune dynamics. Numerical simulations are used to illustrate different dynamical regimes, and to identify basins of attraction of different dynamical states. An important result of the analysis is that not only the parameters of the system, but also the initial level of infection and the initial state of the immune system determine the progress and outcome of the dynamics.


Introduction
An immune system can only be viewed as effective when it can robustly identify and destroy pathogen-infected cells, while also distinguishing such cells from healthy cells. In the case of breakdown of immune tolerance, the immune system fails to discriminate between self-antigens and foreign antigens, which results in autoimmune disease, i.e., undesired destruction of healthy cells. Under normal conditions, once foreign epitopes are presented on antigen presenting cells (APCs) to T cells, this results in the proliferation of T cells and eliciting their effector function. While this mechanism is responsible for a successful clearance of infection, cross-reactivity between epitopes associated with foreign and self-antigens can lead to a T cell response against healthy host cells [1,2].
For many autoimmune diseases, the disease occurs in a specific organ or part of the body, such as retina in uveitis, central nervous system in multiple sclerosis, or pancreatic β-cells in type-1 diabetes [3,4,5]. It is extremely difficult to identify the specific causes of autoimmunity in individual patients, as it usually has contributions from a number of internal and external factors, including a genetic predisposition, age, previous immune challenges, exposure to pathogens etc., [6,7,8,9]. Even though genetic predisposition is known to play a very significant role, it is believed that some additional environmental triggers are required for the onset of autoimmunity, and these are usually represented by infections [10,11]. A very recent work has experimentally identified a gut bacterium that, when present in mice and humans, can migrate to other parts of the body, facilitating subsequent triggering of autoimmune disease in those organs [12]. Various mechanisms of onset of pathogen-induced autoimmune disease have been identified, including bystander activation [13] and molecular mimicry [14,15], which is particularly important in the context of autoimmunity caused by viral infections.
A number of mathematical models have looked into dynamics of onset and development of autoimmune disease. Segel et al. [16] analysed interactions between effector and regulatory T cells in the context of T cell vaccination, without explicitly specifying possible causes of autoimmunity. Similar models were later studied by Borghans et al. [17,18] who demonstrated possible onset of autoimmune state, defined as stable above-threshold oscillations in the number of autoreactive cells, as a result of interactions between regulatory and autoreactive T cells. León et al. [19,20,21] and Carneiro et al. [22] have studied interactions between different T cells, with an emphasis on the suppressing role of regulatory T cells in the dynamics of immune response and control of autoimmunity. Alexander and Wahl [23] have also looked into the role of regulatory T cells, in particular focusing on their interactions with professional APCs and effector cells for the purpose of controlling immune response. Iwami et al. [24,25] explicitly included a separate compartment representing the viral population in their models of immune response, and showed that the functional form of the growth function for susceptible host cells can have a significant effect on the resulting immune dynamics. Despite being able to explain the emergence of autoimmunity as a byproduct of immune response to infection, these models were not able to exhibit another practically important dynamical regime of normal viral clearance. For the case of pathogen-induced autoimmunity arising through bystander activation, Burroughs et al. [26,27,28] have developed a model that investigates interactions between T cells and interleukin-2 (IL-2), an important cytokine, in mediating the onset of autoimmunity.
Among various parts of the immune system involved in coordinating an effective immune response, a particularly significant role is known to be played by the T cells, with experimental evidence suggesting that regulatory T cells are vitally important for controlling autoimmunity [29,30,31,32]. To account for this fact in mathematical models, Alexander and Wahl [23] and Burroughs et al. [26,27] have explicitly included a separate compartment for regulatory T cells that are activated by autoantigens and suppress the activity of autoreactive T cells. Another framework for modelling the effects of T cells on autoimmune dynamics is by using the idea that T cells have different or tunable activation thresholds (TAT), which result in different immune functionality of the same T cells, and also allow T cells to adjust their response to simulation by autoantigens. This approach was proposed for the analysis of the peripheral and central T cell dynamics [33,34,35], it has also been used to study differences in activation/response thresholds that are dependent on the activation state of the T cell [36]. Murine and human experiments have confirmed that activation of T cells can indeed change dynamically during their circulation [37,38,39,40]. To model this feature, van den Berg and Rand [41], and Scherer et al. [42] developed stochastic models for the tuning of activation thresholds.
Blyuss and Nicholson [43,44] have proposed a mathematical model of autoimmunity resulting from immune response to a viral infection through a mechanism of molecular mimicry. This model explicitly includes the virus population and two types of T cells with different activation thresholds, and it also accounts for a biologically realistic scenario where infection and autoimmune response can occur in different organs of the host. Besides the normal viral clearance and chronic infection, in some parameter regime the model also exhibits an autoimmune state characterised by stable oscillations in the amounts of cell populations. From a clinical perspective, such behaviour is to be expected, as it is associated with relapses and remissions that have been observed in a number of autoimmune diseases, such as autoimmune thyroid disease, MS, and uveitis [45,46,47]. One deficiency of this model is the fact that the oscillations associated with autoimmune regime can only occur if the amount of free virus and the number of infected cells are also exhibiting oscillations, while in clinical and laboratory settings, autoimmunity usually occurs after the initial infection has been fully cleared. To overcome this limitation, Fatehi et al. [48] have recently developed a more advanced model that also includes regulatory T cells and cytokines, which has allowed the authors to obtain a more realistic representation of immune response and various dynamical regimes. A particularly important practical insight provided by this model is the observation that it is not only the system parameters, but also the initial level of infection and the initial state of the immune system, that determine whether the host will just successfully clear the infection, or will proceed to develop autoimmunity. Approaching the same problem from another perspective, Fatehi et al. [49] have investigated the role of stochasticity in driving the dynamics of immune response and determining which of the immune states is more likely to be attained. The authors have also determined an experimentally important characterisation of autoimmune state, as provided by the dependence of variance in cell populations on various system parameters.
In this paper, we develop and analyse a model of autoimmune dynamics, with particular focus on the role of time delays associated with different aspects of immune response, as well as an inhibiting effect of regulatory T cells on secretion of IL-2. In the next section, we introduce the model and discuss its basic properties. Section 3 contains systematic analysis of all steady states, including conditions for their feasibility and stability. Section 4 contains a bifurcation analysis of the model and demonstrates various types of behaviour that the system exhibits depending on parameters and initial conditions, which includes identification of attraction basins of various states. The paper concludes in Section 5 with the discussion of results.

Model Derivation
To understand how interactions between different parts of the immune system and the cytokine can lead to autoimmunity, we consider a model illustrated in a diagram shown in Figure 1. In this model, unlike earlier work of Blyuss and Nicholson [43,44], we consider a situation where a viral infection and possible autoimmunity occur in the same organ of the host. The healthy host cells, whose number is denoted by A(t), in the absence of infection are assumed to grow logistically with linear growth rate r and carrying capacity N , and they acquire infection at rate β from the infected cells F (t). Since experimental evidence suggests that antibodies play a secondary role compared to T cells [50], and autoimmunity can develop even in the absence of B cells [51], we do not include antibody response in the model, but focus solely on the dynamics of T cell populations. Naïve (inactivated) T cells T in (t) are assumed to be in homeostasis [43,24,25], and once they are activated through interaction with infected cells, which occurs at rate α, a proportion p 1 of them will go on to differentiate into additional regulatory T cells, a fraction p 2 will become normal activated T cells T nor (t) able to destroy infected cells at rate µ F upon recognition of foreign antigen present on their surface. The remaining proportion of (1 − p 1 − p 2 ) of T cells will become autoreactive T cells T aut (t) that, in light of their lower activation threshold will be eliminating both infected cells and healthy host cells at rate µ a due to the above-mentioned cross-reactivity between some epitopes in self-and foreign antigens. Regulatory T cells T reg (t) are assumed to be in their own homeostasis [52], and their main contribution to immune dynamics lies in suppressing autoreactive T cells at rate δ 1 . To reduce the dimensionality of the model, it is assumed that the process of viral production is occurring very fast compared to other characteristic timescales of the model, thus the viral population can be represented by its quasi-steady-state approximation, i.e., it is taken to be proportional to the number of infected cells, and this eliminates the need for a separate compartment for free virus.
A number of different cytokines mediate immune response to infection, and in the context of T cell dynamics, a particular important role is played by interleukin 2 (IL-2), represented by the variable I(t) in the model, which acts to enhance the proliferation of T cells, which, in turn, secrete further IL-2. One of the actions of regulatory T cells is to suppress the expression of IL-2 [53], which is only produced by the activated T cells, but not by the regulatory T cells [54,55]. To represent this mathematically, we will assume that T nor and T aut produce IL-2 at rates σ 1 and σ 2 , and conversely, IL-2 enhances proliferation of T reg , T nor and T aut at rates ρ 1 , ρ 2 , and ρ 3 . We include in the model suppression of IL-2 by regulatory T cells at rate δ 2 , in a manner similar to Burroughs et al. [28].
While the production of new virus particles by infected cells is assumed to be fast, we explicitly include in the model time delay τ 1 associated with the actual process of infection, which includes multiple stages of the eclipse phase of viral life cycle, such as virus attachment, cell penetration and uncoating [56,57]. We also include the time delay τ 2 associated with simulation and proliferation of T cells by IL-2, and the time delay τ 3 between antigen encounter and resulting T cell expansion [58].
With the above assumptions, the complete model takes the form Introducing non-dimensional variableŝ yields a rescaled model where all hats in variables and parameters have been dropped for simplicity of notation, and all parameters are assumed to be positive. It is easy to show that this system is well-posed, i.e., solutions with non-negative initial conditions remain non-negative for all t ≥ 0. As a first step in the analysis of model (1), we look at its steady states that can be found by equating to zero the right-hand sides of Equation (1) and solving the resulting system of algebraic equations, deferring the discussion of conditionally stable steady states to Section 3. First, we consider a situation where there are no infected cells at a steady state, i.e., F * = 0, which immediately implies T * in = 1. In this case, there are four possible combinations of steady states depending on whether T * nor and T * aut are each equal to zero or positive. If T * nor = T * aut = 0, there are two steady states of which S * 1 is always unstable, and S * 2 is a disease-free conditionally stable steady state. For T * nor = 0 and T * aut = 0, we again have two steady states , but they are both unstable for any values of parameters. In the case when T * nor = 0 and T * aut = 0, we have two further steady states S * 5 and S * 6 , The steady state S * 5 has A * = 0, which implies the death of host cells, whereas the steady state S * 6 corresponds to an autoimmune regime. The steady state S * 7 with T * nor = 0 and T * aut = 0 exists only for a particular combination of parameters, namely, when and is always unstable. Finally, when F * = 0, the system (1) can have a steady state S * 8 with all of its components being positive, but it does not appear possible to find a closed form expression for this state.
In summary, besides the unconditionally unstable steady states, the model (1) has at most for conditionally stable steady states: the disease-free steady state S * 2 , the steady state with the death of host cells S * 5 , the autoimmune steady state S * 6 , and the persistent or chronic steady state S * 8 .

Stability Analysis of the Disease-Free Steady State
Linearising the system (1) near the disease-free steady state S * 2 yields the following equation for If d F < β, the above equation always has a real positive root for any value τ 1 ≥ 0, implying that the disease-free steady state is always unstable for any value of the time delays. If, however, the condition d F > β holds, the disease-free steady state is stable for τ 1 = 0. To find out whether it can lose stability for τ 1 > 0, we look for solutions of Equation (2) in the form λ = iω. Separating real and imaginary parts yields Squaring and adding these two equations gives the following equation for potential Hopf frequency ω ω 2 + d 2 F − β 2 = 0. since d F > β, this equation does not have real roots for ω, suggesting that there can be no roots of the form λ = iω of the characteristic Equation (2). This implies that in the case d F > β the disease-free steady state S * 2 is stable for all values of the time delay τ 1 ≥ 0.

Stability Analysis of the Death, Autoimmune and Chronic Steady States
The steady state S * 5 (respectively, S * 6 ) is stable if and all roots of the following equation have negative real part where This steady state undergoes a steady-state bifurcation if For τ 2 = 0 these steady states are stable if T * reg satisfies (3) and where To investigate whether stability can be lost for τ 2 > 0, we use an iterative procedure described in [59,60] to determine a function F (ω), whose roots give the Hopf frequency associated with purely imaginary roots of Equation (4). Substituting λ = iω into Equation (4), we define ∆ (1) (τ 2 , λ) as and the bar denotes the complex conjugate. If we define The explicit formulae for other coefficients of F (ω) can be found in Appendix 5. Introducing s = ω 2 , the equation F (ω) = 0 can be equivalently rewritten as follows, Without loss of generality, suppose that Equation (7) has six distinct positive roots denoted by s 1 , s 2 , ... , s 6 , which means that the equation F (ω) = 0 has six positive roots Substituting λ k = iω k into Equation (4) gives This allows us to find as the first time delay for which the roots of the characteristic Equation (4) cross the imaginary axis. To determine whether these steady states actually undergo a Hopf bifurcation at τ 2 = τ * , we have to compute the sign of dRe[λ(τ * )]/dτ 2 . For τ = τ * , λ(τ * ) = iω 0 , and we also define s 0 = ω 2 0 .
0 (iω 0 ) = 0. Then the following transversality condition holds where all x j and y j are expressed in terms of system parameters and steady state values of the variables. Substituting these expressions into ∆(τ 2 , iω 0 ) = 0 and ∆ (1) (τ 2 , iω 0 ) = 0, and then separating real and imaginary parts gives Solving this system of equations provides the values of sin(ω 0 τ * ), cos(ω 0 τ * ), sin(2ω 0 τ * ), and cos(2ω 0 τ * ). Taking the derivative of Equation (4) with respect to τ 2 , one finds Hence,  Table  1 (a), and with µ a = 10 (b). Black and red curves indicate the boundaries of feasibility and the steady-state bifurcation, whereas dashed lines (blue/brown) show the boundaries of Hopf bifurcation of the steady states S * 5 and S * 6 , respectively, with 'fH' indicating the fold-Hopf bifurcation. The first digit of the index refers to S * 5 , while the second corresponds to S * 6 , and they indicate that in that parameter region the respective steady state is unfeasible (index is '0'), stable (index is '1'), unstable via Hopf bifurcation with a periodic solution around this steady state (index is '2'), or unstable via a steady-state bifurcation (index is '3'). In all plots, the condition β < d F holds, so the disease-free steady state S * 2 is also stable.
Since T * reg satisfies conditions (3) and (6), the steady state S * 5 /S * 6 is stable for τ 2 = 0. Lemma 3.1 then ensures that τ * is the first positive value of the time delay τ 2 , for which the roots of the characteristic Equation (4) cross the imaginary axis with positive speed. Hence, the steady state S * 5 /S * 6 is stable for 0 ≤ τ 2 < τ * , unstable for τ 2 > τ * , and undergoes a Hopf bifurcation at τ 2 = τ * . Remark. A similar result can be formulated for a subcritical Hopf bifurcation of the steady state S * 5 /S * 6 at some higher value of τ 2 .
The only remaining steady state is the persistent (chronic) equilibrium S * 8 with all of its components being positive. Since it did not prove possible to find a closed form expression for this steady state, its stability also has to be studied numerically.

Numerical Stability Analysis and Simulations
To investigate the role of different parameters in the dynamics of model (1), in this section we perform a detailed numerical bifurcation analysis and simulations of this model. Stability of different steady states is determined numerically by computing the largest real part of the characteristic eigenvalues, which is achieved by using a pseudospectral method implemented in a traceDDE suite in MATLAB [61].
Analytical results from the previous section suggest that at β = d F , the disease-free steady state S * 2 undergoes a transcritical bifurcation. For β < d F , the disease-free steady state S * 2 is stable, and the chronic steady state is infeasible. On the contrary, when β > d F , the disease-free steady state S * 2 is unstable, and in this case it makes sense to investigate stability of the chronic steady state. Therefore, these two cases are considered separately, and as a first step we fix the baseline values as given in Table 1. For this choice of parameters, we have d F − β > 0, implying that S * 2 is always stable, and Figure 2 illustrates how the stability of S * 5 and S * 6 is affected by parameters. This figure indicates that the steady states S * 5 and S * 6 are only biologically feasible if the regulatory T cells do not grow too rapidly and do not clear autoreactive T cells too quickly. Importantly, Figure 2 shows that the value of the rate δ 2 of clearance of IL-2 by regulatory T cells does not have any effect on the thresholds of λ r and δ 1 , where the steady states S * 5 and S * 6 lose their feasibility. Moreover, if λ r and δ 1 are small, then increasing the rate δ 2 at which Tregs inhibit the production of IL-2 makes S * 6 become unfeasible, resulting in a stable steady state S * 5 , which has the zero population of host cells A. On the other hand, the steady state S * 6 associated with autoimmune responses is favoured for higher values of δ 1 and λ r . In the case stable periodic solutions around these steady states, increasing δ 2 results in the disappearance of oscillations and stabilisation of the associated steady state. At the intersection of the lines of Hopf bifurcation and the steady-state bifurcation, as determined by Theorem 3.2 and conditions (5), one has the co-dimension two fold-Hopf (also known as zero-Hopf or saddle-node Hopf) bifurcation [62].
Since our earlier analysis showed that stability of the steady states S * 5 /S * 6 is affected by the time delay τ 2 , in Figure 3 we consider stability of these equilibria depending on τ 2 and the rate δ 2 . For the steady state S * 5 , if the effect of IL-2 on promoting proliferation of T cells is fast (i.e., τ 2 is small), there is a large range of δ 2 , starting with some very low values, for which S * 5 is stable. Increasing  Table 1. White area shows the region where the steady state S * 6 is infeasible. Colour code denotes max[Re(λ)] for the steady states when they are feasible. In all plots the condition d F > β holds, so the disease-free steady state S * 2 is stable. Basins of attraction of different steady states depending on the initial conditions (c), with other initial conditions specified in (8), and parameter values from Table 1, except for τ 2 = 18. Cyan and pink areas are the basins of attraction of S * 2 and S * 6 , respectively.
the time delay τ 2 results in the Hopf bifurcation of this steady state as described in Theorem 3.2.
One should note that for intermediate values of δ 2 , the steady state S * 5 undergoes stability switches, whereby increasing the delay τ 2 further results in a subcritical Hopf bifurcation, which stabilises S * 5 , but after some number of such stability switches eventually the steady state S * 5 is unstable. For higher still values of δ 2 , the steady state S * 5 remains stable for an entire range of τ 2 values, and the only way to lose its stability is via a steady state bifurcation as given by (5). In the case of autoimmune steady state S * 6 , the situation is somewhat different in that increasing δ 2 beyond some critical values makes this steady state biologically infeasible. At the same time, for an entire range of δ 2 values where it is feasible, this steady state exhibits a single loss of stability through a Hopf bifurcation for some critical value of the time delay τ 2 , in agreement with Theorem 3.2.
As mentioned earlier, for parameter values used in Figure 3, the disease-free steady state S * 2 is stable. Hence, the system exhibits a bi-stability between a disease-free state and either stable steady states S * 5 /S * 6 , or periodic solutions around these steady states. To investigate this bi-stability, we choose parameter values as in Table 1 except for τ 2 = 18, which corresponds to a stable steady state S * 6 , and we fix initial conditions for state variables as follows, except for initial amounts of infected cells and regulatory T cells that are allowed to vary. Figure 3c illustrates the bi-stability between S * 2 and S * 6 in terms of their basins of attraction. It is worth noting that recently significant research in approximation theory and meshless interpolation has focused on developing techniques for detection and analysis of attraction basins [63,64,65,66,67,68]. Figure 3c suggests that for very large initial amounts of regulatory T cells, the system converges to the disease-free steady state. It also indicates that if the initial amount of infected cells is very small or is bigger than some specific value, then the infection will be cleared. Interestingly, increasing the initial amount of the regulatory T cells results in a larger range of initial amounts of infection, for which the system tends to a stable autoimmune state S * 6 . In Figure 3b we discovered that increasing τ 2 makes the autoimmune steady state S * 6 undergo a Hopf bifurcation, in which case the system will exhibit a bi-stability between stable S * 2 and a periodic solution around S * 6 . Our numerical investigation suggests that the shape of basins of attraction in this case is qualitatively similar to that shown in Figure 3c, with the basin of attraction of the stable steady state S * 6 being replaced by the basin of attraction of the periodic solution around this steady state. Figure 4 shows temporary evolution of the system (1) in the regime of bi-stability between a stable disease-free steady state and a stable autoimmune steady state S * 6 (similar pattern of be- haviour is exhibited in the case of bi-stability between S * 2 and S * 5 ). It also illustrates how the system develops a periodic solution around the steady state S * 6 for a higher value of τ 2 . Periodic oscillations around the steady state S * 6 biologically correspond to a genuine autoimmune state: after the initial infection is cleared, the system exhibits sustained endogenous oscillations, characterised by periods of significant reduction in the number of organ cells through a negative action of autoreactive T cells, separated by periods of quiescence. This type of behaviour is often observed in clinical manifestations of autoimmune disease [44,45,46,47]. This result has substantial biological significance as effectively it suggests that even for the same kinetic parameters of immune response, the ultimate state of the system, which can be either a successful clearance of infection without lasting consequences, or progression to autoimmunity, also depends on the strength of the initial infection and of the initial state of the immune system, as represented by the initial number of regulatory T cells.
Next we consider a situation where β > d F , so the disease-free steady state is unstable, and the system can have three steady states S * 5 , S * 6 and S * 8 . Our earlier results [48] suggest that in the case where regulatory T cells do not inhibit the production of IL-2, i.e., for δ 2 = 0, the steady state S * 6 is stable. Figure 5 shows regions of feasibility and stability of these steady states depending on δ 2 and τ 2 in this case. One observes that S * 5 and S * 6 , whose stability boundaries are determined by Theorem 3.2, exhibit the same behaviour as in Figure 3, namely, for S * 5 increasing τ 2 causes multiple stability switches for smaller values of δ 2 , and the steady state is unstable for very small δ 2 and stable for large δ 2 ; in contrast, S * 5 exhibits a single loss of stability via Hopf bifurcation at  Table 1, except for β = 1.4 and σ 2 = 1, so that β > d F . White area shows the region where the steady state is infeasible. Colour code denotes max[Re(λ)] for each steady states when it is feasible. (d) Summary of stability results. Green indicates the region where S * 6 and S * 8 are stable, and S * 5 is unstable, whereas red is the area where S * 5 and S * 8 are stable, and S * 6 is infeasible. Yellow is where S * 8 is stable, S * 5 is unstable, and S * 6 is infeasible. Purple shows the region where S * 6 is stable, but S * 5 and S * 8 are unstable. Blue and cyan indicate the regions where S * 5 and S * 6 are unstable, but S * 8 is stable or unstable, respectively.  Figure 5d divides the δ 2 -τ 2 plane into different regions based on feasibility and stability of these steady states and shows that increasing δ 2 makes the autoimmune steady state S * 6 infeasible. In other regions, the system can exhibit a bi-stability between a stable steady state S * 8 and either a stable steady state S * 5 , or a periodic solution around S * 5 . Figure 6 illustrates the basins of attraction of the steady states S * 5 , S * 6 and S * 8 , as well as periodic solutions around S * 8 . Figure 6a shows the basins of attraction of the steady states S * 5 and S * 8 and demonstrates that if the initial number of regulatory T cells or infected cells is sufficiently high, or the initial amount of infected cells is very low, the immune response neither eliminates infection nor clears autoreactive T cells, and the system approaches the stable steady state S * 5 . Figure 6b illustrates bi-stability between the stable steady state S * 6 and a periodic solution around S * 8 , and has a different behaviour to than shown in Figure 6a. This figure suggests that for a specific range of F (0) the system converges to a stable autoimmune state S * 6 for all values of T reg (0). However, if the initial number of infected cells is very high or very low, the system instead develops a periodic solution around the steady state S * 8 associated with chronic infection. Figure 7 illustrates a regime of bi-stability between a stable steady state S * 6 and a periodic solution around S * 8 for combinations of initial conditions indicated by crossed in Figure 6b. It also illustrates how the system develops a stable solution around the steady state S * 8 for a higher value of τ 2 . This figure shows that by increasing the initial number of infected cells the behaviour of the system changes, as it then approaches the autoimmune steady state S * 6 . Interestingly, one can observe that for high values of F (0) the system can eliminate the infection, but it cannot clear the autoreactive T cells, in which case the system converges to S * 6 . On the other hand, for a smaller number of infected cells the system develops a periodic solution around the endemic steady state. Figure 8 shows how the stability of the chronic infection steady state S * 8 changes with respect to time delays. Figure 8a indicates that for small values of τ 2 (i.e., when the influence of IL-2 on proliferation of T cells is occurring quite rapidly), the steady state S * 8 is stable, and increasing the time delay τ 1 associated with viral eclipse phase does not have an effect on its stability. At the same time, if τ 2 exceeds some specific value, by increasing τ 1 the chronic steady state switches between being stable or unstable. Figure 8b demonstrates a different behaviour, suggesting that for each value of τ 1 , there is small range of τ 3 values where S * 8 is stable, but for smaller and larger values of τ 3 it is unstable. For intermediate values of the eclipse phase delay τ 1 , there is an additional narrow range of τ 3 values where S * 8 is stable. Figure 8c illustrates that for very small, respectively very large, values of τ 3 , the chronic infection steady state is stable, respectively unstable for any value of τ 2 ; for intermediate values of τ 3 , this steady state undergoes a finite number of stability switches for increasing values of τ 2 and eventually becomes unstable.
It should be noted Figure 8 shows that unlike τ 1 and τ 2 , once the steady state S * 8 loses stability via Hopf bifurcation due to increasing τ 3 , it cannot regain stability for higher values of τ 3 .

Conclusions
In this paper, we have developed and analysed a time-delayed model of immune response to a viral infection, which accounts for T cells with different activation thresholds, a cytokine mediating T cell proliferation, as well as regulatory T cells. Particular attention is payed to the dual suppressive role of regulatory T cells in terms of reducing the amount of autoreactive T cells, and also inhibiting IL-2. To achieve better biological realism of the model, we have explicitly included time delays associated with the eclipse phase of the virus life cycle, stimulation/proliferation of T cells by IL-2, and suppression of IL-2 by regulatory T cells. Depending on the values of parameters, the system can have four steady states: the disease-free state, the state characterised by the death of host cells, the autoimmune state, and a state of chronic infection. We have established conditions for stability and steady-state or Hopf bifurcations of these steady states in terms of system parameters.
In the case where the natural death rate of infected cells exceeds the infection rate, the immune system is able to clear the infection, and the disease-fee steady state is stable. In this regime, the system can also support the autoimmune steady state or the steady state with the death of host cells, either of which can be stable, or give rise to a periodic solution emerging via Hopf bifurcation. In the opposite case, when the natural death rate of infected cells is smaller than the infection rate, the disease-fee steady state is unstable, but it is possible to have a bi-stability between the other three steady states or periodic solutions around them. To better understand bi-stability between different dynamical regimes, we have used numerical simulations to identify basins of attraction of different steady states and periodic solutions depending on the initial level of infection and the initial number of regulatory T cells. The fact that for the same parameter values the system can exhibit bi-stability between a disease-free steady state and an autoimmune state, represented by sustained periodic oscillations following the clearance of infection, is very important from a clinical point of view, as effectively it suggests that the progress and eventual outcome of viral infection is also determined by the strength of infection and the initial state of the immune system. One counter-intuitive observation is that in the case of bi-stability with a disease-free steady state, for higher initial numbers of regulatory T cells, the autoimmune steady state is actually stable for a wider range of initial levels of infection. In this regime of bi-stability, increasing the time delay associated with the positive impact of IL-2 on proliferation of T cells results in the loss of stability of autoimmune steady state and emergence of autoimmune dynamics, characterised by stable periodic oscillations. On the contrary, in the case where the disease-free steady state is unstable, increasing this time delay results in stabilisation of the chronic infection.
There are several directions in which the work presented in this paper can be extended. One direction is exploration of the contributions from other components of immune response, more specifically, antibodies and memory T cells, to the onset and progress of autoimmunity [69,70]. This is particularly important from the perspective that clinically the onset of autoimmune disease is often taking place on a much longer scale than the timescale of a regular immune response to a viral infection, so memory T cells can be expected to play a more substantial role. While our model has focused on one specific growth cytokine IL-2, a number of other cytokines, such as IL-7 [71], TNFβ and IL-10 [29], are known to significantly affect homeostasis and proliferation of different types of T cells, as well as mediate their efficiency in eliminating the infection. Including these immune mediators explicitly in the model can provide further significant insights into the dynamics of immune response, as has been recently demonstrated on the example of a detailed model of immune response to hepatitis B [72]. Another biologically relevant and mathematically challenging problem is the investigation of the interplay between stochasticity, which is known to be an intrinsic feature of immune response [73,49], and effects of time delays associated with various aspects of immune dynamics.

Appendix
Coefficients of Equation (7) for Hopf frequency are given below.