A near analytic solution of a stochastic immune response model considering variability in virus and T cell dynamics
AA near analytic solution of a stochastic immune response modelconsidering variability in virus and T cell dynamics
Abhilasha Batra and Rati Sharma ∗ Department of Chemistry,Indian Institute of Science Education and Research (IISER) BhopalBhopal Bypass Road, Bhauri, Bhopal 462066 INDIA
Biological processes at the cellular level are stochastic in nature, and the immuneresponse system is no different. Therefore, models that attempt to explain this sys-tem need to also incorporate noise or fluctuations that can account for the observedvariability. In this work, a stochastic model of the immune response system is pre-sented in terms of the dynamics of the T cells and the virus particles. Making useof the Green’s function and the Wilemski-Fixman approximation, this model is thensolved to obtain the analytical expression for the joint probability density functionof these variables in the early and late stages of infection. This is then also usedto calculate the average level of virus particles in the system. Upon comparing thetheoretically predicted average virus levels to those of COVID-19 patients, it is hy-pothesized that the long lived dynamics that are characteristic of such viral infectionsare due to the long range correlations in the temporal fluctuations of the virions.This model therefore provides an insight into the effects of noise on viral dynamics.
I. INTRODUCTION
The immune system of an organism provides the protection that it needs against for-eign bodies such as microbes, viruses, parasites and more. The appearance of these foreignentities inside the body triggers the immune response, which is an agglomeration of cells,tissues and biochemical processes that function in conjunction with each other to protectthe body [1]. To stop the proliferation of such antigens, the immune system proliferatesits own cells and shields the body from foreign attack, thereby enabling it to carry out itsregular functions [2].The immune response in humans can be understood as a system with two levels of increas- ∗ [email protected] a r X i v : . [ q - b i o . M N ] F e b ing complexity. These are the innate and the adaptive immunity. Innate immunity is thefirst line of defense and is non-specific. It has no immunologic memory as it cannot dis-tinguish self from non-self [3]. Therefore, after the pathogens manage to evade the innateimmune system, the second line of defense, i.e., the adaptive immune response gets activated.The adaptive/acquired immunity releases antigen-specific response and therefore providesa targeted defense against foreign particles. It works by retaining the copies of antibodiesproduced in the previous attack, therefore encoding the new memory for future use. Thishelps the immune system in launching a faster and targeted attack on foreign bodies in thefuture. The adaptive immunity is provided by a combined effort of two types of lympho-cytes, which are, antigen-specific T cells (matured in thymus) and B cells (matured in thebone marrow) which divide into plasma to produce antibodies. T cells have surface specificantibody-like receptors that can recognize antigens inside the target cell of the host’s bodycarrying the virus and directly destroys them [4].Since T cells are responsible for the directed attack on virus particles, a study of the dynam-ics between these two is pertinent. Theoretical studies of such systems are helpful in givingan insight into the complex dynamical phenomena involving their interaction. Therefore,several groups in the past have developed theoretical models to analyze and predict the virusand T cell dynamics in the system. These include models developed to look at the inte-grated immunological response to different viral infections such as Human ImmunodeficiencyVirus (HIV), Influenza virus, Zika virus and the virus that caused the ongoing pandemic,severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [5–8], among others. Thesetheoretical studies generally use one of the two modelling approaches, viz. models withoutand with the immune response. The models without immune response generally incorporatekinetic interactions between healthy, susceptible and infected cells with the virus particles.Various standard versions of models using this approach, such as, “Target Cell Limitedmodel” [7] and “Target Cell Eclipse Phase model” [9, 10] have been studied. On the otherhand, the models that include the immune response deal with interactions between immunecells (T cells) and viral particles [8, 11]. The immune response model was also developed forthe Influenza A virus, where, the dynamics of cytotoxic T cells and virus population werecoupled through a set of coupled ODEs [6].All of these models discussed above are deterministic. However, biological processes at thecellular level, including the immune response systems are stochastic in nature [12]. Therefore,one must take care to incorporate fluctuations/noise in these systems. In fact, an analysisof experimental data has shown that the presence of random noise in gene expression leadsto increased variability in viral (HIV) gene products such as RNA, which contributes to thereplication of viral material and the latency period [13]. Stochasticity has also been reportedin the division and death time of lymphocytes [14].Taking this variability and population heterogeneity into consideration, a few stochasticimmune response models have been developed in the recent past [15–17]. Dalal et al. intheir work introduced stochasticity in the deterministic model of immune response in HIVinfection by parameter perturbation [15]. Wang et al. showed that the stability of thestochastic dynamical system of HIV infection is different when fluctuations are introducedin terms of the Gaussian colored noise in contrast to the Gaussian white noise [16]. Anotherrecent stochastic model of T cell dynamics was used to explain the bistability and crossoverdynamics of the immune response [18]. These stochastic studies indicate that modeling fluc-tuations into the system can give a more accurate picture of the immune response dynamics.In this work, we aim to look at the stochastic nature of the immune response with specialfocus on the SARS-CoV-2 virus for the ongoing global COVID-19 pandemic. To do this, wemodel the immune response dynamics in viral infections by incorporating stochastic fluc-tuations in terms of the Gaussian white noise and the fractional Gaussian noise to the setof coupled ODEs of T cells and virus particles. Our stochastic immune response modelprovides a near analytic solution for the time dependent joint probability distribution of Tcells and virus, which we then use to determine the average number of virus particles in thesystem. Considering the ongoing pandemic, we then also compare our results to the experi-mental SARS-CoV-2 viral infection data from Germany [19, 20]. The long mean incubationperiod, which is approximated to be 5-6 days [21, 22] makes it important to analyse theviral dynamics of SARS-CoV-2 using the stochastic immune response model. Although, inthis work, we carry out the numerical analysis for this particular virus, similar analysis canalso be applied to other viral dynamics.This paper is organized as follows: In Section II, we formulate the stochastic immune re-sponse model using a set of coupled stochastic differential equations (SDEs). Section III pro-vides the derivation of the temporal evolution of the joint probability distribution (Fokker-Planck equation) of T cells and virus particles from the set of coupled SDEs. Part A ofthe results section includes a near analytic solution of Fokker-Planck equation in differenttime regimes. The analysis carried out here is generic and is applicable to viral infectionsin general. Part B has two subsections to show the numerical analysis of the stochasticquantities of immune response in SARS-CoV-2 infection. We finally summarize our resultsin Section V. II. STOCHASTIC IMMUNE RESPONSE MODEL
The simplest model of immune response dynamics needs to account for interactions be-tween the immune cells and the virus. The immune response itself is activated when cellsare under attack by foreign bodies such as viruses. This activation is manifested throughan increased proliferation of the immune cells. Therefore, there needs to be a coupled inter-action between the immune cells and the virus. Since the proliferation and death rates ofboth these entities are intrinsically stochastic [14, 23], any realistic model needs to accountfor this inherent variability as well. The inherent intrinsic variability can be accounted forin the model via the Gaussian white noise (GWN), which is delta correlated. The GWN is afast decaying noise, where successive fluctuations are not correlated. We model the inherentnoise in T cells as a GWN. Viruses though are known to show population level fluctuationsbetween active and latent states [13]. This is a manifestation of fluctuations in the geneexpression, which in turn leads to variability in the gene products. We therefore model theviral dynamics to evolve under the action of the fractional Gaussian noise (fGn), which is atype of colored noise. Specifically, the dynamical interaction between T cells and the virusparticles can then be written as a set of coupled stochastic differential equations (SDEs),given by ˙ x ( t ) = β (cid:0) v ( t ) (cid:1) x ( t ) − γx ( t ) + θ ( t ) (1) (cid:90) t dt (cid:48) K ( t − t (cid:48) ) ˙ v ( t (cid:48) ) = pv ( t ) − cv ( t ) + ξ ( t ) (2)Here, x ( t ) and v ( t ) represent the concentration levels of T cells and virus particles inthe body which are cleared at rates γ and c respectively. p is the replication rate of virusparticles. β (cid:0) v ( t ) (cid:1) is a function of virions concentration, which incorporates the effect of viruslevels in T cell dynamics. β (cid:0) v ( t ) (cid:1) is a positive odd integer power law function that accountsfor the dependency between T cells and virus levels. This term therefore couples the viraldynamics to T-cell dynamics. In general, one can consider β (cid:0) v ( t ) (cid:1) = rv m to account forthe fact that higher the value of m , faster is the rate of increase of β (cid:0) v ( t ) (cid:1) . Earlier immuneresponse models [6] also indicate the same, i.e., increase in virus levels lead to proliferationof T cells with a rate r . In our study, we have set the value of m to be 1. A schematicrepresentation of the stochastic immune response model is shown in Fig. 1.The terms θ ( t ) and ξ ( t ) represent GWN and fGn with noise strengths a and λ respectively.These account for the variability in the system and have the following statistical properties. (cid:104) θ ( t ) (cid:105) = 0 (cid:104) ξ ( t ) (cid:105) = 0 (cid:104) θ ( t ) θ ( t (cid:48) ) (cid:105) = aδ ( t − t (cid:48) ) (cid:104) ξ ( t ) ξ ( t (cid:48) ) (cid:105) = λK ( | t − t (cid:48) | ) (3)Here, the angular brackets represent an average over all realizations of the noise. Asspecified in Eq(3), both the noise terms, θ ( t ) and ξ ( t ), have zero mean. θ ( t ) represents r γ p c x Clearance of T cells
Clearance of virus
Virus
T cells
Figure 1.
Schematic representation of stochastic immune response model : r is the prolif-eration rate of T cells, p is the replication rate of virus particles, x is the production rate of basallevel T cells, c and γ are clearance rates of the virus particles and the T cells respectively. [Theseare representative images, not to scale] fluctuations of a Markov process which is delta correlated, i.e., the fluctuation at any time t is uncorrelated with the previous time t (cid:48) , while ξ ( t ) represents correlated fluctuationscharacteristic of fractional Gaussian noise (fGn). fGn is a Non-Markovian process, which istemporally correlated by a memory kernel K ( | t − t (cid:48) | ), which has the following form [24–28] K ( | t − t (cid:48) | ) = 2 H (2 H − | t − t (cid:48) | H − (4)Here H is the Hurst index such that 1 / ≤ H <
1. The value H = 1 / i.e. a and λ respectively,specifies the deviation of noise from its mean. III. TRANSFORMATION TO THE FOKKER-PLANCK EQUATION
The advantage of representing a system through a stochastic model is the possibility ofobtaining a multivariate probability distribution function from the set of SDEs, Eqs. (1)and (2). This distribution function can then be used to obtain the average values of therelevant variables, in this case, the number of virions inside the host’s body at a given time.To begin with, we first define the distribution function P ( x, v, t ) as the probability densityof finding x T cells and v virions at a particular time t . This can be written as P ( x, v, t ) = (cid:104) δ ( x − x ( t )) δ ( v − v ( t )) (cid:105) (5)where x ( t ) and v ( t ) are functionals of the noise θ ( t ) and ξ ( t ) respectively and the angularbrackets represent an average over all realizations of the noise. Now, substituting the solu-tions of the differential equations, Eqs(1) and (2) into Eq (5) and making use of the noiseproperties (Eq(3)), one can obtain the time evolution of the probability density function P ( x, v, t ). This equation, known as the Fokker-Planck equation, is given by ∂∂t P ( x, v, t ) = − β (cid:0) v ( t ) (cid:1) P ( x, v, t ) − L P ( x, v, t ) (6)where the operator L is defined as L = β (cid:0) v ( t ) (cid:1) x ∂∂x − γ ∂∂x x − a ∂ ∂x − η ( t ) ∂∂v v − λ | ( c − p ) | η ( t ) ∂ ∂v (7)Here, η ( t ) is a time-dependent function defined as η ( t ) ≡ − ˙ X ( t ) / X ( t ), where X ( t ) = E − H ( − ( t/τ ) − H ) and τ = (cid:16) Γ(2 H +1) | ( c − p ) | (cid:17) / (2 − H ) . E α,β ( z ) is a Mittag-Leffler function of theform (cid:80) ∞ n =0 z n / Γ( αn + β ), where Γ( αn + β ) is a gamma function.The details of this transformation (from Eqs(1) and (2) to Eq(6)) are shown in AppendixA. The exact joint probability distribution function, P ( x, v, t ), can be obtained from thetime dependent solution of the Fokker-Planck Equation (Eq(6)), which will then be usefulin studying the various statistics of T cells and virus particles. IV. RESULTSA. Solution of the Fokker-Planck Equation: Obtaining P ( x , v , t ) The Fokker-Planck Equation, Eq(6) can be represented through an equivalent form bymaking use of the Green’s function. The system evolves with time from its equilibrium stateand therefore the formal solution of Eq(6) is given by [29] P ( x, v, t ) = P eq ( x, v ) − (cid:90) ∞ dx (cid:48) (cid:90) ∞ dv (cid:48) (cid:90) t dt (cid:48) G ( x, v, t − t (cid:48) | x (cid:48) , v (cid:48) ) β ( v (cid:48) ) P ( x (cid:48) , v (cid:48) , t (cid:48) ) (8)Here the Green’s function, G ( x, v, t − t (cid:48) | x , v ), is the time dependent conditional probabilityof finding the system in the state ( x, v ) at time t given that it was in the state ( x , v ) attime t = 0. The detailed derivation of the Green’s function using the operator L is providedin Appendix B. After some lengthy algebra, the Green’s function (propagator) is given by G ( x, v, t | x , v ,
0) = 12 π (cid:113) a γ λ | ( c − p ) | (1 − e − γt )(1 − X ( t )) exp (cid:18) − (cid:18) a γ ( x − x e − γt ) (1 − e − γt )+ 1 λ | ( c − p ) | ( v − v X ( t )) (1 − X ( t )) (cid:19)(cid:19) (9)This propagator satisfies the condition that when t → ∞ , the time dependent conditionalprobability density , G ( x, v, t − t (cid:48) | x (cid:48) , v (cid:48) ) = P eq ( x, v ). Therefore, P eq ( x, v ) = 12 π (cid:113) a γ λ | ( c − p ) | exp (cid:32) − (cid:32) x ( a γ ) + v λ | ( c − p ) | (cid:33)(cid:33) (10)Eq(8) can provide the implicit solution for P ( x, v, t ), but to determine the analytic expressionfor the time dependent joint probability distribution function, we use the closure schemeintroduced by Wilemski and Fixman [30, 31]. The Wilemski-Fixman (WF) approximationwas originally developed to estimate the rate of the reaction between reactive groups ateither ends of a polymer chain [31]. This was later extended to the systems of catalyticbimolecular reactions [32], non-exponential DNA escape kinetics [33], dynamic disorder inchain unfolding [34] and chain closure in entangled polymer systems [35], among others.In the original work by Wilemski and Fixman [30], a sink term was introduced to providea valid “closure approximation” for the diffusion equation in the many-particle system ofpolymer reactions. Using this approximation, the solution of Eq(8) can be replaced byan approximate expression which involves the product of two terms - (i) the equilibriumprobability density P eq ( x, v ), that is time independent and corresponds to the situationwhen T cells and virus levels in the system are independent of each other and (ii) a self-consistently determined time-dependent term which evolves with time from the equilibriumdistribution ( P eq ( x, v )) as a consequence of the sink term, i.e., β (cid:0) v ( t ) (cid:1) in the case of ourstochastic immune response model.Therefore using this approximation, the probability distribution function P ( x, v, t ) canbe defined by the introduction of two functions w ( t ) and ¯ w such that, P ( x, v, t ) = P eq ( x, v ) w ( t )¯ w (11)where, w ( t ) = (cid:90) ∞ (cid:90) ∞ dxdvβ ( v ) P ( x, v, t ) and ¯ w = (cid:90) ∞ (cid:90) ∞ dxdvβ ( v ) P eq ( x, v ) (12)Multiplying Eq(8) by β ( v ), integrating over x and v and then substituting Eq(11) into itprovides the expression for w ( t ) such that w ( t ) = ¯ w − (cid:90) t dt (cid:48) C ( t − t (cid:48) ) w ( t (cid:48) ) / ¯ w (13)where C ( t − t (cid:48) ) = (cid:90) ∞ dx (cid:48) (cid:90) ∞ dx (cid:90) ∞ dv (cid:48) (cid:90) ∞ dvβ ( v ) G ( x, v, t − t (cid:48) | x (cid:48) , v (cid:48) ) β ( v (cid:48) ) P eq ( x (cid:48) , v (cid:48) ) (14)Following the above steps, one can see that w ( t ) is required in the calculation of the proba-bility distribution for the immune response model, which in turn requires the evaluation of C ( t ). Computation of C ( t ) is carried out by substituting Eq(9), Eq(10) and the expressionfor β ( v ) into Eq [14] and then carrying out the integration. This gives C ( t ) = (cid:16) π (cid:17) r λ | ( c − p ) | (cid:32) (cid:112) − X ( t ) + X ( t ) (cid:32) π + 2ArcTan (cid:32) X ( t ) (cid:112) − X ( t ) (cid:33)(cid:33)(cid:33)(cid:18) π + 2ArcTan (cid:18) e − γt √ − e − γt (cid:19)(cid:19) (15) w ( t ) can simply be obtained by using the method of Laplace transforms. This requires theLaplace transform of C ( t ) as well. Determining the simple algebraic form of the Laplacetransform of Eq(15) is non-trivial due to the presence of X ( t ) i.e. the Mittag-Leffler function.However, viral dynamics are particularly interesting during the early and late stages ofinfection. Viral populations reach their peaks in the early stage of infection and are especiallylong lived and decay gradually during the late stages of infection [20]. Therefore, here, weare primarily interested in viral dynamics, and in turn, X ( t ) in two different time regimesi.e. at short and long times.
1. Short time regime
In the case of short times ( t/τ << γt << X ( t ) = 1 − a t b + O ( t b ), where b = 2 − H , a = 1 /τ b Γ(3 − H ). Thefunctions ArcTan (cid:18) X ( t ) √ −X ( t ) (cid:19) ≈ π − (cid:0) − a t b (cid:1) √ a t b and ArcTan (cid:16) e − γt √ − e − γt (cid:17) ≈ π − (1 − γt ) √ γt , for small value of γ this approximates to π . Ignoring higher order terms of t , theexpression obtained for C ( t ) in the short time regime is, C ( t ) = (cid:18) π (cid:19) r ϑ (cid:0) − a t b (cid:1) (16)where ϑ = λ/ | ( c − p ) | Substituting the Laplace transform of Eq(16) into the Laplace Transform of Eq(13), i.e.,into w ( s ) = ¯ w/ (cid:16) s (cid:16) C ( s )¯ w (cid:17)(cid:17) , one can obtain the complete expression for w ( s ). This isgiven by w ( s ) = 4 ¯ ws b ( r π ϑs b / ¯ w ) + 4 s b − ( r π λ/ ( ¯ w Γ(3 − b ))) (17)The series expansion of the above expression gives w ( s ) = ¯ w (cid:88) ∞ k =0 ( − k (cid:18) r ϑπ w (cid:19) k s b ( k +1) ( s b − A ) k +1 (18)where A = ( rπ ) w λ Γ(3 − b ) . The inverse Laplace transform of Eq(18) provides the expression for w ( t ) in the short time regime. This is given by w ( t ) = ¯ w (cid:88) ∞ k =0 ( − k k ! (cid:18) r ϑπ w (cid:19) k t k E ( k )1+ b, − bk (cid:0) At b (cid:1) (19)where E kα,β ( z ) = d k E α,β ( z ) dz k is the k th derivative of the Mittag-Leffler function with respect toits argument. The series expansion of the Mittag-Leffler function [36] is given by E α,β ( z ) = (cid:80) ∞ n =0 z n Γ( αn + β ) whose k th derivative turns out to be E ( k )1+ b, − bk (cid:0) At b (cid:1) = (cid:88) ∞ n =0 ( n + k )! n ! ( At b ) n Γ( n (1 + b ) + k + 1) (20)Substitution of Eq(20) in Eq(19) results in the expression for w ( t ). For short times (i.e.,taking into account only the k = 0 term as in [36]), the final expression is w ( t ) = ¯ w E b, (cid:0) At b (cid:1) (21)where ¯ w = √ π r √ ϑ . Substituting the above expression into Eq(11) gives the joint proba-bility density expression, P ( x, v, t ), in the short time regime, which is P ( x, v, t ) = 12 π √(cid:61) ϑ exp (cid:18) − (cid:18) x (cid:61) + v ϑ (cid:19)(cid:19) E b, (cid:18) π / √ r Γ(3 − b ) √ λ (cid:112) | ( c − p ) | t b (cid:19) (22)where (cid:61) = a/ γ and ϑ = λ/ | ( c − p ) | . The consequences of the dynamics in the short timeregime will be discussed in Part B of this section.
2. Long time regime
Another regime of interest is the viral dynamics in the long time regime. Therefore,for large values of t , the Mittag-Leffler function is approximated to X ( t ) ≈ a t − b , where a = τ b / Γ(2 H −
1) and b = 2 − H [33]. Further, ignoring higher order terms and usingthe property that ArcTan( f ( t )) ≈ f ( t ) when f ( t ) <<
1, we can approximate the functionsArcTan (cid:18) X ( t ) √ −X ( t ) (cid:19) ≈ a t − b and ArcTan (cid:16) e − γt √ − e − γt (cid:17) ≈
0. After applying these approxima-tions, the closed form expression for C ( t ) in the long time regime is given by C ( t ) = (cid:16) π (cid:17) r ϑ (cid:0) π + π a t − b (cid:1) (23)0where ϑ = λ/ | ( c − p ) | Substituting the Laplace transform of Eq(23) into the Laplace Transform of Eq(13), i.e.,into w ( s ) = ¯ w/ (cid:16) s (cid:16) C ( s )¯ w (cid:17)(cid:17) , we get w ( s ) = 16 ¯ w r ϑπ (cid:16)
16 ¯ ws π r ϑ + π − b ) | ( c − p ) | s b (cid:17) (24)The series expansion of Eq(24) is w ( s ) = ¯ w (cid:88) ∞ k =0 ( − k k ! (cid:18) π r ϑ
16 ¯ w (cid:19) k k ! s − b ( k +1) ( s (1 − b ) + A ) k +1 (25)where A = (cid:16) rπ (cid:17) − b )¯ w λ | ( c − p ) | . The inverse Laplace transform of Eq(25) provides w ( t ) = ¯ w (cid:88) ∞ k =0 ( − k k ! (cid:18) π r ϑ
16 ¯ w (cid:19) k t k E k − b, bk (cid:0) − A t − b (cid:1) (26)Making use of the asymptotic form of the Mittag-Leffler function for larger values of itsargument i.e. E α,β ( z ) = z Γ( β − α ) , the k th derivative turns out to be E k − b, bk ( − A t − b ) = k !Γ( b (1 + k )) (cid:0) A t − b (cid:1) − ( k +1) (27)Upon substituting the above expression into Eq(24), the final expression for w ( t ) in thelong time regime is given by w ( t ) = 16 ¯ w a r ϑπ Γ(1 − b ) t − b E b,b (cid:18) − t b π | ( c − p ) | Γ(3 − b ) (cid:19) (28)where ¯ w = √ π r √ ϑ and E b,b is the Mittag-Leffler function. The joint probability distribu-tion function, P ( x, v, t ) in the long time regime is then given by P ( x, v, t ) = 8 π π ) / ra √(cid:61) ϑ Γ(1 − b ) t b − exp (cid:18) − (cid:18) x (cid:61) + v ϑ (cid:19)(cid:19) E b,b (cid:18) − t b π | ( c − p ) | Γ(3 − b ) (cid:19) (29)where (cid:61) = a/ γ and ϑ = λ/ | ( c − p ) | . The consequences of the dynamics in the long timeregime will be discussed in Part B of this section. Equations [22] and [29] are the mainresults of this work, which we further explore in the next section.1 B. Numerical implementations of stochastic immune response model
The analytical results presented in the previous section, i.e., Eqs(22) and (29), providethe joint probability density function P ( x, v, t ) for T cells and virus particles in different timeregimes, namely, the early (short times) and late (long times) stages of infection. Followingthis calculation, we carry out its numerical analysis to study the temporal evolution of thetime dependent joint probability distribution function. Our analysis applies in general totypical virus and immune response interactions, but, because the ongoing pandemic due toCOVID-19 demands special attention, we analyze our results in the context of the immuneresponse system for T cells and the SARS-CoV-2 virus, a pathogenic RNA virus with a lipidenvelope. SARS-CoV-2 has long mean incubation period, which is approximated to be 5-6days [21]. The latency period i.e. the time lag between the time of infection and the onsetof the initial symptoms, makes it pertinent to analyse the viral dynamics of SARS-CoV-2using the stochastic immune response model. To gain an insight into the applicability ofour model to the real world data, we have carried out a quantitative analysis by comparingour results to those obtained from the clinical data of COVID patients [19, 20].
1. Temporal evolution of the joint probability distribution function
The time dependent joint probability distribution expression for T cells and virus par-ticles in the short and the long time regimes are given by Eqs(22) and (29) respectively.In this work, we have studied the evolution of probability distribution function in two dif-ferent time regimes: (a) from the time when SARS-CoV-2 enters the host’s body to startthe infection and (b) during the extinction of SARS-CoV-2 from a patients’ body. Theparameters required for this study have been determined from the earlier research studies[8, 10] which have looked at the effect of SARS-CoV-2 on the basis of experimental dataof viral load within the patients’ bodies. The strength of the noise parameters, i.e., GWNin T cells and fGn in virus particles are a = 0 . λ = 0 . γ has been estimated on thebasis of the half life of T cells, which is approximated to be 4 - 34 days [37], therefore inthese calculations we set γ to 0 . − . Clearance rate of virus particles ( c ), proliferationrate of T cells ( r ), and replication rate of virus particles ( p ) are set to the mean values ofthese parameters for different patients’ in different time regimes. These are listed in TableI. Fig. 2 shows the temporal evolution of the distribution function in the initial period ofinfection, approximately up to 2 days. Fig. 3, on the other hand provides an insight intothe decrease in the joint probability distribution of T cells and virus particles at larger times2 t= 1.00 c) P ( x , v , t ) P ( x , v , t ) P ( x , v , t ) P ( x , v , t ) t= 0.80 b) t= 0.10 a) t= 1.40 c) Figure 2.
Temporal evolution of joint probability distribution of T cells and virionsduring the early infection period, (short time regime) for H = . a) P ( x, v ) at 0.1days. b) P ( x, v ) at 0.8 days. c) P ( x, v ) at 1 day. d) P ( x, v ) at 1.4 days from the start of infection. i.e. during virus extinction and when T cell levels drop to basal values. These results showthat the joint probability of T cells and virus particles increases with time at the start ofinfection and shows the decrease in probability distribution when virus levels attain a verylow value (virus extinction period) within the host’s body.
2. Average level of virus particles
The next step in the analysis of the stochastic immune response model is to compute theaverage level of virus particles in the system. The analytic expressions for this in differenttime regimes are obtained using Eq(22) and Eq(29). For short times, the average level ofvirus particles is given by, (cid:104) v ( t ) (cid:105) = √ ϑ √ π E b, (cid:18) π / √ r Γ(3 − b ) √ λ (cid:112) | ( c − p ) | t b (cid:19) (30)3 P ( x , v , t ) P ( x , v , t ) P ( x , v , t ) t= 10 a) t= 12 b) t= 15 c) t= 20 c) P ( x , v , t ) Figure 3.
Temporal evolution of joint probability distribution of T cells and virionsduring the late stage of infection, (the long time regime) for H = . a) P ( x, v ) after10 days. b) P ( x, v ) after 12 days. c) P ( x, v ) after 15 days. d) P ( x, v ) after 20 days of infection. where ϑ = λ/ | ( c − p ) | For long times, this is given by (cid:104) v ( t ) (cid:105) = 2 t b − ra π Γ(1 − b ) E b,b (cid:18) − t b π | ( c − p ) | Γ(3 − b ) (cid:19) ; (31)After having obtained the expression for mean levels of virus particles, the next step wouldbe to compare it to experimental data. Fig. 4 illustrates the comparison between thenumerically evaluated average number of virions with that determined clinically in COVIDpatients from Germany [20]. The rate parameters used in the analysis of each patient arelisted in Table I. The parameter values used have been selected on the basis of fits to dataavailable for experimental viral load in COVID patients in published papers [8, 10]. Otherparameters have been set to the same values as used in the analysis of the temporal evolutionof the joint probability distribution for T cells and virus particles, i.e. in Figs. 2 and 3. Inshort time regime, the initial value v is considered to be 100 copies/ml (which is the lowerlimit of detection in experiments [20]). In the long time regime, the initial value v is in4 Patients Regime c ( day − ) p ( day − ) r (ml/cells/day)A Short time 6.89 7.87 0.30Long time 7.50 4.20 0.33B Short time 5.32 6.81 0.33Long time 7.87 4.37 0.40C Short time 7.13 8.53 0.31Long time 7.67 5.19 0.39D Short time 4.93 5.54 0.34Long time 6.92 5.49 0.26E Short time 4.98 5.33 0.31Long time 6.92 5.31 0.26F Short time 5.27 7.90 0.11Long time 7.04 4.22 0.20G Short time 6.74 7.62 0.26Long time 11.71 6.23 0.27H Short time 6.22 9.73 0.35Long time 15.07 9.12 0.39Short time [4.93-7.13] [5.33-9.73] [0.11-0.35]Median 5.77 7.74 0.31Long time [6.92-15.07] [4.20-9.12] [0.20-0.40]Median 7.58 5.25 0.30Table I. Numerical values of rate parameters used in the quantitative analysis of the stochasticimmune response model for SARS-CoV-2 virus. c and p are respectively the clearance rate andthe replication rate of virus particles. r is the proliferation rate of T cells. range of 10 − copies/ml (assuming the peak viral load in patients). Eqs. (30) and (31)are then used to obtain the average virus levels at the early (short time regime) and latestages (long time regime during virus extinction i.e. after the viral peak inside the patient’sbody) of infection. Fig. 4 also shows that change in the value of Hurst index, determinesthe best fit to the experimental data of SARS-CoV-2 load in patients.The models for H values 0.55, 0.60 and 0.75 have been compared to the SARS-CoV-2viral load data [10] from patients. As listed in Table[II], the model with the best fit hasbeen determined by comparison of the mean square error (MSE) and the Akaike informationcriterion (AIC) for individual models. The AIC values for the individual models have been5 Figure 4.
Average level of SARS-CoV-2 virus using stochastic immune response model:
The average number of virions at different times is compared for H values 0 . , . , .
75. Numer-ical results from stochastic immune response model (solid lines) are compared to the experimentaldata (points) of viral load from German patients. Parameters are listed in Table I. Patients H = 0.55 H = 0.65 H = 0.75A AIC -8.51 4.72 12.93MSE 0.48 1.17 2.03B AIC -13.98 0.04 7.45MSE 0.28 0.83 1.47C AIC 6.36 7.44 13.15MSE 1.36 1.48 2.30D AIC -3.72 -3.57 -2.70MSE 0.49 0.50 0.55E AIC -7.37 -5.78 -2.52MSE 0.33 0.39 0.56F AIC -8.91 0.49 4.76MSE 0.28 0.79 1.27G AIC 24.06 19.52 18.37MSE 5.305 3.74 3.42H AIC 20.17 16.95 16.50MSE 4.40 3.36 3.24Table II. Comparison of AIC and MSE of three numerically different stochastic immune responsemodels to the experimental data. calculated using AIC = n log (cid:18) RSSn (cid:19) + 2 mnn − m − n is the number of data points, m is the number of unknown parameters and RSS is the residual sum of squares obtained from the fitting routine [38, 39]. Lower the valueof MSE and AIC, better is the model fit to the experimental data of viral load in patients.Fig. 4 shows that the maximum viral load or viral peak lies between the predicted resultsfor short time and long time regimes.
V. SUMMARY AND CONCLUSION
Understanding the functioning and dynamics of the immune system becomes importantgiven the role that it plays in fighting off infection and disease. However, this becomesnon-trivial because just like all the other biological processes at the cellular level, this tooshows a lot of variability. Therefore, stochastic models of the immune response, most ofwhich primarily focused on the HIV virus, have been more successful in explaining some7of the heterogeneity and variability associated with the system [14, 15, 17, 18, 40]. Earliertheoretical studies that have attempted to explain different aspects of the immune responsehave been based on deterministic models [5–8, 38], which illustrate the mean dynamics butfail to account for the inherent stochasticity and variability of the process. In light of this,in this article, we have developed and analyzed a stochastic version of the immune responsemodel.Our stochastic immune response model is composed of coupled Langevin equations for theT cells and virus particles with two kinds of noise, GWN and fGN, respectively. Thisallowed us to account for stochasticity within the model itself and obtain the temporalprobability distributions of the main variables. In this work, we first derived the Fokker-Planck Equation, which we then used to compute the joint probability distribution of T cellsand virus particles by making use of the Wilemski-Fixman approximation. This approachallowed us to obtain analytical solutions of the probability distribution functions and theaverage virus particles in the limit of short and long times, showing how the infection beginsand ends (see Figs. 2 and 3).A further advantage of an analytical expression is that a direct comparison can be madebetween the predicted theoretical dynamics and the experimental results. We have carriedout such a comparison with the available SARS-CoV-2 virus data from patients in Germany.At short times, i.e., during the early period of infection, the model predicts that there isa steep rise in the virus levels with time, whereas, at long times, the virus levels dropgradually, in accordance with the model’s prediction. As shown in Fig. 4, our StochasticImmune Response model gives a good fit to the experimental data at both short and longtimes.One of the parameters that is crucial in obtaining good fits is the Hurst index, H , whichin the case of fGn takes values between 1/2 and 1. H value between 0.5 and 1 correspondsto a system with long-ranged correlated fluctuations and values between 0 and 0.5 standfor anti-correlated time series [12]. The H value in the case of fGn, that represents theNon-Markovian viral dynamics, indicates the long-ranged time correlation of the noise ξ ( t ).Higher the H value, greater is the correlation between noises at any time t and the previoustime t (cid:48) . The H = 1 / ξ ( t )is completely uncorrelated to previous times t (cid:48) and therefore the two successive times aredelta correlated. In addition to the expression of the average level of virions at arbitrary H ,we have also calculated its expression in the limit of H = 1 /
2. In this limit, the expressionsimplifies to simple exponentials. In the short time regime (cid:104) v ( t ) (cid:105) = √ ϑ √ π cosh (cid:18)(cid:18) π / ra √ ϑ √ (cid:19) / t (cid:19) (33)8and in the long time regime it is given by (cid:104) v ( t ) (cid:105) = 2 rπ | ( c − p ) | Γ(3 − b ) exp (cid:18) − t | ( c − p ) | π (cid:19) (34)In the long time regime, for H = 1 /
2, as evident from Eq(34), there is fast exponentialdecrease in the virus levels. However, as seen from Fig. 4, the long time regime shows slowtemporal decay. This cannot be explained by the GWN limit (Eq(34)) of the average viruslevels. Therefore, our model, which includes long ranged noise correlations through fGnprovides a more accurate picture of the viral dynamics. For most of the plots in Fig. 4, H = 0 .
55 gives a better fit in comparison to other values.We have also looked at the effect of the strength of the noise on T cell dynamics. The averagelevel of T cells in two different time regimes can be derived using Eq(22) and Eq(29). In theshort time regime, this is given by (cid:104) x ( t ) (cid:105) = 14 √ π (cid:114) aγ E b, (cid:18) π / √ r Γ(3 − b ) f ( λ ) (cid:112) | ( c − p ) | t b (cid:19) (35)where f ( λ ) = √ λ . At long times, the average level of T cells is given by (cid:104) x ( t ) (cid:105) = √ rπ | ( c − p ) | Γ(3 − b ) (cid:114) aγ t b − E b,b (cid:18) − t b π | ( c − p ) | Γ(3 − b ) (cid:19) ; (36)From Eq(35), it is clear that the terms (cid:112) a/γ and f ( λ ) account for the strength of thenoise in T cell and viral dynamics respectively. The presence of f ( λ ) in the argument ofthe Mittag-Leffler function leads to a faster increase in the T cells level with increased λ (as clearance rate γ of T cells has been fixed). Thus in the short time regime, increase inthe strength of the noise will cause an increase in the rate of production of T cells. Thisphenomenon may affect the system in a way where T cells attain its peak value before themaximum viral load and thus might not be optimized to clear all the virus.In the long time regime, the expression for the average levels of T cells is a product ofa slowly increasing function (power law in time) and a decreasing Mittag-Leffler function.From the expression in Eq(36), it is clear that it is only the pre-factor (cid:112) a/γ that accountsfor the effect of the strength of noise. Thus, an increase in the value of a will affect theincreasing function, but the net effect on T cell levels will not be significant. The T cellsdynamics will show long lived dynamics in the long time regime and will have similar valuesfor different strengths of the noise. Thus, our model predicts that the noise in the systemmay have a major effect at the start of the infection time. The time at which the populationof T cells reaches the maximum value within patients is an important factor. The analysisof these levels may then help in determining when the immune response modifiers should beadministered to the patients.9The present formulation can also be extended to incorporate increasing complexity by con-sidering the effects of susceptible and infected cells on the immune response system. Thismodel can provide useful insights into the dynamics of various other viral diseases as well,such as measles, influenza, Zika virus, which also have long incubation period as found inSARS-CoV-2 [41, 42]. The colored noise incorporated in the model accounts for the longlived dynamics of the virus and can therefore provide more accurate predictions. Thesestochastic models can therefore help in a better understanding of the immune responsesystem. DATA AVAILABILITY
The data that supports the findings of this study are available within the article.
ACKNOWLEDGEMENTS
This work is supported by the Science and Engineering Research Board (SERB) MATRICSGrant (Ref. No. MSC/2020/000370) awarded by Department of Science and Technology(DST), India.
Appendix A: Derivation of the Fokker-Planck Equation
We have carried out the derivation of the Fokker-Planck equation by using the methodsdescribed in [43]The multivariate probability density distribution of x and v at time t is given by P ( x, v, t ) = (cid:104) δ ( x − x ( t )) δ ( v − v ( t )) (cid:105) (A1)where x ( t ) and v ( t ) are functionals of θ ( t ) and ξ ( t ) respectively. Differentiation of Eq(A1)with respect to time t gives ∂∂t P ( x, v, t ) = − ∂∂x (cid:28) δ ( x − x ( t )) δ ( v − v ( t )) ˙ x ( t ) (cid:29) − ∂∂v (cid:28) δ ( x − x ( t )) δ ( v − v ( t )) ˙ v ( t ) (cid:29) (A2)Solution of Eq(A2) is obtained as follows. Laplace transform of Eq(2) provides ˙ v ( t ) suchthat v ( t ) = v (0) X ( t ) + 1 | ( c − p ) | (cid:90) t dt (cid:48) ξ ( t (cid:48) ) φ ( t − t (cid:48) ) (A3)where X ( t ) and φ ( t ) are inverse Laplace transforms of˜ X ( s ) = ˜ K ( s ) | ( c − p ) | + s ˜ K ( s ) and ˜Φ( s ) = 1 − s ˜ X ( s ) (A4)0respectively. By making use of the definition X (0) = 1 and eliminating v (0) in Eq(A3), weget, ˙ v ( t ) = ˙ X ( t ) X ( t ) v ( t ) + 1 | ( c − p ) | X ( t ) ddt (cid:18)(cid:90) t dt (cid:48) φ ( t − t (cid:48) ) ξ ( t (cid:48) ) X ( t ) (cid:19) (A5)Substituting ˙ x ( t ) from Eq(1)into Eq(A2) and taking an average over all realizations of thenoise, we obtain, ∂∂t P ( x, v, t ) = (cid:18) − β (cid:0) v ( t ) (cid:1) ∂∂x x + γ ∂∂x x + 12 a ∂ ∂x (cid:19) P ( x, v, t ) − ∂∂v (cid:104) δ ( x − x ( t )) δ ( v − v ( t )) ˙ v ( t ) (cid:105) (A6)Substitution of Eq(A5) into Eq(A6) gives ∂∂t P ( x, v, t ) = (cid:18) − β (cid:0) v ( t ) (cid:1) ∂∂x x + γ ∂∂x x + 12 a ∂ ∂x (cid:19) P ( x, v, t ) + η ( t ) ∂∂v vP ( x, v, t ) − | ( c − p ) | ∂∂v (cid:28) δ ( x − x ( t )) δ ( v − v ( t )) ξ ( t ) (cid:29) (A7)where, η ( t ) ≡ − ˙ X ( t ) X ( t ) and ξ ( t ) ≡ X ( t ) ddt (cid:18)(cid:90) t dt (cid:48) φ ( t − t (cid:48) ) ξ ( t (cid:48) ) X ( t ) (cid:19) (A8) ξ ( t ) is linearly related to ξ ( t ), which is a Gaussian random function. Therefore, by Novikov’stheorem [44], we get, (cid:10) δ ( x − x ( t )) δ ( v − v ( t )) ξ ( t ) (cid:11) = (cid:90) t dt (cid:48) (cid:10) ξ ( t ) ξ ( t (cid:48) ) (cid:11)(cid:28) δδξ ( t (cid:48) ) δ ( x − x ( t )) δ ( v − v ( t )) (cid:29) = − ∂∂v (cid:90) t dt (cid:48) (cid:10) ξ ( t ) ξ ( t (cid:48) ) (cid:11) × (cid:28) δ ( x − x ( t )) δ ( v − v ( t )) δv ( t ) δξ ( t (cid:48) ) (cid:29) (A9)To find the functional derivative in Eq(A9), Eq(A5) is solved by using the integrating factor,such that v ( t ) = exp (cid:18) − (cid:90) t dt (cid:48) η ( t (cid:48) ) (cid:19) (cid:34) v + 1 | ( c − p ) | (cid:90) t dt (cid:48) exp (cid:32)(cid:90) t (cid:48) dt (cid:48)(cid:48) η ( t (cid:48)(cid:48) ) (cid:33) ξ ( t (cid:48) ) (cid:35) (A10)Here the term v exp (cid:16) − (cid:82) t dt (cid:48) η ( t (cid:48) ) (cid:17) is a complementary function which is obtained by satis-fying the initial conditions, and the other term is the particular integral. Thus the functionalderivative is given by δv ( t ) δξ ( t (cid:48) ) = 1 | ( c − p ) | exp (cid:18) − (cid:90) tt (cid:48) dt η ( t ) (cid:19) (A11)1Therefore, (cid:10) δ ( x − x ( t )) δ ( v − v ( t )) ξ ( t ) (cid:11) = − ∂∂v | ( c − p ) | P ( x, v, t ) D ( t ) (A12)where D ( t ) = (cid:90) t dt (cid:48) ξ ( t ) ξ ( t (cid:48) ) exp (cid:18) − (cid:90) tt (cid:48) dt η ( t ) (cid:19) (A13)Substitution of Eq(A12) and Eq(A13) in Eq(A7), gives the expression, ∂∂t P ( x, v, t ) = − β (cid:0) v ( t ) (cid:1) ∂∂x xP ( x, v, t ) + γ ∂∂x xP ( x, v, t ) + 12 a ∂ ∂x P ( x, v, t )+ η ( t ) ∂∂v vP ( x, v, t ) + 1 | ( c − p ) | ∂ ∂v P ( x, v, t ) D ( t ) (A14)Substituting ξ ( t ) from Eq(A8) into Eq(A13), we get D ( t ) = 12 X ( t ) X ( t (cid:48) ) ddt X ( t ) X ( t (cid:48) ) (cid:90) t (cid:48) dt (cid:90) t dt φ ( t − t ) φ ( t − t ) ξ ( t ) ξ ( t ) (A15)Solution of Eq [A15] is obtained by making use of double Laplace transforms and performinga lengthy algebra as mentioned in [45], which gives D ( t ) = 12 λ | ( c − p ) | X ( t ) ddt X ( t ) (cid:18) − (cid:0) X ( t ) (cid:1) (cid:19) (A16)Differentiation of Eq(A16) and substitution of Eq(A8) in Eq(A14) gives us the desiredFokker-Planck equation. ∂∂t P ( x, v, t ) = (cid:18) − β (cid:0) v ( t ) (cid:1) − β (cid:0) v ( t ) (cid:1) x ∂∂x + γ ∂∂x x + 12 a ∂ ∂x + η ( t ) ∂∂v v + λ | ( c − p ) | η ( t ) ∂ ∂v (cid:19) P ( x, v, t ) (A17) Appendix B: Derivation of the propagator, G ( x, v, t | x , v , The Green’s function follows the equation ∂∂t G ( x, v, t | x , v ,
0) = − L G ( x, v, t | x , v ,
0) (B1)with the initial condition given by G ( x, v, | x , v ) = δ ( x − x ) δ ( v − v ) (B2)Hee, the operator L (from Eq(7)) is given by − L = − β ( v ) x ∂∂x + γ ∂∂x x + 12 a ∂ ∂x + η ( t ) ∂∂v v + λ | ( c − p ) | η ( t ) ∂ ∂v (B3)2where β ( v ) = rv m ; m = 1. Green’s function can be found explicitly by using the methodof Fourier transforms, where,ˆ G ( k , k , t | x , v ,
0) = 12 π (cid:90) ∞−∞ dk (cid:90) ∞−∞ dk exp( ιk x ) exp( ιk v ) G ( x, v, t | x , v ,
0) (B4)Carrying out the Fourier transform of Eq [B1], we get, ∂∂t ˆ G ( k , k , t | x , v ,
0) = (cid:18) − r (cid:18) ι ∂∂k (cid:19) (cid:18) ι ∂∂k (cid:19) ( ιk ) + γ + γ (cid:18) ι ∂∂k (cid:19) ( ιk ) + 12 a ( ιk ) + η ( t ) + η ( t ) (cid:18) ι ∂∂k (cid:19) ( ιk ) + λ | ( c − p ) | η ( t )( ιk ) (cid:19) ˆ G ( k , k , t | x , v , G ( k , k , t | x , v ,
0) gives ∂∂t ln ˆ G ( k , k , t | x , v ,
0) = ιrk ∂∂k ∂∂k ln ˆ G ( k , k , t | x , v ) + γ − γk ∂∂k ln ˆ G ( k , k , t | x , v , − ak + η ( t ) − η ( t ) k ∂∂k ln ˆ G ( k , k , t | x , v , − λ | ( c − p ) | η ( t ) k (B6)Following the Gaussian ansatz,ln ˆ G ( k , k , t | x , v ,
0) = ιk m ( t ) + ιk b ( t ) − k s ( t ) − k y ( t ) − k k z ( t ) (B7)Differentiating Eq(B7) with time and equating terms of corresponding powers of k , k , k , k and k k with those of Eq(B6), we get, ddt m ( t ) = − γm ( t ) − rz ( t ) ddt b ( t ) = − η ( t ) b ( t ) ddt s ( t ) = − γs ( t ) + addt y ( t ) = − η ( t ) (cid:18) y ( t ) − λ | ( c − p ) | (cid:19) ddt z ( t ) = − γz ( t ) + η ( t ) z ( t ) (B8)Using the conditions m (0) = x , b (0) = v , s (0) = 0, y (0) = 0 and z (0) = 0, the solution ofdifferential equations in Eq(B8) turns out to be m ( t ) = x e − γt ; b ( t ) = X ( t ) v ; s ( t ) = a γ (1 − e − γt ); y ( t ) = λ | ( c − p ) | (cid:18) − X ( t ) (cid:19) ; z ( t ) = 0 (B9)3The double inverse Fourier transform of Eq(B7) gives the expression for G ( x, v, t | x , v , G ( x, v, t | x , v ,
0) = 12 π (cid:113) y ( t ) s ( t ) − z ( t ) exp (cid:18) − (cid:16) − z ( t ) s ( t ) y ( t ) (cid:17)(cid:20) ( x − m ( t )) s ( t ) + ( v − b ( t )) y ( t ) − z ( t )( x − m ( t ))( v − b ( t )) s ( t ) y ( t ) (cid:21) (cid:19) (B10)Eq(B10) is in the form of a bi-variate Gaussian distribution for two random variables. Sub-stitution of Eq(B9) in Eq(B10), provides the expression for the Green’s function mentionedin Eq(9), which is, G ( x, v, t | x , v ,
0) = 12 π (cid:113) a γ λ | ( c − p ) | (1 − e − γt )(1 − X ( t )) exp (cid:18) − (cid:18) a γ ( x − x e − γt ) (1 − e − γt )+ 1 λ | ( c − p ) | ( v − v X ( t )) (1 − X ( t )) (cid:19)(cid:19) (B11) [1] Jacqueline Parkin and Bryony Cohen. An overview of the immune system. The Lancet ,357(9270):1777–1789, 2001.[2] Alan S. Perelson and G´erard Weisbuch. Immunology for physicists.
Reviews of ModernPhysics , 69(4):1219–1267, 1997.[3] David D. Chaplin. Overview of the immune response.
Journal of Allergy and Clinical Im-munology , 125(2, Supplement 2):S3–S23, 2010. Primer on Allergic and Immunologic Diseases.[4] Jean S. Marshall, Richard Warrington, Wade Watson, and Harold L. Kim. An introduction toimmunology and immunopathology.
Allergy, Asthma and Clinical Immunology , 14(s2):1–10,2018.[5] Alan S. Perelson. Modelling viral and immune system dynamics.
Nature Reviews Immunology ,2(1):28–36, 2002.[6] Alessandro Boianelli, Van Kinh Nguyen, Thomas Ebensen, Kai Schulze, Esther Wilk, NiharikaSharma, Sabine Stegemann-Koniszewski, Dunja Bruder, Franklin R. Toapanta, Carlos A.Guzm´an, Michael Meyer-Hermann, and Esteban A. Hernandez-Vargas.
Modeling influenzavirus infection: A roadmap for influenza research , volume 7. 2015.[7] Katharine Best and Alan S. Perelson.
Mathematical modeling of within-host Zika virus dy-namics , volume 285. 2018.[8] Esteban A Hernandez-vargas and Jorge X Velasco-hernandez. In-host Modelling of COVID-19in Humans. pages 1–19, 2020. [9] Prasith Baccam, Catherine Beauchemin, Catherine A. Macken, Frederick G. Hayden, andAlan S. Perelson. Kinetics of influenza a virus infection in humans. Journal of Virology ,80(15):7590–7599, 2006.[10] Sunpeng Wang, Yang Pan, Quanyi Wang, Hongyu Miao, Ashley N. Brown, and Libin Rong.Modeling the viral dynamics of SARS-CoV-2 infection.
Mathematical Biosciences , 328(Au-gust):108438, 2020.[11] Alexis Erich S. Almocera, Van Kinh Nguyen, and Esteban A. Hernandez-Vargas. Multiscalemodel within-host and between-host for viral infectious diseases.
Journal of MathematicalBiology , 77(4):1035–1057, 2018.[12] Ester L´azaro, Cristina Escarm´ıs, Juan P´erez-Mercader, Susanna C. Manrubia, and EstebanDomingo. Resistance of virus to extinction on bottleneck passages: Study of a decayingand fluctuating pattern of fitness loss.
Proceedings of the National Academy of Sciences ,100(19):10830–10835, 2003.[13] Abhyudai Singh, Brandon Razooky, Chris D. Cox, Michael L. Simpson, and Leor S. Wein-berger. Transcriptional bursting from the HIV-1 promoter is a significant source of stochasticnoise in HIV-1 gene expression.
Biophysical Journal , 98(8):L32–L34, 2010.[14] E. D. Hawkins, M. L. Turner, M. R. Dowling, C. Van Gend, and P. D. Hodgkin. A model ofimmune regulation as a consequence of randomized lymphocyte division and death times.
Pro-ceedings of the National Academy of Sciences of the United States of America , 104(12):5032–5037, 2007.[15] Nirav Dalal, David Greenhalgh, and Xuerong Mao. A stochastic model for internal HIVdynamics.
Journal of Mathematical Analysis and Applications , 341(2):1084–1101, 2008.[16] Xiying Wang, Yuanxiao Li, and Xiaomei Wang. The Stochastic Stability of Internal HIVModels with Gaussian White Noise and Gaussian Colored Noise.
Discrete Dynamics in Natureand Society , 2019, 2019.[17] Farzad Fatehi, Sergey N. Kyrychko, Aleksandra Ross, Yuliya N. Kyrychko, and Konstantin B.Blyuss. Stochastic effects in autoimmune dynamics.
Frontiers in Physiology , 9(FEB):1–14,2018.[18] Susmita Roy and Biman Bagchi. Fluctuation theory of immune response: A statistical me-chanical approach to understand pathogen induced T-cell population dynamics.
Journal ofChemical Physics , 153(4), 2020.[19] Merle M B¨ohmer et al.
Investigation of a covid-19 outbreak in germany resulting from a singletravel-associated primary case: a case series.
The Lancet Infectious Diseases , 20(8):920–928,2020. [20] Roman et al. , W¨olfel. Virological assessment of hospitalized patients with covid-2019. Nature ,581(7809):465–469, May 2020.[21] Roy M Anderson, Hans Heesterbeek, Don Klinkenberg, and T D´eirdre Hollingsworth. Howwill country-based mitigation measures influence the course of the covid-19 epidemic?
TheLancet , 395(10228):931–934, 2020.[22] The incubation period of coronavirus disease 2019 (covid-19) from publicly reported confirmedcases: Estimation and application.
Annals of Internal Medicine , 172(9):577–582, 2020. PMID:32150748.[23] Sudha B. Singh, Wojciech Ornatowski, Isabelle Vergne, John Naylor, Monica Delgado, EstebanRoberts, Marisa Ponpuak, Sharon Master, Manohar Pilli, Eileen White, Masaaki Komatsu,and Vojo Deretic. Human IRGM regulates autophagy and cell-autonomous immunity functionsthrough mitochondria.
Nature Cell Biology , 12(12):1154–1165, 2010.[24] S. C. Kou and X. Sunney Xie. Generalized langevin equation with fractional gaussian noise:Subdiffusion within a single protein molecule.
Physical Review Letters , 93(18):1–4, 2004.[25] Wei Min, Guobin Luo, Binny J. Cherayil, S. C. Kou, and X. Sunney Xie. Observation of apower-law memory Kernel for fluctuations within a single protein molecule.
Physical ReviewLetters , 94(19):1–4, 2005.[26] Benoit B. Mandelbrot and John W. Van Ness. Fractional brownian motions, fractional noisesand applications.
SIAM Review , 10(4):422–437, 1968.[27] S. C. Lim and S. V. Muniandy. Self-similar Gaussian processes for modeling anomalous diffu-sion.
Physical Review E - Statistical Physics, Plasmas, Fluids, and Related InterdisciplinaryTopics , 66(2), 2002.[28] Kwok Sau Fa and E. K. Lenzi. Time-fractional diffusion equation with time dependent diffu-sion coefficient.
Physical Review E - Statistical, Nonlinear, and Soft Matter Physics , 72(1):0–3,2005.[29] Gerald Wilemski and Marshall Fixman. General theory of diffusion-controlled reactions.4009(1973), 2001.[30] Gerald Wilemski and Marshall Fixman. Diffusion-controlled intrachain reactions of polymers.I Theory.
The Journal of Chemical Physics , 60(3):866–877, 1974.[31] Gerald Wilemski and Marshall Fixman. Diffusion-controlled intrachain reactions of polymers.II Results for a pair of terminal reactive groups.
The Journal of Chemical Physics , 60(3):878–890, 1974.[32] O. B´enichou, M. Coppey, M. Moreau, and G. Oshanin. Kinetics of diffusion-limited cat-alytically activated reactions: An extension of the Wilemski-Fixman approach.
Journal ofChemical Physics , 123(19), 2005. [33] Debarati Chatterjee and Binny J. Cherayil. Anomalous reaction-diffusion as a model of non-exponential DNA escape kinetics. Journal of Chemical Physics , 132(2), 2010.[34] Debarati Chatterjee and Binny J. Cherayil. The stretching of single poly-ubiquitin molecules:Static versus dynamic disorder in the non-exponential kinetics of chain unfolding.
Journal ofChemical Physics , 134(16):1–6, 2011.[35] Pinaki Bhattacharyya, Rati Sharma, and Binny J. Cherayil. Confinement and viscoelasticeffects on chain closure dynamics.
Journal of Chemical Physics , 136(23), 2012.[36] A. D. Vi˜nales and M. A. Desp´osito. Anomalous diffusion: Exact solution of the general-ized Langevin equation for harmonically bounded particle.
Physical Review E - Statistical,Nonlinear, and Soft Matter Physics , 73(1):5–8, 2006.[37] M. McDonagh and E. B. Bell. The survival and turnover of mature and immature cd8 t cells.
Immunology , 84(4):514–520, Apr 1995. 7790023.[38] Kasia A. Pawelek, Giao T. Huynh, Michelle Quinlivan, Ann Cullinane, Libin Rong, andAlan S. Perelson. Modeling within-host dynamics of influenza virus infection including immuneresponses.
PLoS Computational Biology , 8(6), 2012.[39] Kenneth P. Burnham and David R. Anderson. Multimodel inference: Understanding AIC andBIC in model selection.
Sociological Methods and Research , 33(2):261–304, 2004.[40] Dao Guang Wang, Shaobing Wang, Bo Huang, and Feng Liu. Roles of cellular heterogeneity,intrinsic and extrinsic noise in variability of p53 oscillation.
Scientific Reports , 9(1):1–11, 2019.[41] Justin Lessler, Nicholas Reich, Ron Brookmeyer, Trish Perl, Kenrad Nelson, and Derek Cum-mings. Incubation periods of acute respiratory viral infections: A systematic review.
TheLancet infectious diseases , 9:291–300, 06 2009.[42] Elisabeth Krow-Lucal, Brad Biggerstaff, and Jessica Staples. Estimated incubation period forzika virus disease.
Emerging Infectious Diseases , 23, 05 2017.[43] Srabanti Chaudhury and Binny J. Cherayil. Complex chemical kinetics in single enzymemolecules: Kramers’s model with fractional Gaussian noise.
Journal of Chemical Physics ,125(2), 2006.[44] E A Novikov. Novikov1965Functionals. 20(5):1290–1294, 1965.[45] Ronald Forrest Fox. The generalized Langevin equation with Gaussian fluctuations.