Effects of viral and cytokine delays on dynamics of autoimmunity
EEffects of Viral and Cytokine Delays on Dynamics of Autoimmunity
F. Fatehi Chenar, Y.N. Kyrychko, K.B. Blyuss ∗ Department of Mathematics, University of Sussex, Falmer,Brighton, BN1 9QH, United KingdomMay 21, 2018
Abstract
A major contribution to the onset and development of autoimmune disease is known to comefrom infections. An important practical problem is identifying the precise mechanism by which thebreakdown 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, whichincludes T cells with different activation thresholds, regulatory T cells (Tregs), and a cytokinemediating immune dynamics. Particular emphasis is made on the role of time delays associatedwith the processes of infection and mounting the immune response. Stability analysis of varioussteady states of the model allows us to identify parameter regions associated with different typesof immune behaviour, such as, normal clearance of infection, chronic infection, and autoimmunedynamics. Numerical simulations are used to illustrate different dynamical regimes, and to identifybasins of attraction of different dynamical states. An important result of the analysis is that notonly the parameters of the system, but also the initial level of infection and the initial state of theimmune system determine the progress and outcome of the dynamics.
An immune system can only be viewed as effective when it can robustly identify and destroypathogen-infected cells, while also distinguishing such cells from healthy cells. In the case ofbreakdown of immune tolerance, the immune system fails to discriminate between self-antigensand foreign antigens, which results in autoimmune disease, i.e., undesired destruction of healthycells. 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 betweenepitopes associated with foreign and self-antigens can lead to a T cell response against healthy hostcells [1, 2].For many autoimmune diseases, the disease occurs in a specific organ or part of the body, suchas retina in uveitis, central nervous system in multiple sclerosis, or pancreatic β -cells in type-1 dia-betes [3, 4, 5]. It is extremely difficult to identify the specific causes of autoimmunity in individualpatients, as it usually has contributions from a number of internal and external factors, includinga 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 thatsome additional environmental triggers are required for the onset of autoimmunity, and these areusually represented by infections [10, 11]. A very recent work has experimentally identified a gutbacterium that, when present in mice and humans, can migrate to other parts of the body, facilitat-ing subsequent triggering of autoimmune disease in those organs [12]. Various mechanisms of onset ∗ Corresponding author. Email: [email protected] a r X i v : . [ q - b i o . Q M ] M a y f pathogen-induced autoimmune disease have been identified, including bystander activation [13]and molecular mimicry [14, 15], which is particularly important in the context of autoimmunitycaused by viral infections.A number of mathematical models have looked into dynamics of onset and development of au-toimmune disease. Segel et al. [16] analysed interactions between effector and regulatory T cellsin 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 onsetof autoimmune state, defined as stable above-threshold oscillations in the number of autoreactivecells, as a result of interactions between regulatory and autoreactive T cells. Le´on et al. [19, 20, 21]and Carneiro et al. [22] have studied interactions between different T cells, with an emphasis onthe suppressing role of regulatory T cells in the dynamics of immune response and control ofautoimmunity. Alexander and Wahl [23] have also looked into the role of regulatory T cells, inparticular focusing on their interactions with professional APCs and effector cells for the purposeof controlling immune response. Iwami et al. [24, 25] explicitly included a separate compartmentrepresenting the viral population in their models of immune response, and showed that the func-tional form of the growth function for susceptible host cells can have a significant effect on theresulting immune dynamics. Despite being able to explain the emergence of autoimmunity as a by-product of immune response to infection, these models were not able to exhibit another practicallyimportant dynamical regime of normal viral clearance. For the case of pathogen-induced autoim-munity arising through bystander activation, Burroughs et al. [26, 27, 28] have developed a modelthat investigates interactions between T cells and interleukin-2 (IL-2), an important cytokine, inmediating the onset of autoimmunity.Among various parts of the immune system involved in coordinating an effective immune re-sponse, a particularly significant role is known to be played by the T cells, with experimentalevidence 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] andBurroughs et al. [26, 27] have explicitly included a separate compartment for regulatory T cells thatare activated by autoantigens and suppress the activity of autoreactive T cells. Another frameworkfor modelling the effects of T cells on autoimmune dynamics is by using the idea that T cells havedifferent or tunable activation thresholds (TAT), which result in different immune functionality ofthe same T cells, and also allow T cells to adjust their response to simulation by autoantigens. Thisapproach was proposed for the analysis of the peripheral and central T cell dynamics [33, 34, 35], ithas also been used to study differences in activation/response thresholds that are dependent on theactivation state of the T cell [36]. Murine and human experiments have confirmed that activationof T cells can indeed change dynamically during their circulation [37, 38, 39, 40]. To model thisfeature, van den Berg and Rand [41], and Scherer et al. [42] developed stochastic models for thetuning of activation thresholds.Blyuss and Nicholson [43, 44] have proposed a mathematical model of autoimmunity resultingfrom immune response to a viral infection through a mechanism of molecular mimicry. This modelexplicitly 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 responsecan 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 stableoscillations in the amounts of cell populations. From a clinical perspective, such behaviour is tobe expected, as it is associated with relapses and remissions that have been observed in a numberof autoimmune diseases, such as autoimmune thyroid disease, MS, and uveitis [45, 46, 47]. Onedeficiency of this model is the fact that the oscillations associated with autoimmune regime can onlyoccur 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 hasbeen fully cleared. To overcome this limitation, Fatehi et al. [48] have recently developed a moreadvanced model that also includes regulatory T cells and cytokines, which has allowed the authorsto obtain a more realistic representation of immune response and various dynamical regimes. Aparticularly important practical insight provided by this model is the observation that it is not onlythe system parameters, but also the initial level of infection and the initial state of the immune ystem, that determine whether the host will just successfully clear the infection, or will proceedto 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 anddetermining which of the immune states is more likely to be attained. The authors have alsodetermined an experimentally important characterisation of autoimmune state, as provided by thedependence of variance in cell populations on various system parameters.In this paper, we develop and analyse a model of autoimmune dynamics, with particular focus onthe role of time delays associated with different aspects of immune response, as well as an inhibitingeffect of regulatory T cells on secretion of IL-2. In the next section, we introduce the model anddiscuss its basic properties. Section 3 contains systematic analysis of all steady states, includingconditions for their feasibility and stability. Section 4 contains a bifurcation analysis of the modeland demonstrates various types of behaviour that the system exhibits depending on parametersand initial conditions, which includes identification of attraction basins of various states. The paperconcludes in Section 5 with the discussion of results. To understand how interactions between different parts of the immune system and the cytokinecan lead to autoimmunity, we consider a model illustrated in a diagram shown in Figure 1. Inthis model, unlike earlier work of Blyuss and Nicholson [43, 44], we consider a situation where aviral infection and possible autoimmunity occur in the same organ of the host. The healthy hostcells, whose number is denoted by A ( t ), in the absence of infection are assumed to grow logisticallywith linear growth rate r and carrying capacity N , and they acquire infection at rate β from theinfected cells F ( t ). Since experimental evidence suggests that antibodies play a secondary rolecompared to T cells [50], and autoimmunity can develop even in the absence of B cells [51], we donot 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 areactivated through interaction with infected cells, which occurs at rate α , a proportion p of them willgo on to differentiate into additional regulatory T cells, a fraction p will become normal activatedT cells T nor ( t ) able to destroy infected cells at rate µ F upon recognition of foreign antigen presenton their surface. The remaining proportion of (1 − p − p ) of T cells will become autoreactiveT cells T aut ( t ) that, in light of their lower activation threshold will be eliminating both infectedcells and healthy host cells at rate µ a due to the above-mentioned cross-reactivity between someepitopes in self- and foreign antigens. Regulatory T cells T reg ( t ) are assumed to be in their ownhomeostasis [52], and their main contribution to immune dynamics lies in suppressing autoreactiveT cells at rate δ . To reduce the dimensionality of the model, it is assumed that the process ofviral 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 istaken to be proportional to the number of infected cells, and this eliminates the need for a separatecompartment for free virus.A number of different cytokines mediate immune response to infection, and in the context ofT cell dynamics, a particular important role is played by interleukin 2 (IL-2), represented by thevariable I ( t ) in the model, which acts to enhance the proliferation of T cells, which, in turn, secretefurther 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]. Torepresent this mathematically, we will assume that T nor and T aut produce IL-2 at rates σ and σ , and conversely, IL-2 enhances proliferation of T reg , T nor and T aut at rates ρ , ρ , and ρ . Weinclude in the model suppression of IL-2 by regulatory T cells at rate δ , in a manner similar toBurroughs et al. [28].While the production of new virus particles by infected cells is assumed to be fast, we explicitlyinclude in the model time delay τ associated with the actual process of infection, which includesmultiple stages of the eclipse phase of viral life cycle, such as virus attachment, cell penetration anduncoating [56, 57]. We also include the time delay τ associated with simulation and proliferation A schematic diagram of immune response to infection. Blue indicates host cells (susceptible andinfected), red denotes different T cells (na¨ıve, regulatory, normal activated, and autoreactive T cells), yellowshows cytokines (interleukin 2). τ i inside each of the subnetworks shows the time delay associated with thatprocess. f T cells by IL-2, and the time delay τ between antigen encounter and resulting T cell expansion[58].With the above assumptions, the complete model takes the form dAdt = rA (cid:18) − AN (cid:19) − βAF − µ a T aut A,dFdt = βA ( t − τ ) F ( t − τ ) − 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 ( t − τ ) F ( t − τ ) + ρ I ( t − τ ) T reg ( t − τ ) ,dT nor dt = p αT in ( t − τ ) F ( t − τ ) − d n T nor + ρ I ( t − τ ) T nor ( t − τ ) ,dT aut dt = (1 − p − p ) αT in ( t − τ ) F ( t − τ ) − d a T aut − δ T reg T aut + ρ I ( t − τ ) T aut ( t − τ ) ,dIdt = σ T nor + σ T aut − δ T reg I − d i I. Introducing non-dimensional variablesˆ t = rt, A = N ˆ A, F = N ˆ F , T in = λ in d in ˆ T in , T reg = λ in d in ˆ T reg ,T nor = λ in d in ˆ T nor , T aut = λ in d in ˆ T aut , I = λ in d in ˆ I, where ˆ β = βNr , ˆ µ a = µ a λ in rd in , ˆ d F = d F r , ˆ µ F = µ F λ in rd in , ˆ d in = d in r , ˆ α = αNr , ˆ λ r = λ r d in λ in r , ˆ d r = d r r , ˆ d n = d n r , ˆ d a = d a r , ˆ ρ i = ρ i λ in rd in , i = 1 , , , ˆ δ = δ λ in rd in , ˆ δ = δ λ in rd in , ˆ σ = σ r , ˆ σ = σ r , ˆ d i = d i r , yields a rescaled model dAdt = A (1 − A ) − βAF − µ a T aut A,dFdT = βA ( T − τ ) F ( T − τ ) − d F F − µ F T nor F − µ a T aut F,dT in dT = d in (1 − T in ) − αT in F,dT reg dT = λ r − d r T reg + p αT in ( T − τ ) F ( T − τ ) + ρ I ( T − τ ) T reg ( T − τ ) ,dT nor dT = p αT in ( T − τ ) F ( T − τ ) − d n T nor + ρ I ( T − τ ) T nor ( T − τ ) ,dT aut dT = (1 − p − p ) αT in ( T − τ ) F ( T − τ ) − d a T aut − δ T reg T aut + ρ I ( T − τ ) T aut ( T − τ ) ,dIdT = σ T nor + σ T aut − δ T reg I − d i I, (1)where all hats in variables and parameters have been dropped for simplicity of notation, and allparameters 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 ≥ S ∗ = (cid:0) A ∗ , F ∗ , T ∗ in , T ∗ reg , T ∗ nor , T ∗ aut , I ∗ (cid:1) , hat can be found by equating to zero the right-hand sides of Equation (1) and solving the resultingsystem of algebraic equations, deferring the discussion of conditionally stable steady states toSection 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 ofsteady 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 S ∗ = (cid:18) , , , λ r d r , , , (cid:19) , S ∗ = (cid:18) , , , λ r d r , , , (cid:19) , of which S ∗ is always unstable, and S ∗ is a disease-free conditionally stable steady state.For T ∗ nor (cid:54) = 0 and T ∗ aut = 0, we again have two steady states S ∗ = (cid:18) , , , λ r ρ ρ d r − ρ d n , T ∗ nor , , d n ρ (cid:19) , S ∗ = (cid:18) , , , λ r ρ ρ d r − ρ d n , T ∗ nor , , d n ρ (cid:19) , where T ∗ nor = d n ( λ r δ ρ + d i d r ρ − d i d n ρ ) ρ σ ( ρ d r − ρ d n ) , but they are both unstable for any values of pa-rameters. In the case when T ∗ nor = 0 and T ∗ aut (cid:54) = 0, we have two further steady states S ∗ and S ∗ , S ∗ = (cid:32) , , , T ∗ reg , , (cid:0) d i + δ T ∗ reg (cid:1) (cid:0) d a + δ T ∗ reg (cid:1) ρ σ , d a + δ T ∗ reg ρ (cid:33) ,S ∗ = (cid:32) A ∗ , , , T ∗ reg , , (cid:0) d i + δ T ∗ reg (cid:1) (cid:0) d a + δ T ∗ reg (cid:1) ρ σ , d a + δ T ∗ reg ρ (cid:33) , where A ∗ = 1 − µ a (cid:0) d i + δ T ∗ reg (cid:1) (cid:0) d a + δ T ∗ reg (cid:1) ρ σ , and T ∗ reg = d r ρ − ρ d a ± (cid:113) ( d r ρ − ρ d a ) − ρ δ λ r ρ ρ δ . The steady state S ∗ has A ∗ = 0, which implies the death of host cells, whereas the steady state S ∗ corresponds to an autoimmune regime. The steady state S ∗ with T ∗ nor (cid:54) = 0 and T ∗ aut (cid:54) = 0 existsonly for a particular combination of parameters, namely, when δ ρ λ r = ( ρ d n − ρ d a )( ρ d r − ρ d n ) , and is always unstable. Finally, when F ∗ (cid:54) = 0, the system (1) can have a steady state S ∗ with allof its components being positive, but it does not appear possible to find a closed form expressionfor this state.In summary, besides the unconditionally unstable steady states, the model (1) has at most forconditionally stable steady states: the disease-free steady state S ∗ , the steady state with the deathof host cells S ∗ , the autoimmune steady state S ∗ , and the persistent or chronic steady state S ∗ . Linearising the system (1) near the disease-free steady state S ∗ yields the following equation forcharacteristic roots λ λ + d F − βe − λτ = 0 . (2)If d F < β , the above equation always has a real positive root for any value τ ≥
0, implying thatthe disease-free steady state is always unstable for any value of the time delays. If, however, the ondition d F > β holds, the disease-free steady state is stable for τ = 0. To find out whether itcan lose stability for τ >
0, we look for solutions of Equation (2) in the form λ = iω . Separatingreal and imaginary parts yields d F = β cos( ωτ ) ,ω = − β sin( ωτ ) . Squaring and adding these two equations gives the following equation for potential Hopf fre-quency ω ω + d F − β = 0 . since d F > β , this equation does not have real roots for ω , suggesting that there can be no rootsof the form λ = iω of the characteristic Equation (2). This implies that in the case d F > β thedisease-free steady state S ∗ is stable for all values of the time delay τ ≥ The steady state S ∗ (respectively, S ∗ ) is stable if P < d a + δ T ∗ reg ρ < d n ρ , (3)and all roots of the following equation have negative real part∆( τ , λ ) = p ( λ ) e − λτ + p ( λ ) e − λτ + p ( λ ) = 0 , (4)where p ( λ ) = ρ (cid:0) d a + δ T ∗ reg (cid:1) ρ (cid:0) λ + 2 d i + δ T ∗ reg (cid:1) ,p ( λ ) = − (cid:0) d a + δ T ∗ reg (cid:1) ρ (cid:110) ( ρ + ρ ) λ + (cid:2) ρ (cid:0) d i + d a + δ T ∗ reg (cid:1) + ρ (cid:0) d r + 2 d i + 2 δ T ∗ reg (cid:1)(cid:3) λ + d i ( ρ d a + 2 d r ρ ) + δ T ∗ reg (cid:0) − ρ δ T ∗ reg + 2 d r ρ (cid:1) (cid:111) ,p ( λ ) = ( λ + d r ) (cid:0) λ + d i + δ T ∗ reg (cid:1) (cid:0) λ + d a + δ T ∗ reg (cid:1) , and P = σ µ a (cid:0) d i + δ T ∗ reg (cid:1) , for S ∗ ,σ ( β − d F ) µ a (1 + β ) (cid:0) d i + δ T ∗ reg (cid:1) , for S ∗ . This steady state undergoes a steady-state bifurcation if d a + δ T ∗ reg ρ = P, or d a + δ T ∗ reg ρ = d n ρ , or δ ρ (cid:0) T ∗ reg (cid:1) = λ r ρ . (5)For τ = 0 these steady states are stable if T ∗ reg satisfies (3) and δ ρ (cid:0) T ∗ reg (cid:1) > λ r ρ ,a (cid:0) T ∗ reg (cid:1) + a (cid:0) T ∗ reg (cid:1) + a (cid:0) T ∗ reg (cid:1) + a (cid:0) T ∗ reg (cid:1) + a T ∗ reg + a > , (6)where a = − δ δ ( δ ρ − δ ρ + δ ρ ) , a = d a δ ( δ ρ − δ ρ − δ ρ ) − d i δ ( δ ρ − δ ρ + 2 δ ρ ) ,a = − d i δ ( d a ρ + d i ρ ) + d a d i δ ( ρ − ρ ) + λ r δ ( δ ρ + δ ρ ) ,a = − d a d i ρ + λ r δ ( d a ρ + 2 d i ρ ) , a = λ r ρ ( d i + δ λ r ) , a = d i ρ λ r . o investigate whether stability can be lost for τ >
0, we use an iterative procedure described in[59, 60] to determine a function F ( ω ), whose roots give the Hopf frequency associated with purelyimaginary roots of Equation (4). Substituting λ = iω into Equation (4), we define ∆ (1) ( τ , λ ) as∆ (1) ( τ , λ ) = p ( iω )∆( τ , iω ) − p ( iω ) e − iωτ ∆( τ , iω ) = p (1)0 ( iω ) + p (1)1 ( iω ) e − iωτ , where p (1)0 ( iω ) = | p ( iω ) | − | p ( iω ) | ,p (1)1 ( iω ) = p ( iω ) p ( iω ) − p ( iω ) p ( iω ) , and the bar denotes the complex conjugate. If we define F ( ω ) = (cid:12)(cid:12)(cid:12) p (1)0 ( iω ) (cid:12)(cid:12)(cid:12) − (cid:12)(cid:12)(cid:12) p (1)1 ( iω ) (cid:12)(cid:12)(cid:12) , then ∆( τ , iω ) = 0 whenever ω is a root of F ( ω ) = 0. The function F ( ω ) has the explicit form F ( ω ) = ω + b ω + b ω + b ω + b ω + b ω + b , with b = (cid:0) δ d a + T ∗ reg (cid:1) ρ (cid:0) d i + δ T ∗ reg (cid:1) (cid:0) T ∗ reg δ ρ + d a ρ − d r ρ (cid:1)(cid:104) (cid:0) d i + δ T ∗ reg (cid:1) ( d a ρ + 3 d r ρ ) + 2 d i ρ (cid:0) d a + δ T ∗ reg (cid:1) (cid:3)(cid:104) ρ (cid:0) d a + δ T ∗ reg (cid:1) (cid:0) d i + δ T ∗ reg (cid:1) − d r ρ (cid:0) d i + δ T ∗ reg (cid:1) (cid:105) . The explicit formulae for other coefficients of F ( ω ) can be found in Appendix 5. Introducing s = ω , the equation F ( ω ) = 0 can be equivalently rewritten as follows, h ( s ) = s + b s + b s + b s + b s + b s + b = 0 . (7)Without loss of generality, suppose that Equation (7) has six distinct positive roots denoted by s , s , ... , s , which means that the equation F ( ω ) = 0 has six positive roots ω i = √ s i , i = 1 , , ..., . Substituting λ k = iω k into Equation (4) gives τ k,j = 1 ω k (cid:34) arctan (cid:32) ω k (cid:0) ( ρ + ρ ) ω k + f ω k + f (cid:1)(cid:0) ρ Z − d r ρ − ρ I ∗ − ρ δ T ∗ reg (cid:1) ω k + g ω k + g (cid:33) + jπ (cid:35) , for k = 1 , , ..., j = 0 , , , ... , where f = − ρ ρ I ∗ Z − ρ ρ (2 ρ + 3 ρ ) I ∗ Z + ρ ρ T ∗ reg ( − δ ρ + 3 δ ρ + δ ρ ) I ∗ Z − T ∗ reg δ ρ ρ I ∗ − T ∗ reg δ ρ ρ I ∗ Z + d r ρ (cid:0) − δ ρ T ∗ reg + d r ρ (cid:1) I ∗ Z + d r (cid:0) − T ∗ reg δ ρ + 2 d r ρ (cid:1) Z ,f = − ρ ρ I ∗ + ρ I ∗ Z + ( ρ + 2 ρ ) Z + ρ T ∗ reg ( δ − δ ) Z − d r (cid:0) T ∗ reg δ ρ − d r ρ (cid:1) ,g = ρ ρ I ∗ (cid:16) Z − T ∗ reg δ Z + T ∗ reg δ (cid:17) + ρ ρ (cid:0) − T ∗ reg δ ρ + 3 d r ρ (cid:1) I ∗ Z + ρ ρ δ T ∗ reg (cid:0) T ∗ reg δ ρ − d r ρ (cid:1) I ∗ Z + d r ρ (cid:0) T ∗ reg δ ρ − d r ρ (cid:1) I ∗ Z ,g = ρ I ∗ (cid:16) ρ ρ I ∗ − ρ I ∗ Z − ρ Z − T ∗ reg δ ρ Z − d r ρ (cid:17) − ρ (cid:0) d r + δ T ∗ reg (cid:1) Z + d r (cid:0) − T ∗ reg δ ρ + T ∗ reg δ ρ + d r ρ (cid:1) Z, nd I ∗ = d a + δ T ∗ reg ρ , Z = d i + δ T ∗ reg . This allows us to find τ ∗ = τ k , = min ≤ k ≤ { τ k, } , ω = ω k , as the first time delay for which the roots of the characteristic Equation (4) cross the imaginaryaxis. To determine whether these steady states actually undergo a Hopf bifurcation at τ = τ ∗ , wehave to compute the sign of d Re[ λ ( τ ∗ )] /dτ . For τ = τ ∗ , λ ( τ ∗ ) = iω , and we also define s = ω . Lemma 3.1.
Suppose h (cid:48) ( s ) (cid:54) = 0 and p (1)0 ( iω ) (cid:54) = 0 . Then the following transversality conditionholds sgn (cid:40) d Re( λ ) d τ (cid:12)(cid:12)(cid:12)(cid:12) τ = τ ∗ (cid:41) = sgn[ p (1)0 ( iω ) h (cid:48) ( s )] . Proof.
Considering p j ( iω ) = x j ( ω ) + iy j ( ω ) for j = 0 , ,
2, we have p (1)0 ( iω ) = x + y − x − y ,p (1)1 ( iω ) =( x x + y y − x x − y y ) + ( x y + x y − x y − x y ) i, where all x j and y j are expressed in terms of system parameters and steady state values of thevariables. Substituting these expressions into ∆( τ , iω ) = 0 and ∆ (1) ( τ , iω ) = 0, and thenseparating real and imaginary parts gives x cos(2 ω τ ∗ ) + y sin(2 ω τ ∗ ) + x cos( ω τ ∗ ) + y sin( ω τ ∗ ) = − x ,y cos(2 ω τ ∗ ) − x sin(2 ω τ ∗ ) + y cos( ω τ ∗ ) − x sin( ω τ ∗ ) = − y , ( x x + y y − x x − y y ) cos( ω τ ∗ ) + ( x y + x y − x y − x y ) sin( ω τ ∗ )= − x − y + x + y , ( x y + x y − x y − x y ) cos( ω τ ∗ ) − ( x x + y y − x x − y y ) sin( ω τ ∗ ) = 0 . Solving this system of equations provides the values of sin( ω τ ∗ ), cos( ω τ ∗ ), sin(2 ω τ ∗ ), and cos(2 ω τ ∗ ).Taking the derivative of Equation (4) with respect to τ , one finds (cid:18) d λd τ (cid:19) − = p (cid:48) ( λ ) e − λτ + p (cid:48) ( λ ) e − λτ + p (cid:48) ( λ ) λ (2 p ( λ ) e − λτ + p ( λ ) e − λτ ) − τ λ . Hence, (cid:32) d Re( λ ) d τ (cid:12)(cid:12)(cid:12)(cid:12) τ = τ ∗ (cid:33) − = Re (cid:26) p (cid:48) ( λ ) e − λτ + p (cid:48) ( λ ) e − λτ + p (cid:48) ( λ ) λ (2 p ( λ ) e − λτ + p ( λ ) e − λτ ) (cid:27) τ = τ ∗ − Re (cid:110) τ λ (cid:111) τ = τ ∗ = Re (cid:26) p (cid:48) ( iω ) e − iω τ + p (cid:48) ( iω ) e − iω τ + p (cid:48) ( iω ) iω (2 p ( iω ) e − iω τ + p ( iω ) e − iω τ ) (cid:27) = 1 ω Im (cid:26) p (cid:48) ( iω ) e − iω τ + p (cid:48) ( iω ) e − iω τ + p (cid:48) ( iω )2 p ( iω ) e − iω τ + p ( iω ) e − iω τ (cid:27) = 1Λ ω (cid:104) − x x (cid:48) − y y (cid:48) + x x (cid:48) + y y (cid:48) + ( x y (cid:48) − y x (cid:48) + x y (cid:48) − x (cid:48) y ) sin( ω τ ∗ )+ ( x x (cid:48) + y y (cid:48) − x (cid:48) x − y (cid:48) y ) cos( ω τ ∗ ) + ( x y (cid:48) − x (cid:48) y + x y (cid:48) − x (cid:48) y ) sin(2 ω τ ∗ )+ ( x x (cid:48) + y y (cid:48) − x (cid:48) x − y (cid:48) y ) cos(2 ω τ ∗ ) (cid:105) , where Λ = (cid:12)(cid:12) p ( iω ) e − iω τ + p ( iω ) e − iω τ (cid:12)(cid:12) . Regions of feasibility and stability of the steady states S ∗ and S ∗ with parameter values from Table1 ( a ), and with µ a = 10 ( b ). Black and red curves indicate the boundaries of feasibility and the steady-statebifurcation, whereas dashed lines (blue/brown) show the boundaries of Hopf bifurcation of the steady states S ∗ and S ∗ , respectively, with ‘fH’ indicating the fold-Hopf bifurcation. The first digit of the index refers to S ∗ ,while the second corresponds to S ∗ , and they indicate that in that parameter region the respective steady stateis unfeasible (index is ‘0’), stable (index is ‘1’), unstable via Hopf bifurcation with a periodic solution aroundthis 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 ∗ is also stable.Substituting the values of sin( ω τ ∗ ), cos( ω τ ∗ ), sin(2 ω τ ∗ ), and cos(2 ω τ ∗ ) found earlier gives (cid:32) d Re( λ ) d τ (cid:12)(cid:12)(cid:12)(cid:12) τ = τ ∗ (cid:33) − = 1Λ ω F (cid:48) ( ω )2 p (1)0 ( iω ) = h (cid:48) ( s )Λ p (1)0 ( iω ) . Therefore sgn (cid:40) d Re( λ ) d τ (cid:12)(cid:12)(cid:12)(cid:12) τ = τ ∗ (cid:41) = sgn (cid:32) d Re( λ ) d τ (cid:12)(cid:12)(cid:12)(cid:12) τ = τ ∗ (cid:33) − = sgn (cid:40) h (cid:48) ( s )Λ p (1)0 ( iω ) (cid:41) = sgn[ p (1)0 ( iω ) h (cid:48) ( s )] , which completes the proof.We can now formulate the main result concerning stability of the steady states S ∗ and S ∗ . Theorem 3.2.
Suppose the value of T ∗ reg satisfies conditions (3) and (6). If Equation (7) hasat least one positive root s , and p (1)0 ( iω ) h (cid:48) ( s ) > with ω = √ s , then the steady state S ∗ (respectively, S ∗ ) is stable for ≤ τ < τ ∗ , unstable for τ > τ ∗ , and undergoes a Hopf bifurcationat τ = τ ∗ . Since T ∗ reg satisfies conditions (3) and (6), the steady state S ∗ /S ∗ is stable for τ = 0. Lemma3.1 then ensures that τ ∗ is the first positive value of the time delay τ , for which the roots of thecharacteristic Equation (4) cross the imaginary axis with positive speed. Hence, the steady state S ∗ /S ∗ is stable for 0 ≤ τ < τ ∗ , unstable for τ > τ ∗ , and undergoes a Hopf bifurcation at τ = τ ∗ . Remark.
A similar result can be formulated for a subcritical Hopf bifurcation of the steady state S ∗ /S ∗ at some higher value of τ . Parameter Value Parameter Value β ρ µ a d n d F d a µ F δ d in δ α σ λ r σ d r d i p τ p τ ρ τ ρ The only remaining steady state is the persistent (chronic) equilibrium S ∗ with all of its com-ponents being positive. Since it did not prove possible to find a closed form expression for thissteady state, its stability also has to be studied numerically. To investigate the role of different parameters in the dynamics of model (1), in this section we per-form a detailed numerical bifurcation analysis and simulations of this model. Stability of differentsteady states is determined numerically by computing the largest real part of the characteristiceigenvalues, which is achieved by using a pseudospectral method implemented in a traceDDE suitein MATLAB [61].Analytical results from the previous section suggest that at β = d F , the disease-free steady state S ∗ undergoes a transcritical bifurcation. For β < d F , the disease-free steady state S ∗ is stable, andthe chronic steady state is infeasible. On the contrary, when β > d F , the disease-free steady state S ∗ 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 valuesas given in Table 1. For this choice of parameters, we have d F − β >
0, implying that S ∗ is alwaysstable, and Figure 2 illustrates how the stability of S ∗ and S ∗ is affected by parameters. This figureindicates that the steady states S ∗ and S ∗ are only biologically feasible if the regulatory T cells donot grow too rapidly and do not clear autoreactive T cells too quickly. Importantly, Figure 2 showsthat the value of the rate δ of clearance of IL-2 by regulatory T cells does not have any effect onthe thresholds of λ r and δ , where the steady states S ∗ and S ∗ lose their feasibility. Moreover,if λ r and δ are small, then increasing the rate δ at which Tregs inhibit the production of IL-2makes S ∗ become unfeasible, resulting in a stable steady state S ∗ , which has the zero populationof host cells A . On the other hand, the steady state S ∗ associated with autoimmune responses isfavoured for higher values of δ and λ r . In the case stable periodic solutions around these steadystates, increasing δ results in the disappearance of oscillations and stabilisation of the associatedsteady 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 (alsoknown as zero-Hopf or saddle-node Hopf) bifurcation [62].Since our earlier analysis showed that stability of the steady states S ∗ /S ∗ is affected by the timedelay τ , in Figure 3 we consider stability of these equilibria depending on τ and the rate δ . Forthe steady state S ∗ , if the effect of IL-2 on promoting proliferation of T cells is fast (i.e., τ is small),there is a large range of δ , starting with some very low values, for which S ∗ is stable. Increasing Stability of the steady states S ∗ ( a ), and S ∗ ( b ) with parameter values from Table 1. White areashows the region where the steady state S ∗ is infeasible. Colour code denotes max[Re( λ )] for the steady stateswhen they are feasible. In all plots the condition d F > β holds, so the disease-free steady state S ∗ is stable.Basins of attraction of different steady states depending on the initial conditions ( c ), with other initial conditionsspecified in (8), and parameter values from Table 1, except for τ = 18. Cyan and pink areas are the basins ofattraction of S ∗ and S ∗ , respectively.the time delay τ results in the Hopf bifurcation of this steady state as described in Theorem 3.2.One should note that for intermediate values of δ , the steady state S ∗ undergoes stability switches,whereby increasing the delay τ further results in a subcritical Hopf bifurcation, which stabilises S ∗ , but after some number of such stability switches eventually the steady state S ∗ is unstable.For higher still values of δ , the steady state S ∗ remains stable for an entire range of τ values,and the only way to lose its stability is via a steady state bifurcation as given by (5). In the caseof autoimmune steady state S ∗ , the situation is somewhat different in that increasing δ beyondsome critical values makes this steady state biologically infeasible. At the same time, for an entirerange of δ values where it is feasible, this steady state exhibits a single loss of stability through aHopf bifurcation for some critical value of the time delay τ , in agreement with Theorem 3.2.As mentioned earlier, for parameter values used in Figure 3, the disease-free steady state S ∗ isstable. Hence, the system exhibits a bi-stability between a disease-free state and either stable steadystates S ∗ /S ∗ , or periodic solutions around these steady states. To investigate this bi-stability, wechoose parameter values as in Table 1 except for τ = 18, which corresponds to a stable steadystate S ∗ , and we fix initial conditions for state variables as follows,( A (0) , T in (0) , T nor (0) , T aut (0) , I (0)) = (0 . , . , , , , (8)except for initial amounts of infected cells and regulatory T cells that are allowed to vary. Figure 3cillustrates the bi-stability between S ∗ and S ∗ in terms of their basins of attraction. It is worth notingthat recently significant research in approximation theory and meshless interpolation has focusedon 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 tothe disease-free steady state. It also indicates that if the initial amount of infected cells is very smallor is bigger than some specific value, then the infection will be cleared. Interestingly, increasingthe 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 ∗ . In Figure 3b we discovered thatincreasing τ makes the autoimmune steady state S ∗ undergo a Hopf bifurcation, in which casethe system will exhibit a bi-stability between stable S ∗ and a periodic solution around S ∗ . Ournumerical investigation suggests that the shape of basins of attraction in this case is qualitativelysimilar to that shown in Figure 3c, with the basin of attraction of the stable steady state S ∗ beingreplaced 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 astable disease-free steady state and a stable autoimmune steady state S ∗ (similar pattern of be- Numerical solutions of the model with parameters values from Table 1, except for τ = 18. ( a , b ) Sta-ble disease-free steady state S ∗ for F (0) = 0 .
18, and T reg (0) = 100. ( c , d ) Transient oscillations settling on astable steady state S ∗ for F (0) = 0 .
18, and T reg (0) = 10. ( e , f ) Autoimmune dynamics represented by periodicoscillations around the steady state S ∗ for τ = 32, F (0) = 0 .
18, and T reg (0) = 10.haviour is exhibited in the case of bi-stability between S ∗ and S ∗ ). It also illustrates how thesystem develops a periodic solution around the steady state S ∗ for a higher value of τ . Periodicoscillations around the steady state S ∗ biologically correspond to a genuine autoimmune state:after the initial infection is cleared, the system exhibits sustained endogenous oscillations, charac-terised by periods of significant reduction in the number of organ cells through a negative action ofautoreactive T cells, separated by periods of quiescence. This type of behaviour is often observed inclinical manifestations of autoimmune disease [44, 45, 46, 47]. This result has substantial biologicalsignificance 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 withoutlasting consequences, or progression to autoimmunity, also depends on the strength of the initialinfection and of the initial state of the immune system, as represented by the initial number ofregulatory T cells.Next we consider a situation where β > d F , so the disease-free steady state is unstable, and thesystem can have three steady states S ∗ , S ∗ and S ∗ . Our earlier results [48] suggest that in the casewhere regulatory T cells do not inhibit the production of IL-2, i.e., for δ = 0, the steady state S ∗ is stable. Figure 5 shows regions of feasibility and stability of these steady states depending on δ and τ in this case. One observes that S ∗ and S ∗ , whose stability boundaries are determinedby Theorem 3.2, exhibit the same behaviour as in Figure 3, namely, for S ∗ increasing τ causesmultiple stability switches for smaller values of δ , and the steady state is unstable for very small δ and stable for large δ ; in contrast, S ∗ exhibits a single loss of stability via Hopf bifurcation at Stability of S ∗ ( a ), S ∗ ( b ), and S ∗ ( c ), with parameter values from Table 1, except for β = 1 . σ = 1, so that β > d F . White area shows the region where the steady state is infeasible. Colour code denotesmax[Re( λ )] for each steady states when it is feasible. ( d ) Summary of stability results. Green indicates theregion where S ∗ and S ∗ are stable, and S ∗ is unstable, whereas red is the area where S ∗ and S ∗ are stable, and S ∗ is infeasible. Yellow is where S ∗ is stable, S ∗ is unstable, and S ∗ is infeasible. Purple shows the region where S ∗ is stable, but S ∗ and S ∗ are unstable. Blue and cyan indicate the regions where S ∗ and S ∗ are unstable,but S ∗ is stable or unstable, respectively. Bi-stability analysis of the steady states S ∗ , S ∗ , and S ∗ with the same parameter values as inFigure 5, except for ( a ) δ = 0 .
1, ( b ) δ = 0 .
02. Yellow indicates the basin of attraction of the chronic steadystate S ∗ , purple is the basin of attraction of periodic solutions around S ∗ . Red and pink are the basins ofattraction of the steady states S ∗ and S ∗ , respectively.some critical value of the time delay τ , which itself increases with δ . Behaviour of S ∗ is similarto that of S ∗ in that there are multiple stability switches for increasing value of τ and small tointermediate values of δ , while for high values of δ , the chronic steady state S ∗ is stable forall values of τ . Figure 5d divides the δ - τ plane into different regions based on feasibility andstability of these steady states and shows that increasing δ makes the autoimmune steady state S ∗ infeasible. In other regions, the system can exhibit a bi-stability between a stable steady state S ∗ and either a stable steady state S ∗ , or a periodic solution around S ∗ .Figure 6 illustrates the basins of attraction of the steady states S ∗ , S ∗ and S ∗ , as well as periodicsolutions around S ∗ . Figure 6a shows the basins of attraction of the steady states S ∗ and S ∗ anddemonstrates that if the initial number of regulatory T cells or infected cells is sufficiently high, orthe initial amount of infected cells is very low, the immune response neither eliminates infectionnor clears autoreactive T cells, and the system approaches the stable steady state S ∗ . Figure 6billustrates bi-stability between the stable steady state S ∗ and a periodic solution around S ∗ , andhas a different behaviour to than shown in Figure 6a. This figure suggests that for a specific rangeof F (0) the system converges to a stable autoimmune state S ∗ for all values of T reg (0). However, ifthe initial number of infected cells is very high or very low, the system instead develops a periodicsolution around the steady state S ∗ associated with chronic infection.Figure 7 illustrates a regime of bi-stability between a stable steady state S ∗ and a periodicsolution around S ∗ for combinations of initial conditions indicated by crossed in Figure 6b. It alsoillustrates how the system develops a stable solution around the steady state S ∗ for a higher valueof τ . This figure shows that by increasing the initial number of infected cells the behaviour ofthe system changes, as it then approaches the autoimmune steady state S ∗ . Interestingly, one canobserve that for high values of F (0) the system can eliminate the infection, but it cannot clear theautoreactive T cells, in which case the system converges to S ∗ . On the other hand, for a smallernumber 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 ∗ changes with respectto time delays. Figure 8a indicates that for small values of τ (i.e., when the influence of IL-2 on Numerical solutions of the model with same parameters values as Figure 6b. ( a , b ) Stable steady state S ∗ for F (0) = 0 .
002 and T reg (0) = 200. ( c , d ) Periodic oscillations around the steady state S ∗ for F (0) = 0 . T reg (0) = 200. ( e , f ) Transient oscillations settling on a stable steady state S ∗ for τ = 25, F (0) = 0 . T reg (0) = 200.proliferation of T cells is occurring quite rapidly), the steady state S ∗ is stable, and increasing thetime delay τ associated with viral eclipse phase does not have an effect on its stability. At the sametime, if τ exceeds some specific value, by increasing τ the chronic steady state switches betweenbeing stable or unstable. Figure 8b demonstrates a different behaviour, suggesting that for eachvalue of τ , there is small range of τ values where S ∗ is stable, but for smaller and larger valuesof τ it is unstable. For intermediate values of the eclipse phase delay τ , there is an additionalnarrow range of τ values where S ∗ is stable. Figure 8c illustrates that for very small, respectivelyvery large, values of τ , the chronic infection steady state is stable, respectively unstable for anyvalue of τ ; for intermediate values of τ , this steady state undergoes a finite number of stabilityswitches for increasing values of τ and eventually becomes unstable.It should be noted Figure 8 shows that unlike τ and τ , once the steady state S ∗ loses stabilityvia Hopf bifurcation due to increasing τ , it cannot regain stability for higher values of τ . Colour code denotes max[Re( λ )] for the endemic steady state S ∗ depending on τ − τ ( a ), τ − τ ( b ), and τ − τ ( c ), with the parameter values taken from Table 1, except for β = 1 . σ = 1, and δ = 0 . In this paper, we have developed and analysed a time-delayed model of immune response to a viralinfection, which accounts for T cells with different activation thresholds, a cytokine mediating T cellproliferation, as well as regulatory T cells. Particular attention is payed to the dual suppressive roleof regulatory T cells in terms of reducing the amount of autoreactive T cells, and also inhibitingIL-2. To achieve better biological realism of the model, we have explicitly included time delaysassociated 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 systemcan 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 stabilityand 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 immunesystem is able to clear the infection, and the disease-fee steady state is stable. In this regime, thesystem can also support the autoimmune steady state or the steady state with the death of hostcells, 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 otherthree steady states or periodic solutions around them. To better understand bi-stability betweendifferent dynamical regimes, we have used numerical simulations to identify basins of attractionof different steady states and periodic solutions depending on the initial level of infection and theinitial number of regulatory T cells. The fact that for the same parameter values the system canexhibit bi-stability between a disease-free steady state and an autoimmune state, represented bysustained periodic oscillations following the clearance of infection, is very important from a clinicalpoint of view, as effectively it suggests that the progress and eventual outcome of viral infectionis also determined by the strength of infection and the initial state of the immune system. Onecounter-intuitive observation is that in the case of bi-stability with a disease-free steady state, forhigher initial numbers of regulatory T cells, the autoimmune steady state is actually stable for awider range of initial levels of infection. In this regime of bi-stability, increasing the time delayassociated with the positive impact of IL-2 on proliferation of T cells results in the loss of stability ofautoimmune steady state and emergence of autoimmune dynamics, characterised by stable periodicoscillations. On the contrary, in the case where the disease-free steady state is unstable, increasingthis 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 di-rection is exploration of the contributions from other components of immune response, more specif-ically, antibodies and memory T cells, to the onset and progress of autoimmunity [69, 70]. This isparticularly important from the perspective that clinically the onset of autoimmune disease is often aking place on a much longer scale than the timescale of a regular immune response to a viralinfection, so memory T cells can be expected to play a more substantial role. While our model hasfocused 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 ofT cells, as well as mediate their efficiency in eliminating the infection. Including these immunemediators explicitly in the model can provide further significant insights into the dynamics of im-mune response, as has been recently demonstrated on the example of a detailed model of immuneresponse to hepatitis B [72]. Another biologically relevant and mathematically challenging problemis the investigation of the interplay between stochasticity, which is known to be an intrinsic featureof immune response [73, 49], and effects of time delays associated with various aspects of immunedynamics. Appendix
Coefficients of Equation (7) for Hopf frequency are given below. b = (cid:16) − δ δ ρ Y + 2 d r δ δ (cid:17) T ∗ reg + (cid:16) − δ δ ρ (cid:16) d a δ ρ + d i δ ρ − d r δ ρ (cid:17) I ∗ +8 d r δ δ (cid:16) d a δ + d i δ (cid:17)(cid:17) T ∗ reg + (cid:16) − δ δ ρ ρ (cid:16) δ ρ + δ ρ + 5 δ ρ (cid:17) I ∗ − δ δ (cid:16) d a δ ρ +2 d a d i δ δ ρ − d a d r δ ρ ρ + d i δ ρ − d i d r δ δ ρ ρ +2 d r δ ρ +2 d r δ δ ρ +6 d r δ δ ρ ρ + d r δ ρ +8 d r δ ρ (cid:17) I ∗ +2 d r δ δ (cid:16) d a δ +16 d a d i δ δ +6 d i δ + d r δ + d r δ (cid:17)(cid:17) T ∗ reg + (cid:16) − ρ δ ρ δ (cid:16) d a δ δ ρ + 6 d a δ ρ ρ + 9 d i δ ρ + 5 d i δ δ ρ + 22 d i δ δ ρ ρ + 6 d r δ ρ +2 d r δ δ ρ − d r δ ρ (cid:17) I ∗ +2 δ δ (cid:16) d a d i δ δ ρ +2 d a d r δ ρ ρ + d a d i δ δ ρ +10 d a d i d r δ δ ρ ρ − d a d r δ δ ρ − d a d r δ δ ρ − d a d r δ δ ρ ρ − d a d r δ ρ +6 d i d r δ δ ρ ρ − d i d r δ ρ − d i d r δ δ ρ − d i d r δ δ ρ ρ − d i d r δ δ ρ − d i d r δ δ ρ ρ − d i d r δ δ ρ + d r δ δ ρ ρ +2 d r δ δ ρ ρ +2 d r δ ρ ρ (cid:17) I ∗ +4 d r δ δ (cid:16) d a δ +12 d a d i δ δ +12 d a d i δ δ +2 d a d r δ δ + d a d r δ +2 d i δ + d i d r δ +2 d i d r δ δ (cid:17)(cid:17) T ∗ reg + (cid:16) − δ ρ ρ (cid:16) δ ρ + δ δ ρ +3 δ δ ρ ρ +2 δ ρ (cid:17) I ∗ − ρ δ ρ (cid:16) d a δ δ ρ + d a δ ρ ρ +9 d a d i δ δ ρ +24 d a d i δ δ ρ ρ +12 d a d r δ δ ρ − d a d r δ ρ +4 d i δ ρ +8 d i δ δ ρ +33 d i δ δ ρ ρ +9 d i d r δ ρ +9 d i d r δ δ ρ +2 d i d r δ δ ρ − d i d r δ δ ρ − d r δ δ ρ ρ − d r δ δ ρ ρ − d r δ δ ρ + d r δ ρ ρ − d r δ ρ (cid:17) I ∗ + (cid:16) d a d i δ δ ρ + 6 d a d i δ δ ρ + 4 d a d i d r δ δ ρ ρ − d a d r δ δ ρ − d a d r δ δ ρ − d a d r δ δ ρ ρ − d a d r δ ρ + 2 d a d i δ δ ρ + 12 d a d i d r δ δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ ρ +2 d a d i d r δ δ ρ − d a d i d r δ δ ρ ρ − d a d i d r δ δ ρ +2 d a d r δ δ ρ ρ +8 d a d r δ δ ρ ρ +4 d i d r δ δ ρ ρ − d i d r δ ρ − d i d r δ δ ρ − d i d r δ δ ρ ρ − d i d r δ δ ρ − d i d r δ δ ρ ρ − d i d r δ δ ρ +10 d i d r δ δ ρ ρ +12 d i d r δ δ ρ ρ − d r δ δ ρ − d r δ ρ (cid:17) I ∗ + 2 d r (cid:16) d a δ + 16 d a d i δ δ + 36 d a d i δ δ + 6 d a d r δ δ + d a d r δ +16 d a d i δ δ +8 d a d i d r δ δ +8 d a d i d r δ δ + d i δ + d i d r δ +6 d i d r δ δ (cid:17)(cid:17) T ∗ reg + (cid:16) − δ ρ ρ (cid:16) d a δ δ ρ + 2 d a δ ρ ρ + 2 d i δ ρ + 5 d i δ δ ρ + 14 d i δ δ ρ ρ + 2 d i δ ρ ρ +12 d i δ ρ − d r δ δ ρ ρ − d r δ ρ ρ (cid:17) I ∗ +2 ρ ρ (cid:16) d a d i δ δ ρ − d a d i δ ρ ρ − d a d r δ δ ρ +4 d a d i δ δ ρ +2 d a d i δ δ ρ − d a d i δ δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ +20 d a d i d r δ ρ +8 d a d r δ δ ρ ρ +4 d a d r δ ρ − d i δ δ ρ − d i δ δ ρ ρ − d i d r δ ρ − d i d r δ δ ρ − d i d r δ δ ρ +36 d i d r δ δ ρ +6 d i d r δ δ ρ ρ +20 d i d r δ δ ρ ρ +12 d i d r δ δ ρ − d i d r δ ρ ρ +20 d i d r δ ρ +3 d r δ δ ρ − d r δ ρ (cid:17) I ∗ + (cid:16) d a d i δ δ ρ − d a d i d r δ ρ ρ − d a d r δ δ ρ − d a d r δ ρ ρ +2 d a d i δ δ ρ − d a d i d r δ δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ ρ − d a d i d r δ ρ ρ − d a d i d r δ ρ − d a d r δ δ ρ ρ +4 d a d r δ ρ ρ − d a d i d r δ δ ρ ρ − d a d i d r δ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ ρ +2 d a d i d r δ δ ρ − d a d i d r δ δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ ρ +20 d a d i d r δ δ ρ ρ − d a d i d r δ ρ ρ − d a d r δ δ ρ − d i d r δ ρ − d i d r δ ρ ρ − d i d r δ δ ρ − d i d r δ δ ρ ρ − d i d r δ δ ρ − d i d r δ ρ ρ + 8 d i d r δ δ ρ ρ + 12 d i d r δ δ ρ ρ − d i d r δ δ ρ − d i d r δ ρ (cid:17) I ∗ +8 d r (cid:16) d a d i δ + 6 d a d i δ δ + d a d r δ δ + 6 d i δ δ d a + 3 d a d i d r δ δ + d a d i d r δ + d a d i δ + d a d i δ d r +3 d a d i δ d r δ + d i δ δ d r (cid:17)(cid:17) T ∗ reg + (cid:16) δ ρ ρ I ∗ − ρ ρ (cid:16) d a δ ρ +6 d a d i δ δ ρ + 18 d a d i δ ρ ρ + 2 d a d r δ ρ ρ + 4 d i δ ρ + 16 d i δ δ ρ + 40 d i δ δ ρ ρ + d i δ ρ + 20 d i δ ρ ρ + 52 d i δ ρ + 4 d i d r δ δ ρ ρ − d i d r δ ρ ρ + 5 d r δ ρ (cid:17) I ∗ +2 ρ ρ (cid:16) d a d i δ ρ − d a d r δ ρ +12 d a d i δ δ ρ + d a d i δ ρ +3 d a d i δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ ρ +4 d a d r δ ρ ρ +4 d a d i δ ρ +6 d a d i δ δ ρ − d a d i δ δ ρ ρ − d a d i d r δ δ ρ + d a d i d r δ ρ +36 d a d i d r δ ρ +14 d a d i d r δ δ ρ ρ +12 d a d i d r δ ρ +3 d a d r δ ρ − d i δ ρ ρ − d i d r δ ρ − d i d r δ δ ρ +28 d i d r δ δ ρ +3 d i d r δ ρ ρ +22 d i d r δ δ ρ ρ +12 d i d r δ δ ρ +3 d i d r δ ρ ρ +36 d i d r δ ρ +9 d i d r δ δ ρ − d i d r δ ρ (cid:17) I ∗ + (cid:16) − d a d i δ ρ − d a d r δ ρ − d a d i δ δ ρ − d a d i d r δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ ρ ρ − d a d r δ ρ ρ − d a d i δ ρ − d a d i d r δ δ ρ ρ − d a d i d r δ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ ρ − d a d i d r δ ρ − d a d i d r δ ρ ρ − d a d i d r δ ρ − d a d i d r δ δ ρ ρ +10 d a d i d r δ ρ ρ − d a d r δ ρ − d a d i d r δ ρ ρ − d a d i d r δ ρ − d a d i d r δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ ρ ρ +16 d a d i d r δ δ ρ ρ − d a d i d r δ ρ ρ − d a d i d r δ δ ρ − d i d r δ ρ − d i d r δ ρ ρ − d i d r δ ρ +2 d i d r δ ρ ρ +4 d i d r δ δ ρ ρ − d i d r δ ρ − d i d r δ ρ (cid:17) I ∗ +2 d r (cid:16) d a d i δ + d a d r δ +16 d a d i δ δ +8 d a d i d r δ δ +6 d a d i δ +6 d a d i d r δ +6 d a d i d r δ +8 d a d i d r δ δ + d i d r δ (cid:17)(cid:17) T ∗ reg + (cid:16) d i δ ρ ρ I ∗ − d i ρ ρ (cid:16) d a δ ρ +2 d a d i δ ρ +6 d a d i δ ρ ρ +3 d a d r δ ρ ρ +2 d i δ ρ +4 d i δ ρ ρ + d i δ ρ +8 d i δ ρ ρ + 12 d i δ ρ + 2 d i d r δ ρ ρ − d i d r δ ρ ρ + 3 d r δ ρ (cid:17) I ∗ + 2 d i ρ ρ (cid:16) d a d i δ ρ + d a d r δ ρ + 5 d a d i δ ρ + 3 d a d i δ ρ + 8 d a d i δ ρ ρ + 5 d a d i d r δ ρ − d a d i d r δ ρ +8 d a d r δ ρ ρ +4 d a d i δ ρ − d a d i d r δ ρ +3 d a d i d r δ ρ +28 d a d i d r δ ρ +8 d a d i d r δ ρ ρ +12 d a d i d r δ ρ + 9 d a d r δ ρ − d i d r δ ρ + 8 d i d r δ ρ + 8 d i d r δ ρ ρ + 4 d i d r δ ρ +8 d i d r δ ρ ρ +28 d i d r δ ρ +6 d i d r δ ρ − d i d r δ ρ (cid:17) I ∗ − d i (cid:16) d a d i δ ρ + d a d r δ ρ + d a d i δ ρ +6 d a d i d r δ ρ ρ +3 d a d i d r δ ρ +4 d a d i d r δ ρ ρ +4 d a d r δ ρ ρ +4 d a d i d r δ ρ ρ + d a d i d r δ ρ + 4 d a d i d r δ ρ ρ + 2 d a d i d r δ ρ + 6 d a d i d r δ ρ ρ + 16 d a d i d r δ ρ +7 d a d i d r δ ρ ρ − d a d i d r δ ρ ρ +5 d a d r δ ρ + d a d i d r δ ρ +4 d a d i d r δ ρ ρ +8 d a d i d r δ ρ − d a d i d r δ ρ ρ + 6 d a d i d r δ ρ ρ + 5 d a d i d r δ ρ + 8 d i d r δ ρ (cid:17) I ∗ + 4 d a d i d r (cid:16) d a d i δ + d a d r δ +2 d a d i δ +2 d a d i d r δ +2 d a d i d r δ + d i d r δ (cid:17)(cid:17) T ∗ reg +8 d i ρ ρ I ∗ − d i ρ ρ (cid:16) d a ρ +8 d a d i ρ ρ +12 d a d r ρ ρ +4 d i ρ +16 d i ρ ρ +16 d i ρ − d i d r ρ ρ +8 d r ρ (cid:17) I ∗ +2 d i ρ ρ (cid:16) d a d i ρ +3 d a d r ρ + 2 d a d i ρ + 4 d a d i ρ ρ − d a d i d r ρ + 5 d a d r ρ ρ + 2 d a d i d r ρ + 8 d a d i d r ρ +4 d a d i d r ρ +6 d a d r ρ +4 d i d r ρ ρ +8 d i d r ρ − d i d r ρ (cid:17) I ∗ − d i (cid:16) d a d i ρ +2 d a d r ρ +4 d a d i d r ρ ρ +2 d a d i d r ρ ρ +6 d a d r ρ ρ +2 d a d i d r ρ +4 d a d i d r ρ ρ +8 d a d i d r ρ − d a d i d r ρ ρ + 5 d a d r ρ + 4 d a d i d r ρ ρ + 4 d i d r ρ (cid:17) I ∗ + 2 d a d i d r (cid:16) d a d i + d a d r + d i d r (cid:17) . b = T ∗ reg δ δ + 4 δ δ ( d a δ + d i δ ) T ∗ reg + (cid:16) − δ δ (cid:16) δ ρ + 2 δ δ ρ + 6 δ δ ρ ρ + δ ρ + 4 δ ρ (cid:17) I ∗ + 2 δ δ (cid:16) d a δ + 8 d a d i δ δ + 3 d i δ + 2 d r δ + 2 d r δ (cid:17)(cid:17) T ∗ reg + (cid:16) − δ δ (cid:16) d a δ δ ρ + 2 d a δ δ ρ + 8 d a δ δ ρ ρ + 4 d a δ ρ + d i δ ρ + 3 d i δ δ ρ +8 d i δ δ ρ ρ + d i δ δ ρ +2 d i δ δ ρ ρ +8 d i δ δ ρ − d r δ δ ρ ρ − d r δ δ ρ ρ − d r δ ρ ρ (cid:17) I ∗ + δ δ (cid:16) d a δ + 6 d a d i δ δ + 6 d a d i δ δ + 4 d a d r δ δ + 2 d a d r δ + d i δ + 2 d i d r δ +4 d i d r δ δ (cid:17)(cid:17) T ∗ reg + (cid:16) ρ δ ρ (cid:16) δ ρ + 4 δ δ ρ ρ − δ δ ρ + 6 δ δ ρ ρ + 4 δ δ ρ − δ ρ ρ + 4 δ ρ (cid:17) Y + (cid:16) − d a δ δ ρ − d a δ δ ρ − d a δ δ ρ ρ − d a δ ρ − d a d i δ δ ρ − d a d i δ δ ρ − d a d i δ δ ρ ρ +2 d a d i δ δ ρ − d a d i δ δ ρ ρ − d a d i δ δ ρ +2 d a d r δ δ ρ ρ + 8 d a d r δ δ ρ ρ − d i δ ρ − d i δ δ ρ − d i δ δ ρ ρ − d i δ δ ρ − d i δ δ ρ ρ − d i δ δ ρ +10 d i d r δ δ ρ ρ +12 d i d r δ δ ρ ρ − d r δ ρ − d r δ δ ρ − d r δ δ ρ ρ − d r δ δ ρ − d r δ δ ρ ρ − d r δ δ ρ − d r δ δ ρ − d r δ δ ρ ρ − d r δ ρ (cid:17) I ∗ + d a δ +16 d a d i δ δ +36 d a d i δ δ +24 d a d r δ δ +4 d a d r δ +16 d a d i δ δ +32 d a d i d r δ δ + 32 d a d i d r δ δ + d i δ + 4 d i d r δ + 24 d i d r δ δ + d r δ + 4 d r δ δ + d r δ (cid:17) T ∗ reg + (cid:16) ρ ρ (cid:16) d a δ δ ρ − d a δ δ ρ + 8 d a δ δ ρ ρ + 4 d a δ ρ + 3 d i δ ρ +3 d i δ δ ρ +6 d i δ δ ρ ρ − d i δ δ ρ +20 d i δ δ ρ ρ +12 d i δ δ ρ − d i δ ρ ρ +20 d i δ ρ + d r δ ρ − d r δ δ ρ − d r δ δ ρ +3 d r δ δ ρ − d r δ ρ (cid:17) I ∗ + (cid:16) − d a δ δ ρ − d a δ ρ ρ − d a d i δ δ ρ − d a d i δ δ ρ − d a d i δ δ ρ ρ − d a d i δ ρ ρ − d a d i δ ρ − d a d r δ δ ρ ρ +4 d a d r δ ρ ρ − d a d i δ ρ − d a d i δ δ ρ − d a d i δ δ ρ ρ +2 d a d i δ δ ρ − d a d i δ δ ρ ρ − d a d i δ δ ρ − d a d i d r δ δ ρ ρ +20 d a d i d r δ δ ρ ρ − d a d i d r δ ρ ρ − d a d r δ ρ − d a d r δ δ ρ − d a d r δ δ ρ ρ − d a d r δ δ ρ − d a d r δ δ ρ ρ − d a d r δ δ ρ − d a d r δ ρ ρ − d i δ ρ − d i δ ρ ρ − d i δ δ ρ − d i δ δ ρ ρ − d i δ δ ρ − d i d r δ ρ ρ + 8 d i d r δ δ ρ ρ +12 d i d r δ δ ρ ρ − d i d r δ ρ − d i d r δ ρ ρ − d i d r δ δ ρ − d i d r δ δ ρ ρ − d i d r δ δ ρ − d i d r δ δ ρ − d i d r δ δ ρ ρ − d i d r δ ρ ρ − d i d r δ ρ − d r δ ρ ρ +4 d r δ δ ρ ρ +2 d r δ δ ρ ρ + 4 d r δ ρ ρ (cid:17) I ∗ + 4 d a d i δ + 24 d a d i δ δ + 16 d a d r δ δ + 24 d i δ δ d a +48 d a d i d r δ δ + 16 d a d i d r δ + 4 d a d i δ + 16 d a d i δ d r + 48 d a d i δ d r δ + 4 d a d r δ +8 d a d r δ δ + 16 d i δ δ d r + 8 d i d r δ δ + 4 d i d r δ (cid:17) T ∗ reg + (cid:16) − ρ ρ (cid:16) δ ρ + 2 δ δ ρ +6 δ δ ρ ρ + δ ρ + 2 δ ρ ρ + 5 δ ρ (cid:17) I ∗ + 2 ρ ρ (cid:16) d a δ δ ρ − d a δ ρ + 4 d a δ ρ ρ +7 d a d i δ ρ +2 d a d i δ δ ρ +14 d a d i δ δ ρ ρ +2 d a d i δ ρ +12 d a d i δ ρ +3 d a d r δ ρ − d a d r δ δ ρ − d a d r δ ρ +3 d a d r δ ρ +4 d i δ ρ +3 d i δ ρ ρ +4 d i δ δ ρ +22 d i δ δ ρ ρ +12 d i δ δ ρ + d i δ ρ + 3 d i δ ρ ρ + 36 d i δ ρ − d i d r δ ρ − d i d r δ δ ρ + 9 d i d r δ δ ρ − d i d r δ ρ − d i d r δ ρ + d r δ ρ ρ + 2 d r δ δ ρ ρ + 2 d r δ δ ρ + 4 d r δ ρ ρ + 3 d r δ ρ (cid:17) I ∗ + (cid:16) − d a δ ρ − d a d i δ δ ρ − d a d i δ ρ ρ − d a d r δ ρ ρ − d a d i δ ρ − d a d i δ δ ρ − d a d i δ δ ρ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i δ ρ − d a d i d r δ δ ρ ρ +10 d a d i d r δ ρ ρ − d a d r δ ρ − d a d r δ δ ρ − d a d r δ δ ρ ρ − d a d r δ ρ − d a d r δ ρ ρ − d a d r δ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i δ δ ρ − d a d i δ δ ρ ρ − d a d i δ δ ρ − d a d i d r δ ρ ρ + 16 d a d i d r δ δ ρ ρ − d a d i d r δ ρ ρ − d a d i d r δ ρ − d a d i d r δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ δ ρ ρ − d a d i d r δ δ ρ − d a d i d r δ ρ ρ − d a d r δ ρ ρ +8 d a d r δ δ ρ ρ − d a d r δ ρ ρ − d i δ ρ − d i δ ρ ρ − d i δ ρ +2 d i d r δ ρ ρ +4 d i d r δ δ ρ ρ − d i d r δ ρ − d i d r δ ρ ρ − d i d r δ ρ − d i d r δ δ ρ − d i d r δ δ ρ ρ − d i d r δ ρ − d i d r δ ρ ρ − d i d r δ ρ +2 d i d r δ ρ ρ +10 d i d r δ ρ ρ − d r δ ρ − d r δ ρ (cid:17) I ∗ +6 d a d i δ +4 d a d r δ +16 d a d i δ δ +32 d a d i d r δ δ +6 d a d i δ +24 d a d i d r δ +24 d a d i d r δ +6 d a d r δ +4 d a d r δ +32 d a d i d r δ δ +16 d a d i d r δ δ +4 d i d r δ +4 d i d r δ +6 d i d r δ (cid:17) T ∗ reg + (cid:16) − ρ ρ (cid:16) d a δ ρ +2 d a δ ρ ρ + d i δ ρ +2 d i δ ρ ρ +2 d i δ ρ +6 d i δ ρ ρ +6 d i δ ρ + d r δ ρ ρ − d r δ ρ ρ (cid:17) I ∗ + 2 ρ ρ (cid:16) d a δ ρ + 5 d a d i δ ρ + d a d i δ ρ + 8 d a d i δ ρ ρ + 3 d a d r δ ρ − d a d r δ ρ +6 d a d i δ ρ +8 d a d i δ ρ ρ +4 d a d i δ ρ +12 d a d i δ ρ − d a d i d r δ ρ + d a d i d r δ ρ +9 d a d i d r δ ρ + 2 d a d r δ ρ ρ + 2 d a d r δ ρ + 3 d i δ ρ + 8 d i δ ρ ρ + 4 d i δ ρ + 3 d i δ ρ +8 d i δ ρ ρ + 28 d i δ ρ − d i d r δ ρ + 6 d i d r δ ρ − d i d r δ ρ − d i d r δ ρ + 2 d i d r δ ρ ρ + d i d r δ ρ +8 d i d r δ ρ ρ +9 d i d r δ ρ + d r δ ρ − d r δ ρ (cid:17) I ∗ + (cid:16) − d a d i δ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i d r δ ρ ρ − d a d r δ ρ − d a d r δ ρ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i δ ρ − d a d i d r δ ρ ρ +8 d a d i d r δ ρ ρ − d a d i d r δ ρ − d a d i d r δ ρ ρ − d a d i d r δ ρ − d a d i d r δ ρ ρ − d a d i d r δ ρ − d a d r δ ρ ρ +4 d a d r δ ρ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i δ ρ + 4 d a d i d r δ ρ ρ − d a d i d r δ ρ ρ − d a d i d r δ ρ − d a d i d r δ ρ ρ − d a d i d r δ ρ − d a d i d r δ ρ ρ +4 d a d i d r δ ρ ρ − d a d i d r δ ρ ρ − d a d r δ ρ − d i d r δ ρ − d i d r δ ρ ρ − d i d r δ ρ − d i d r δ ρ ρ − d i d r δ ρ − d i d r δ ρ ρ +8 d i d r δ ρ ρ − d i d r δ ρ (cid:17) I ∗ +4 d a d i δ +8 d a d i d r δ +4 d a d i δ +16 d a d i d r δ +4 d a d r δ +16 d a d i d r δ +8 d a d i d r δ +8 d a d i d r δ +8 d a d i d r δ +4 d i d r δ (cid:17) T ∗ reg + I ∗ ρ ρ − ρ ρ (cid:16) d a ρ +2 d a d i ρ ρ +2 d a d r ρ ρ +5 d i ρ +12 d i ρ ρ +8 d i ρ − d i d r ρ ρ + d r ρ (cid:17) I ∗ +2 ρ ρ (cid:16) d a d i ρ + d a d r ρ +3 d a d i ρ +5 d a d i ρ ρ − d a d i d r ρ + d a d r ρ ρ + d a d i ρ +4 d a d i ρ +3 d a d i d r ρ +6 d a d i d r ρ + d a d i d r ρ + d a d r ρ + 2 d i ρ + 4 d i ρ ρ + 8 d i ρ − d i d r ρ − d i d r ρ +5 d i d r ρ ρ +6 d i d r ρ − d i d r ρ (cid:17) I ∗ + (cid:16) − d a d i ρ − d a d r ρ − d a d i ρ ρ − d a d i d r ρ ρ − d a d i d r ρ ρ − d a d r ρ ρ − d a d i ρ − d a d i ρ ρ − d a d i ρ +2 d a d i d r ρ ρ − d a d i d r ρ − d a d i d r ρ ρ − d a d i d r ρ + 2 d a d i d r ρ ρ − d a d r ρ − d a d i d r ρ ρ − d a d i d r ρ ρ − d a d i d r ρ ρ − d i d r ρ − d i d r ρ ρ − d i d r ρ + 2 d i d r ρ ρ − d i d r ρ (cid:17) I ∗ + d a d i +4 d a d i d r + d a d r + 4 d a d i d r + 4 d a d i d r + d i d r . b = 2 δ δ (cid:16) δ + δ (cid:17) T ∗ reg +4 δ δ (cid:16) d a δ δ + d a δ + d i δ +2 d i δ δ (cid:17) T ∗ reg + (cid:16)(cid:16) − δ ρ − δ δ ρ − δ δ ρ ρ − δ δ ρ − δ δ ρ ρ − δ δ ρ − δ δ ρ − δ δ ρ ρ − δ ρ (cid:17) I ∗ +12 d a δ δ +2 d a δ +16 d a d i δ δ +16 d a d i δ δ +2 d i δ +12 d i δ δ +2 d r δ +8 d r δ δ + 2 d r δ (cid:17) T ∗ reg + (cid:16)(cid:16) − d a δ ρ − d a δ δ ρ − d a δ δ ρ ρ − d a δ δ ρ − d a δ δ ρ ρ − d a δ δ ρ − d a δ ρ ρ − d i δ ρ − d i δ ρ ρ − d i δ δ ρ − d i δ δ ρ ρ − d i δ δ ρ − d i δ δ ρ − d i δ δ ρ ρ − d i δ ρ ρ − d i δ ρ − d r δ ρ ρ +4 d r δ δ ρ ρ +2 d r δ δ ρ ρ +4 d r δ ρ ρ (cid:17) I ∗ +8 d a δ δ +24 d a d i δ δ +8 d a d i δ +8 d a d i δ +24 d a d i δ δ +8 d a d r δ +16 d a d r δ δ +8 d i δ δ +16 d i d r δ δ +8 d i d r δ (cid:17) T ∗ reg + (cid:16) ρ ρ (cid:16) δ ρ + δ ρ ρ +4 δ δ ρ +2 δ δ ρ ρ +2 δ δ ρ − δ ρ +4 δ ρ ρ +3 δ ρ (cid:17) I ∗ + (cid:16) − d a δ ρ − d a δ δ ρ − d a δ δ ρ ρ − d a δ ρ − d a δ ρ ρ − d a δ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i δ δ ρ − d a d i δ δ ρ ρ − d a d i δ δ ρ − d a d i δ ρ ρ − d a d r δ ρ ρ +8 d a d r δ δ ρ ρ − d a d r δ ρ ρ − d i δ ρ − d i δ ρ ρ − d i δ ρ − d i δ δ ρ − d i δ δ ρ ρ − d i δ ρ − d i δ ρ ρ − d i δ ρ +2 d i d r δ ρ ρ +10 d i d r δ ρ ρ − d r δ ρ − d r δ ρ ρ − d r δ ρ − d r δ δ ρ − d r δ δ ρ ρ − d r δ ρ − d r δ ρ ρ − d r δ ρ (cid:17) I ∗ +2 d a δ +16 d a d i δ δ +12 d a d i δ +12 d a d i δ + 12 d a d r δ + 8 d a d r δ + 16 d a d i δ δ + 32 d a d i d r δ δ + 2 d i δ + 8 d i d r δ +12 d i d r δ + 2 d r δ + 2 d r δ (cid:17) T ∗ reg + (cid:16) ρ ρ (cid:16) d a δ ρ + 2 d a δ ρ ρ + 2 d a δ ρ + 2 d a δ ρ +3 d i δ ρ + 2 d i δ ρ ρ + d i δ ρ + d i δ ρ + 8 d i δ ρ ρ + 9 d i δ ρ + d r δ ρ + d r δ ρ − d r δ ρ − d r δ ρ (cid:17) I ∗ + (cid:16) − d a δ ρ − d a δ ρ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i δ ρ − d a d r δ ρ ρ + 4 d a d r δ ρ ρ − d a d i δ ρ − d a d i δ ρ ρ − d a d i δ ρ − d a d i δ ρ ρ + 4 d a d i d r δ ρ ρ − d a d i d r δ ρ ρ − d a d r δ ρ − d a d r δ ρ ρ − d a d r δ ρ − d a d r δ ρ ρ − d i δ ρ − d i δ ρ ρ − d i δ ρ − d i δ ρ ρ − d i δ ρ − d i d r δ ρ ρ +8 d i d r δ ρ ρ − d i d r δ ρ − d i d r δ ρ ρ − d i d r δ ρ − d i d r δ ρ ρ − d i d r δ ρ − d r δ ρ ρ +4 d r δ ρ ρ (cid:17) Y + 4 d a d i δ + 8 d a d i δ + 8 d a d r δ + 8 d a d i δ + 16 d a d i d r δ + 4 d a d i δ +16 d a d i d r δ + 4 d a d r δ + 8 d i d r δ + 4 d i d r δ (cid:17) T ∗ reg − ρ ρ (cid:16) ρ + ρ (cid:17) I ∗ + 2 ρ ρ (cid:16) d a ρ + d a ρ ρ + d a d i ρ + d a d i ρ + d a d r ρ + d a d r ρ + 3 d i ρ + 5 d i ρ ρ + 6 d i ρ − d i d r ρ − d i d r ρ + d r ρ ρ + d r ρ (cid:17) Y + (cid:16) − d a ρ − d a d i ρ ρ − d a d r ρ ρ − d a d i ρ − d a d i ρ ρ − d a d i ρ +2 d a d i d r ρ ρ − d a d r ρ − d a d r ρ ρ − d a d r ρ − d a d i ρ ρ − d a d i d r ρ ρ − d a d i d r ρ ρ − d a d r ρ ρ − d i ρ − d i ρ ρ − d i ρ +2 d i d r ρ ρ − d i d r ρ − d i d r ρ ρ − d i d r ρ + 2 d i d r ρ ρ − d r ρ (cid:17) I ∗ + 2 d a d i + 2 d a d r + 2 d a d i + 8 d a d i d r + 2 d a d r +2 d i d r + 2 d i d r . b = (cid:16) δ +4 δ δ + δ (cid:17) T ∗ reg + (cid:16) d a δ +8 d a δ δ +8 d i δ δ +4 d i δ (cid:17) T ∗ reg + (cid:16)(cid:16) − δ ρ − δ ρ ρ − δ ρ − δ δ ρ − δ δ ρ ρ − δ ρ − δ ρ ρ − δ ρ (cid:17) I ∗ + 6 d a δ + 4 d a δ +16 d a d i δ δ +4 d i δ +6 d i δ +4 d r δ +4 d r δ (cid:17) T ∗ reg + (cid:16)(cid:16) − d a δ ρ − d a δ ρ ρ − d a δ ρ − d a δ ρ ρ − d i δ ρ − d i δ ρ ρ − d i δ ρ − d i δ ρ ρ − d i δ ρ − d r δ ρ ρ +4 d r δ ρ ρ (cid:17) I ∗ +4 d a δ + 8 d a d i δ + 8 d a d i δ + 8 d a d r δ + 4 d i δ + 8 d i d r δ (cid:17) T ∗ reg + 2 ρ ρ (cid:16) ρ + ρ ρ + ρ (cid:17) I ∗ + (cid:16) − d a ρ − d a ρ ρ − d a ρ − d a d i ρ ρ − d a d r ρ ρ − d i ρ − d i ρ ρ − d i ρ +2 d i d r ρ ρ − d r ρ − d r ρ ρ − d r ρ (cid:17) I ∗ + d a + 4 d a d i + 4 d a d r + d i + 4 d i d r + d r . b = (cid:16) δ + 2 δ (cid:17) T ∗ reg + (cid:16) d a δ + 4 d i δ (cid:17) T reg ∗ − (cid:16) ρ + ρ (cid:17) I ∗ + 2 d a + 2 d i + 2 d r . References [1] Mason, D., A very high level of crossreactivity is an essential feature of the T-cell receptor,
Immunol. Today , , 395–404.[2] Anderson, A.C.; Waldner, H.; Turchin, V.; Jabs, C.; Prabhu Das, M.; Kuchroo, V.K.;Nicholson, L.B., Autoantigen responsive T cell clones demonstrate unfocused TCR cross-reactivity towards multiple related ligands: Implications for autoimmunity, Cell. Immunol. , , 88–96.[3] Kerr, E.C.; Copland, D.A.; Dick, A.D.; Nicholson, L.B., The dynamics of leukocyte infiltra-tion in experimental autoimmune uveoretinitis, Prog. Retin. Eye Res. , , 527–535.[4] Prat, E.; Martin, R., The immunopathogenesis of multiple sclerosis, J. Rehabil. Res. Dev. , , 187–200.[5] Santamaria, P., The long and winding road to understanding and conquering type 1 diabetes, Immunity , , 437–445.[6] Root-Bernstein, R.; Fairweather, D., Unresolved issues in theories of autoimmune diseaseusing myocarditis as a framework, J. Theor. Biol. , , 101–123.[7] Caforio, A.L.P.; Iliceto, S., Genetically determined myocarditis: Clinical presentation andimmunological characteristics, Curr. Opin. Cardiol. , , 219–226.[8] Li, H.S.; Ligons, D.L.; Rose, N.R., Genetic complexity of autoimmune myocarditis, Autoim-mun. Rev. , , 168–173.[9] Guilherme, L.; K¨ohler, K.F.; Postol, E.; Kalil, J., Genes, autoimmunity and pathogenesisof rheumatic heart disease, Ann. Pediatr. Cardiol. , , 13–21.[10] Germolic, D.; Kono, D.H.; Pfau, J.C.; Pollard, K.M., Animal models used to examine therole of environment in the development of autoimmune disease: Findings from an NIEHSExpert Panel Workshop, J. Autoimmun. , , 285–293.[11] Mallampalli, M.P.; Davies, E.; Wood, D.; Robertson, H.; Polato, F.; Carter, C.L., Role ofenvironment and sex differences in the development of autoimmune disease: A roundtablemeeting report, J. Women’s Health , , 578–586.[12] Manfredo Vieira, S.; Hiltensperger, M.; Kumar, V.; Zegarra-Ruiz, D.; Dehne, C.; Khan, N.;Costa, F.R.C.; Tiniakou, E.; Greiling, T.; et al., Translocation of a gut pathobiont drivesautoimmunity in mice and humans, Science , , 1156–1161.[13] Fujinami, R.S., Can virus infections trigger autoimmune disease? J. Autoimmun. , , 229–234.
14] Von Herrath, M.G.; Oldstone, M.B.A., Virus-induced autoimmune disease,
Curr. Opin.Immunol. , , 878–885.[15] Ercolini, A.M.; Miller, S.D.,The role of infections in autoimmune disease, Clin. Exp. Im-munol. , , 1–15.[16] Segel, L.A.; J¨ager, E.; Elias, D.; Cohen, I.R., A quantitative model of autoimmune diseaseand T-cell vaccination: Does more mean less? Immunol. Today , , 80–84.[17] Borghans, J.A.M.; De Boer, R.J., A minimal model for T-cell vaccination, Proc. R. Soc.Lond. B Biol. Sci. , , 173–178.[18] Borghans, J.A.M.; De Boer, R.J.; Sercarz, E.; Kumar, V., T cell vaccination in experimentalautoimmune encephalomyelitis: A mathematical model, J. Immunol. , , 1087–1093.[19] Le´on, K.; Perez, R.; Lage, A.; Carneiro, J., Modelling T-cell-mediated suppression depen-dent on interactions in multicellular conjugates, J. Theor. Biol. , , 231–254.[20] Le´on, K.; Lage, A.; Carneiro, J., Tolerance and immunity in a mathematical model of T-cellmediated suppression, J. Theor. Biol. , , 107–126.[21] Le´on, K.; Faro, J.; Lage, A.; Carneiro, J., Inverse correlation between the incidences ofautoimmune disease and infection predicted by a model of T cell mediated tolerance, J.Autoimmun. , , 31–42.[22] Carneiro, J.; Paix˜ao, T.; Milutinovic, D.; Sousa, J.; Leon, K.; Gardner, R.; Faro, J., Im-munological self-tolerance: Lessons from mathematical modeling, J. Comput. Appl. Math. , , 77–100.[23] Alexander, H.K.; Wahl, L.M., Self-tolerance and autoimmunity in a regulatory T cell model, Bull. Math. Biol. , , 33–71.[24] Iwami, S.; Takeuchi, Y.; Miura, Y.; Sasaki, T.; Kajiwara, T., Dynamical properties ofautoimmune disease models: Tolerance, flare-up, dormancy, J. Theor. Biol. , , 646–659.[25] Iwami, S.; Takeuchi, Y.; Iwamoto, K.; Naruo, Y.; Yasukawa, M., A mathematical design ofvector vaccine against autoimmune disease, J. Theor. Biol. , , 382–392.[26] Burroughs, N.J.; Ferreira, M.; Oliveira, B.M.P.M.; Pinto, A.A., Autoimmunity arising frombystander proliferation of T cells in an immune response model, Math. Comput. Model. , , 1389–1393.[27] Burroughs, N.J.; Ferreira, M.; Oliveira, B.M.P.M.; Pinto, A.A., A transcritical bifurcationin an immune response model, J. Differ. Equ. Appl. , , 1101–1106.[28] Burroughs, N.J.; de Oliveira, B.M.P.M.; Pinto, A.A., Regulatory T cell adjustment of quo-rum growth thresholds and the control of local immune responses, J. Theor. Biol. , , 134–141.[29] Sakaguchi, S., Naturally arising CD4 + regulatory T cells for immunologic self-tolerance andnegative control of immune responses, Ann. Rev. Immunol. , , 531–562.[30] Josefowicz, S.Z.; Lu, L.F.; Rudensky, A.Y., Regulatory T cells: Mechanisms of differentia-tion and function, Ann. Rev. Immunol. , , 531–564.[31] Fontenot, J.D.; Gavin, M.A.; Rudensky, A.Y., Foxp3 programs the development and func-tion of CD4 + CD25 + regulatory T cells, Nat. Immunol. , , 330–336.[32] Khattri, R.; Cox, T.; Yasayko, S.A.; Ramsdell, F., An essential role for Scurfin inCD4 + CD25 + T regulatory cells,
Nat. Immunol. , , 337–342.[33] Grossman, Z.; Paul, W.E., Adaptive cellular interactions in the immune system: The tunableactivation threshold and the significance of subthreshold responses, Proc. Natl. Acad. Sci.USA , , 10365–10369.[34] Grossman, Z.; Paul, W.E., Self-tolerance: Context dependent tuning of T cell antigen recog-nition, Semin. Immunol. , , 197–203.
35] Grossman, Z.; Singer, A., Tuning of activation thresholds explains flexibility in the selectionand development of T cells in the thymus,
Proc. Natl. Acad. Sci. USA , , 14747–14752.[36] Altan-Bonnet, G.; Germain, R.N., Modeling T cell antigen discrimination based on feedbackcontrol of digital ERK responses, PLoS Biol. , , e356.[37] Bitmansour, A.D.; Douek, D.C.; Maino, V.C.; Picker, L.J., Direct ex vivo analysis of humanCD4 + memory T cell activation requirements at the single clonotype level, J. Immunol. , , 1207–1218.[38] Nicholson, L.B.; Anderson, A.C.; Kuchroo, V.K., Tuning T cell activation threshold andeffector function with cross-reactive peptide ligands, Int. Immunol. , , 205–213.[39] R¨omer, P.S.; Berr, S.; Avota, E.; Na, S.Y.; Battaglia, M.; ten Berge, I.; Einsele, H.; H¨unig,T., Preculture of PBMC at high cell density increases sensitivity of T-cell responses, reveal-ing cytokine release by CD28 superagonist TGN1412, Blood , , 6772–6782.[40] Stefanov´a, I.; Dorfman, J.R.; Germain, R.N., Self-recognition promotes the foreign antigensensitivity of naive T lymphocytes, Nature , , 429–434.[41] van den Berg, H.A.; Rand, D.A., Dynamics of T cell activation threshold tuning, J. Theor.Biol. , , 397–416.[42] Scherer, A.; Noest, A.; de Boer, R.J., Activation-threshold tuning in an affinity model forthe T-cell repertoire, Proc. R. Soc. B , , 609–616.[43] Blyuss, K.B.; Nicholson, L.B., The role of tunable activation thresholds in the dynamics ofautoimmunity, J. Theor. Biol. , , 45–55.[44] Blyuss, K.B.; Nicholson, L.B., Understanding the roles of activation threshold and infectionsin the dynamics of autoimmune disease, J. Theor. Biol. , , 13–20.[45] Ben Ezra, D.; Forrester, J.V., Fundal white dots: The spectrum of a similar pathologicalprocess, Br. J. Ophthalmol. , , 856–860.[46] Davies, T.F.; Evered, D.C.; Rees Smith, B.; Yeo, P.P.B.; Clark, F.; Hall, R., Value ofthyroid-stimulating-antibody determinations in predicting the short-term thyrotoxic relapsein Graves’ disease, Lancet , , 1181–1182.[47] Nylander, A.; Hafler, D.A., Multiple sclerosis, J. Clin. Investig. , , 1180–1188.[48] Fatehi, F.; Kyrychko, Y.N.; Blyuss, K.B., Interactions between cytokines and T cells in thedynamics of autoimmunity, 2018, submitted.[49] Fatehi, F.; Kyrychko, S.N.; Ross, A.; Kyrychko, Y.N.; Blyuss, K.B., Stochastic effects inautoimmune dynamics, Front. Physiol. , , 45.[50] Wu, H.J.; Ivanov, I.I.; Darce, J.; Hattori, K.; Shima, T.; Umesaki, Y.; Littman, D.R.;Benoist, C.; Mathis, D., Gut-residing segmented filamentous bacteria drive autoimmunearthritis via T helper 17 cells, Immunity , , 815–827.[51] Wolf, S.D.; Dittel, B.N.; Hardardottir, F.; Janeway, C.A., Experimental autoimmune en-cephalomyelitis induction in genetically B cell-deficient mice, J. Exp. Med. , , 2271–2278.[52] Baltcheva, I.; Codarri, L.; Pantaleo, G.; Le Boudec, J.Y., Lifelong dynamics of humanCD4 + CD25 + regulatory T cells: Insights from in vivo data and mathematical modeling, J.Theor. Biol. , , 307–322.[53] Shevach, E.M.; McHugh, R.S.; Piccirillo, C.A.; Thornton, A.M., Control of T-cell activationby CD4 + CD25 + suppressor T cells, Immunol. Rev. , , 58–67.[54] Abbas, A.K.; Lichtman, A.H.; Pillai, S., Cellular and Molecular Immunology ; Elsevier:Philadelphia, USA, 2015.[55] Thornton, A.M.; Shevach, E.M., CD4 + CD25 + immunoregulatory T cells suppress polyclonalT cell activation in vitro by inhibiting interleukin 2 production, J. Exp. Med. , , 287–296.
56] Baccam, P.; Beauchemin, C.; Macken, C.A.; Hayden, F.G.; Perelson, A.S, Kinetics of in-fluenza A virus infection in humans,
J. Virol. , , 7590–7599.[57] Pawelek, K.A.; Huynh, G.T.; Quinlivan, M.; Cullinane, A.; Rong, L.; Perelson, A.S., Mod-eling within-host dynamics of influenza virus infection including immune responses, PLoSComput. Biol. , , e1002588.[58] Kim, P.S.; Lee, P.P.; Levy, D., Modeling regulation mechanisms in the immune system, J.Theor. Biol. , , 33–69.[59] Li, J.; Zhang, L.; Wang, Z., Two effective stability criteria for linear time-delay systemswith complex coefficients, J. Syst. Sci. Complex. , , 835.[60] Rahman, B.; Blyuss, K.B.; Kyrychko, Y.N., Dynamics of neural systems with discrete anddistributed time delays, SIAM J. Appl. Dyn. Syst. , , 2069–2095.[61] Breda, D.; Maset, S.; Vermiglio, R., Numerical computation of characteristic multipliers forlinear time periodic coefficients delay differential equations, IFAC Proc. Vol. , , 163–168.[62] Kuznetsov, Y.A., Elements of Applied Bifurcation Theory ; Springer-Verlag: New York, NY,USA, 1995.[63] Cavoretto, R.; Rossi, A.D.; Perracchione, E.; Venturino, E., Reliable approximation ofseparatrix manifolds in competition models with safety niches,
Int. J. Comput. Math. , , 1826–1837.[64] Cavoretto, R.; De Rossi, A.; Perracchione, E.; Venturino, E., Robust approximation algo-rithms for the detection of attraction basins in dynamical systems, J. Sci. Comput. , , 395–415.[65] De Rossi, A.; Perracchione, E.; Venturino, E., Fast strategy for PU interpolation: Anapplication for the reconstruction of separatrix manifolds, Dolomites Res. Notes Approx. , , 3–12.[66] Francomano, E.; Hilker, F.M.; Paliaga, M.; Venturino, E., On basins of attraction for apredator-prey model via meshless approximation, In AIP Conference Proceedings ; AIP Pub-lishing: Melville, NY, USA, 2016; Volume 1776, p. 070007.[67] Cavoretto, R.; De Rossi, A.; Perracchione, E.; Venturino, E., Graphical representation ofseparatrices of attraction basins in two and three-dimensional dynamical systems,
Int. J.Comput. Methods , , 1750008.[68] Francomano, E.; Hilker, F.M.; Paliaga, M.; Venturino, E., Separatrix reconstruction toidentify tipping points in an eco-epidemiological model, Appl. Math. Comput. , , 80–91.[69] Skapenko, A.; Leipe, J.; Lipsky, P.E.; Schulze-Koops, H., The role of the T cell in autoim-mune inflammation, Arthritis Res. Ther. , (Suppl. 2), S4–S14.[70] Antia, R.; Ganusov, V.V.; Ahmed, R., The role of models in understanding CD8 + T-cellmemory,
Nat. Rev. Immunol. , , 101–111.[71] Schluns, K.S.; Kieper, W.C.; Jameson, S.C.; Lefrancois, L., Interleukin-7 mediates thehomeostasis of na¨ıve and memory CD8 T cells in vivo, Nat. Immunol. , , 426–432.[72] Fatehi Chenar, F.; Kyrychko, Y.N.; Blyuss, K.B., Mathematical model of immune responseto hepatitis B, J. Theor. Biol. , , 98–110.[73] Perelson, A.S.; Weisbuch, G., Immunology for physicists, Rev. Mod. Phys. , , 1219–1267., 1219–1267.