A deterministic model is proposed to describe the interaction between an immune system and an invading virus whose target cells circulate in the blood. The model is a system of two ordinary first order quadratic delay-differential equations with stipulated initial conditions, whose coefficients are eventually constant, so that the system becomes autonomous. The long-term behavior of the solution is investigated with some success. In particular, we find two simple functions of the parameters of the model, whose signs often, but not always, determine whether the virus persists above a nonzero threshold in the circulation or heads toward extinction.
Citation: Joseph E. Carroll. A two-dimensional discrete delay-differential system model of viremia[J]. Mathematical Biosciences and Engineering, 2022, 19(11): 11195-11216. doi: 10.3934/mbe.2022522
[1] | Shaoli Wang, Jianhong Wu, Libin Rong . A note on the global properties of an age-structured viral dynamic model with multiple target cell populations. Mathematical Biosciences and Engineering, 2017, 14(3): 805-820. doi: 10.3934/mbe.2017044 |
[2] | Mohammed Meziane, Ali Moussaoui, Vitaly Volpert . On a two-strain epidemic model involving delay equations. Mathematical Biosciences and Engineering, 2023, 20(12): 20683-20711. doi: 10.3934/mbe.2023915 |
[3] | Xin-You Meng, Yu-Qian Wu . Bifurcation analysis in a singular Beddington-DeAngelis predator-prey model with two delays and nonlinear predator harvesting. Mathematical Biosciences and Engineering, 2019, 16(4): 2668-2696. doi: 10.3934/mbe.2019133 |
[4] | Jiawei Deng, Ping Jiang, Hongying Shu . Viral infection dynamics with mitosis, intracellular delays and immune response. Mathematical Biosciences and Engineering, 2023, 20(2): 2937-2963. doi: 10.3934/mbe.2023139 |
[5] | Tyler Cassidy, Morgan Craig, Antony R. Humphries . Equivalences between age structured models and state dependent distributed delay differential equations. Mathematical Biosciences and Engineering, 2019, 16(5): 5419-5450. doi: 10.3934/mbe.2019270 |
[6] | Bao-Zhu Guo, Li-Ming Cai . A note for the global stability of a delay differential equation of hepatitis B virus infection. Mathematical Biosciences and Engineering, 2011, 8(3): 689-694. doi: 10.3934/mbe.2011.8.689 |
[7] | Sukhitha W. Vidurupola, Linda J. S. Allen . Basic stochastic models for viral infection within a host. Mathematical Biosciences and Engineering, 2012, 9(4): 915-935. doi: 10.3934/mbe.2012.9.915 |
[8] | Ning Bai, Rui Xu . Mathematical analysis of an HIV model with latent reservoir, delayed CTL immune response and immune impairment. Mathematical Biosciences and Engineering, 2021, 18(2): 1689-1707. doi: 10.3934/mbe.2021087 |
[9] | Edward J. Allen . Derivation and computation of discrete-delayand continuous-delay SDEs in mathematical biology. Mathematical Biosciences and Engineering, 2014, 11(3): 403-425. doi: 10.3934/mbe.2014.11.403 |
[10] | Heping Ma, Hui Jian, Yu Shi . A sufficient maximum principle for backward stochastic systems with mixed delays. Mathematical Biosciences and Engineering, 2023, 20(12): 21211-21228. doi: 10.3934/mbe.2023938 |
A deterministic model is proposed to describe the interaction between an immune system and an invading virus whose target cells circulate in the blood. The model is a system of two ordinary first order quadratic delay-differential equations with stipulated initial conditions, whose coefficients are eventually constant, so that the system becomes autonomous. The long-term behavior of the solution is investigated with some success. In particular, we find two simple functions of the parameters of the model, whose signs often, but not always, determine whether the virus persists above a nonzero threshold in the circulation or heads toward extinction.
Over the years, many mathematical models describing the interaction between the human immune system and invading organisms have been proposed. Although there has been interest in nonspecific viruses [1], a great number of these models, probably most, have focused on the dynamics of human HIV infection. As expected, COVID-19 models are now multiplying [2, 3].
Many of the models previously advanced have been systems of first order nonlinear ordinary differential equations [4, 5, 6, 7, 8, 9, 10, 11, 12]. However, one of the key features of viral dynamics is a lag between target cell infection and subsequent lysis with expulsion of virus into the circulation. To depict this, other models employing delay-differential equations have been proposed [3, 13, 14, 15, 16], even with distributed delays [17], usually to mimic this lag but, in one case [18], for the time gap before immune activation. In [19], a delay effect was achieved by adding the time since infection as a new parameter of infected cells and appending a partial differential equation to the system describing them further. In [2] a distributed delay was handled by use of the so-called linear chain trick, which interprets the delay as a sum of exponential distributions (equivalently, an Erlang distribution). These methods allowed sidestepping the problem in stability analysis of having to consider the infinitely many roots of transcendental polynomials [20, 21].
The most general of these differential and delay-differential models have not been autonomous, but analysis has commonly been done on autonomous simplifications. Often, both the number of parameters and dimension of these systems have been substantial, permitting a very granular description of the biologic process. Some [10, 13, 15, 19] have added in pharmacological treatment and the response to it.
The present model is simpler than most. It describes the natural history of viral infection without intervention in a host animal. It is a very basic two-dimensional delay-differential system, largely generalizable to a model with distributed delays which is messier and not included here. It assumes that interactions take place in a single, well-mixed compartment, e.g., the circulation, and as such, is limited in application to the human viral infections about which there is currently most concern including, besides HIV, hepatitis B and C, as well as COVID-19. The simplicity leads to an easier investigation in some ways but, as noted above, delay-differential equations tend to be much less tractable than differential equations.
The first aim for these kinds of models is to delineate the circumstances, in terms of their parameter values, which guarantee extinction versus persistence of both the virus and infected cells. For this model, sufficient conditions for both extinction and for persistence above a positive level are found, but neither is shown to be necessary. There is an "interval" of parameter states for which either of these outcomes might be possible, or even persistence with extinction as a limit point, an unlikely scenario physiologically. Nevertheless, the results may point to which parameters and what particular functions of their values are important in real life and, secondarily, which ones are the best to aim at tweaking with pharmacological treatments, including vaccines, in order to provide the most effective effort.
From the mathematical perspective, many of the arguments employed here are atypical for the analysis of delay-differential systems and might find use for other problems described by similar systems. The question of local stability is approached in a standard, but partially ad hoc, fashion.
Because our system is two-dimensional, we can employ some Bendixson type reasoning to eliminate some very special kinds of cycles. In addition, the results on long term behavior of the virus and of the infected cell densities separately allow reasonable speculations about trajectories in two-dimensional phase space.
To simplify notation where possible, definite integrals will be written without the dummy variable or differential when not potentially confusing. First derivatives with respect to time will be indicated by a "prime" superscript. Equations or inequalities can appear consecutively on a line to save space. They usually form a logical sequence and are separated by semicolons. When subsequently cited, they are indicated by section and numeral enclosed in parentheses, as opposed to propositions, definitions, etc., for which parentheses are not used.
We wish to model infection by a virus which attacks a particular circulating host cell line, using infected cells to manufacture more virus. To counter the attack, we endow the host with several immune mechanisms to clear both virus and the infected cells from the circulation. By choosing a unit of volume, we can assume that the density of host cells in the blood volume is 1. We normalize time so that it takes one unit of time for all virions to confront a host cell and also for immune processes to function. We assume that the body immediately replaces infected host cells which it has cleared or which have been destroyed by the virus, so that the total number of host cells, including both those infected and those not, remains constant. (Clearly, this assumption is not reflective of chronic HIV.) At time t=0, an inoculum of virus with density ν0 is introduced into the formerly virus-free blood.
We denote the rate at which virions penetrate attacked cells by p, where 0<p≤1. The virus does not discriminate with regard to the infection status of the host cells it encounters and, once infected, a cell's status is unchanged by further viral breaches. Since the density of all host cells remains constant at 1, virus is lost to successful attack at the rate p. If an infected cell is not cleared sooner, it is lysed exactly T units of time after infection, releasing b virus particles back into the circulation, for positive constants T and b. We also assume a constant rate, l≥0, of leakage of virus from infected host cells into the circulation.
Since uninfected host cells may have immune function, we assume that the clearance of virus has two components, one dependent on the density of uninfected host cells with relative rate η(t) and another, from circulating antibodies and other immune mechanisms, with relative rate κ(t). For the infected cells, we assume a simple clearance independent of host cells at a relative rate μ(t). These functions, η, κ, and μ defined on (−T,∞), are continuous, take values in [0,1), and eventually become constant at values h, k, and m respectively. Although it is not necessary for the analysis, we might think of them as increasing to their plateau values as in the activation of immunity.
We let V(t) be the density of free circulating virions at time t (equivalently, the number of virions per host cell) and W(t) the density of those host cells which have been infected (warped?). We define E by E(t)=pV(t)(1−W(t)), the rate at which virus infects previously uninfected host cells. Virions are cleared by intact host cells at the rate ηV(1−W)=ηpE, and by host cell independent mechanisms at the rate κV. Their number is increased by viral leakage at the rate lW along with the expulsion of virus at infected cell lysis. The following pair of first order delay-differential equations along with the initial conditions for V and W on [−T,0] describe the model.
V'(t)=−η(t)pE(t)−(p+κ(t))V(t)+lW(t)+be−∫tt−TμE(t−T) | (2.1a) |
W'(t)=E(t)−μ(t)W(t)−e−∫tt−TμE(t−T) | (2.1b) |
V(t)≡0on(−T,0),V(0)=ν0>0,andW(t)≡0on(−T,0] | (2.1c) |
Note that although V has a jump discontinuity at 0, both V and W are constant on (−T,0). The unique solution to the system with these initial conditions will become progressively smoother at the values t=nT for positive integers n.
Let Σ be the closed strip in two-dimensional phase space defined by V(t)≥0 and 0≤W(t)≤1. The following proposition establishes bounds on the region in which the trajectory of (V(t),W(t)) lies.
Proposition 2.1: Let γ=lT+beT. Then, for all t>0,
(i) (V(t),W(t))∈Σ°, the interior of Σ.
(ii) V(t)≤v0+γ+max{γ,γpT}.
(iii) W(t)≤1−e−pvST1+pvST, where vS=supV(t) (which exists by ii).
Proof: (i) The initial conditions (2.1c) state that V(0)>0 and W(0)=0. By (2.1b), we also have W'(0)=E(0)=pV(0)>0. Therefore, by continuity and these inequalities, (V(t),W(t))∈Σ° on an interval (0,u), for some u>0. Suppose u is maximal with this property implying, by continuity and the fact that Σ° is open, that (V(u),W(u)) lies on the boundary of Σ. Let νu=maxt∈[0,u]V(t). Then, by (2.1b), for t∈[0,u),
W'(t)≤E(t)≤pνu(1−W(t));−W'(t)1−W(t)≥−pνu |
Integrating from 0 to u, using continuity of W for positive t, and then exponentiating, we get
log(1−W(u))≥−pνuu;1−W(u)≥e−pνuu;W(u)≤1−e−pνuu<1 |
Next, for t∈(0,u), by (2.1a), the definition of E(t), and what we have just shown, and letting H and K be upper bounds for η and κ respectively,
V'(t)≥−(η(t)+p+κ(t))V(t)≥−(H+p+K)V(t);V'(t)V(t)≥−(H+p+K) |
Integrating from 0 to u, and then exponentiating as above, we get
V(u)≥v0e−(H+p+K)u>0 |
Finally, we take a derivative and use (2.1b) again, yielding for t>0,
(e∫t0μW(t))'=e∫t0μ(W'(t)+μ(t)W(t))=e∫t0μ(E(t)−e−∫tt−TμE(t−T)) |
=e∫t0μE(t)−e∫t−T0μE(t−T) |
If we let g(t) be the first term on the right-hand side of the last equality (so the second is g(t−T)), note that g(t)≡0 on (−T,0], and integrate from 0 to t, using simple substitution on the right, we get
e∫t0μW(t)=∫t0g(x)dx−∫t0g(x−T)dx=∫t0g−∫t−T−Tg=∫tt−Tg−∫0−Tg=∫tt−TgW(t)=e−∫t0μ∫tt−Te∫x0μE(x)dx=∫tt−Te−∫txμE(x)dx | (2.2) |
and, since E>0 on (0,u), W(u)>0. Therefore, (V(u),W(u))∈Σ°, a contradiction, so (V(t),W(t))∈Σ° for all t>0.
(ii) Next, we construct a special sequence of times, 0=t0<t1<⋯<tn<… which increases without bound. Assuming we have determined tn, if V(tn)≤V(t) for t∈[tn,tn+T], take tn+1=tn+T and, otherwise, let tn+1 be any value in [tn,tn+T] for which V(t) attains its minimum on that closed interval. In either case tn+1−tn≤T and tn+2≥tn+T so, indeed, tn→∞. We first wish to show by induction that for all n, V(tn)<v0+max{γ,γpT}. Since V(t0)=V(0)=v0, this is true for n=0, so suppose it true for some n≥0. If tn+1 was determined as a minimum for V(t) on [tn,tn+T], then V(tn+1)<V(tn), so the statement is true for n+1. Otherwise, tn+1=tn+T and tn is a minimum for V(t) on [tn,tn+T]. The following inequalities are consequences of (2.1a) and (2.2) and the facts that W(t)<1 and μ(t)<1 for all t.
V'(t)≤−pV(t)+l+bE(t−T) | (2.3) |
1>W(t)≥e−T∫tt−TE;eT>∫tt−TE=∫t+TtE(x−T)dx | (2.4) |
Integrating (2.3) from tn to tn+1 and using (2.4), we get
V(tn+T)−V(tn)≤−p∫tn+TtnV+∫tn+Ttnl+b∫tn+TtnE(x−T)dx |
≤−p∫tn+TtnV(tn)+lT+beT. |
V(tn+1)≤(1−pT)V(tn)+γ |
If pT≥1, V(tn+1)≤γ, and the statement is true for n+1. If pT<1, then, γ<γpT, and by the induction assumption,
V(tn+1)≤(1−pT)(v0+γpT)+γ=(1−pT)v0+(1−pT)γpT+γ<v0+γpT |
and, again, the statement is true. Now any positive t∈ [tn,tn+1] for some n≥0. Since tn+1−tn≤T, weakening (2.3) above to V'(t)≤l+bE(t−T), integrating from tn to t, and then using (2.4), we get
V(t)−V(tn)≤l(t−tn)+b∫ttnE(x−T)dx<lT+beT=γ |
which finishes the proof of ii.
(iii) First, we claim that for any positive t, there is an t0∈[t−T,t] such that W(t0)≤pvST1+pvST. Suppose this were not true. Then, we would have for all x∈[t−T,t],
W(x)>pvST1+pvST;1−W(x)<11+pvST |
so, by (2.2) and the fact that t∈[t−T,t], implying that W(t)>pvST1+pvST,
pvST1+pvST<W(t)<∫tt−TE≤pvS∫tt−T(1−W)<pvST1+pvST |
a contradiction, so there is, indeed, an t0∈[t−T,t] such that 1−W(t0)≥11+pvST. Now, integrating the inequality −W'(t)1−W(t)≥−pνS, which follows from (2.1b) as in the proof of part i, this time from t0 to t, and exponentiating, we have
1−W(t)1−W(t0)≥e−pνS(t−t0)≥e−pνST |
1−W(t)≥(1−W(t0))e−pνST≥e−pvST1+pvST;W(t)≤1−e−pvST1+pvST |
Proposition 2.1 places V and W and, because they satisfy Eqs (2.1), their derivatives as well, in the real algebra of bounded C1 functions, so they are all Lipschitz. (A function f is Lipschitz if there is an M>0 such that for all t,u in its domain |f(t)−f(u)|≤M|t−u|).
Let τ be so large that for t≥τ−T, μ(t)≡m, η(t)≡h, and κ(t)≡k, making the system autonomous with constant coefficients. We also assume that τ>2T to guarantee that V and W are of class C2 for t≥τ. It will be convenient to define a=p+k>0, and c=pbe−mT−h. Also, let F(t)=e−mT(E(t)−E(t−T)), q=p(1−e−mT)m for m>0, and q=pT if m=0. (Interestingly, limm→01−e−mTm=T, by L'Hôspital's Rule and, as a function of m, q is decreasing, so q≤pT). Let J(t)=qV(t)(1−W(t))−W(t). Then, for t≥τ, the system becomes
V'(t)=−hpE(t)−aV(t)+lW(t)+be−mTE(t−T)=cV(t)(1−W(t))−aV(t)+lW(t)−be−mT(E(t)−E(t−T))=cV(t)(1−W(t))−aV(t)+lW(t)−bF(t)=(c+lq)V(t)(1−W(t))−aV(t)+l(W(t)−qV(t)(1−W(t)))−bF(t) |
=V(t)(c+lq−a−(c+lq)W(t))−lJ(t)−bF(t) | (3.1a) |
W'(t)=E(t)−mW(t)−e−mTE(t−T) |
=(1−e−mT)E(t)−mW(t)+F(t) |
=m(qV(t)(1−W(t))−W(t))+F(t)=mJ(t)+F(t) | (3.1b) |
Although the last form of (3.1b) was derived with the case m>0 in mind, it remains valid for m=0, in which case W'(t)=E(t)−E(t−T)=F(t). The purpose of the term J=qV(1−W)−W in the final forms of defining Eqs (3.1) will become apparent shortly. We can also rewrite (2.2) as
W(t)=∫tt−Te−m(t−x)E(x)dx=p∫tt−Te−m(t−x)V(x)(1−W(x))dx | (3.2) |
By a fixed point, (v∗,w∗), of the system, we mean a constant solution V(t)≡v∗, W(t)≡w∗ to Eqs (3.1) on (−T,∞). We are interested in them because any point in Σ to which (V(t),W(t)) might converge must be a fixed point.
Proposition 3.1: Fixed points for the system include (0,0) and, if and only if c+lq>a, (v∗,w∗)=(c+lq−aqa,c+lq−ac+lq).
Proof: By the third form of (3.1a), the last form of (3.1b) for m>0, and (3.2) for m=0, (v,w)∈Σ is a fixed point of the autonomous system if and only if
v(c(1−w)−a)+lw=0 |
qv(1−w)−w=0ifm>0 |
w=pTv(1−w)=qv(1−w)ifm=0 |
But the second and third equations are identical, making m moot. Now, (0,0) is always a fixed point and, if either v or w is zero, the other is as well by the second equation. Suppose there is another (nonzero) fixed point (v∗,w∗)∈Σ. Then by the first, respectively second, equation,
c(1−w∗)+lw∗v∗=a;c−a=w∗v∗(cv∗−l) | (3.3a) |
w∗v∗=q(1−w∗) | (3.3b) |
Substituting (3.3b) into the first form of (3.3a), we get (c+lq)(1−w∗)=a so, since w∗>0, a second fixed point exists if and only if c+lq>a and, in that case,
1−w∗=ac+lq;w∗=1−ac+lq=c+lq−ac+lq | (3.3c) |
v∗=w∗q(1−w∗)=c+lq−aqa;qv∗+1=c+lq−aa+aa=c+lqa=11−w∗ | (3.3d) |
From now on, whenever c+lq>a, (v∗,w∗) will denote the nonzero fixed point given by this proposition.
Lemma 3.2: There are positive constants d1 and d2 such that for all t≥τ,
d1(V(t))2≤W(t)≤d2V(t) |
Proof: By the first form of (3.1a), and the fact that W≤1, V'(t)≥−(h+a)V(t). For x∈[t−T,t], integrating from x to t and exponentiating, we get
V(t)V(x)≥e−(h+a)(t−x)≥e−(h+a)T;V(x)≤e(h+a)TV(t) |
Therefore, by (3.2), we have
W(t)≤p∫tt−TV≤pTe(h+a)TV(t) |
so we can take d2=pTe(h+a)T. To prove the left-hand inequality, let d be the larger of vST and a Lipschitz constant for V, and observe that for all t≥τ, the line passing through (t,V(t)) with slope d must cross the t-axis at some u≥t−V(t)vST and also lie below the plot of the function V on [t−T,t]. Therefore, ∫tt−TV exceeds the area of the triangle with height V(t) and base V(t)d, i.e., (V(t))22d. Again, employing (3.2) along with Proposition 2.1iii, we have
W(t)≥pe−mT∫tt−TV(1−W)≥pe−(m+pvS)T1+pvST∫tt−TV≥pe−(m+pvS)T2d(1+pvST)(V(t))2 |
Definition 3.3: A method we shall use several times in the proofs of this section is the exploitation of a different kind of bound on integrable functions defined on [τ,∞). We shall say such a function, f, is integrally bounded if there is a Cf>0 such that |∫tuf|≤Cf for all τ≤u<t and, in that case, we write f∼0, a definition relevant only for this paper. For example, f(t)=sint∼0 and its terminal behavior is cyclic, while f(t)=1t even converges to zero, but is not integrally bounded. It is straightforward to demonstrate that if f∼0 and g∼0, then αf+βg∼0 for any α,β∈R. Even further, if f∼0 and g is any bounded function, then gf∼0, and if g is bounded away from zero, then fg∼0. The derivative of any bounded function is integrally bounded, as we have ∫tuf'=f(t)−f(u). We shall also say that f∼g if f−g∼0, so that if f∼g and g∼h, then f∼h.
Lemma 3.4: If f:[τ,∞)→R is bounded, then f(t)∼f(t−T). In particular, F(t)=e−mT(E(t)−E(t−T))∼0.
Proof: For any τ≤u<t, we have
|∫tu(f(x)−f(x−T))dx|=|∫tuf−∫t−Tu−tf|=|∫tt−Tf−∫uu−Tf|≤2Tsup|f| |
Lemma 3.5: (Barbãlat's Lemma): Let t0≥τ, g:[t0,∞)→R be differentiable with uniformly continuous derivative. Then, if limt→∞g(t) exists (and is finite), limt→∞g'(t)=0.
Proof (included for completeness): Suppose g'(t) does not converge to 0 as t→∞, so that there is an ϵ>0 such that |g'(t)|≥ϵ for arbitrarily large values of t. Since g' is uniformly continuous, there is a δ>0 such that |g'|≥ϵ2 on the interval [t,t+δ] for all such t and, consequently, |g(t)−g(t+δ)|≥ϵδ2, which precludes the convergence of g(t).
The next corollary will be applied repeatedly to functions which are both integrally bounded and eventually do not change sign. It is the main reason for defining the equivalence ∼.
Corollary 3.6: Let t0≥τ and g:[τ,∞)→R be Lipschitz and g∼0. Then either limt→∞g(t)=0 or g changes sign for arbitrarily large values of t.
Proof: Suppose that t0≥τ is such that g takes only nonnegative or nonpositive values on [t0,∞) and let G:[t0,∞)→R be defined by G(t)=∫tt0g. Then G is monotone and, since g∼0, it is bounded, so G converges. But, G'(t)=g(t) and, since g is Lipschitz, it is uniformly continuous. Lemma 3.5 applies.
Lemma 3.7: Let f:[τ,∞)→R be continuous and bounded. Then, ∫xx−Tf∼Tf(x).
Proof: This lemma states that the difference between f and its rolling uniform mean over the prior T time units is integrally bounded. Let t−2T>u≥τ and consider ∫tu∫xx−Tf(y)dydx. It will be helpful to visualize the region in the (x,y)-plane which corresponds to this iterated integral. It is a parallelogram whose bases are the horizontal line segments extending to the right from (u−T,u) to (u,u) and (t−T,t) to (t,t) and whose oblique sides are segments of the lines y=x and y=x−T. If we lop off triangles on the top right and bottom left which have height and base T, we are left with another parallelogram bounded by the same oblique lines but with vertical bases extending up from (u,u) to (u,u+T) and (t−T,t−T) to and (t,t−T). The integral of f(y) over those two small removed triangles is clearly bounded regardless of u and t and the integral of f(y) over the smaller parallelogram as an iterated integral is
∫t−Tu+T∫y+Tyf(y)dxdy=∫t−Tu+TTf=∫tuTf−∫uu+TTf−∫tt−TTf |
but the last two terms on the right are also bounded irrespective of u and t.
Remark 3.8: The last form of (3.1b) along with Definition 3.3 and Lemma 3.4 immediately yield for m>0
(i) J∼0
But, if m=0, by (3.2) and Lemma 3.7, W=p∫tt−TV(1−W)∼pTV(1−W)=qV(1−W), so i remains valid. Applying it and Lemma 3.4 to the last form of (3.1a), we get
(ii) V(c+lq−a−(c+lq)W)∼0
Definition 3.9: When c+lq>a, define, for all t≥τ, A(t)=V(t)(w∗−W(t)), B(t)=W(t)(v∗−V(t)), and D(t)=w∗v∗V(t)−W(t). Notice that A−B=v∗D.
Lemma 3.10: If c+lq>a, then
(i) A∼0 and B∼0
(ii) V'(t)=1v∗((cv∗−l)A+lB)−bF(t)
Proof: (i) Remark 3.8ii and the second form of (3.3c) tell us that A∼0. By Remark 3.8i and (3.3b), we have
0∼J=qV(w∗−W+1−w∗)−W=qA+q(1−w∗)V−W=qA+D |
so D∼0. But B=A−v∗D, so B∼0.
(ii) We have, by the second form of (3.1a), and the first form of (3.3a),
V'(t)=cV(t)(w∗−W(t)+1−w∗)−aV(t)+lW(t)−bF(t)=cA(t)+(c(1−w∗)−a)V(t)+lW(t)−bF(t)=cA(t)−lw∗v∗V(t)+lW(t)−bF(t)=cA(t)+l(W(t)−w∗v∗V(t))−bF(t)=cA(t)−lD(t)−bF(t) |
But, we also have
(cv∗−l)A+lB=cv∗A−l(A−B)=cv∗A−lv∗D |
completing the proof.
We are now in a position to exhibit a criterion for the persistence versus eventual extinction of the virus.
Theorem 3.11: (i) If c+lq≤a, then (V(t),W(t))→(0,0) as t→∞.
(ii) If c+lq>a and c≥a (either c>a or c=a and l>0), then V(t) and W(t) are bounded above zero.
Proof: (i) Suppose that 0<c+lq≤a. By the right-hand side of Lemma 3.2, it is adequate to show that V(t)→0. By Remark 3.8ii, V(c+lq−ac+lq−W)∼0. If c+lq<a, then c+lq−ac+lq−W is bounded below zero which, by Definition 3.3, implies that V∼0 and, by Corollary 3.6, V(t)→0. If c+lq=a, then VW∼0 and, this time, V(t)W(t)→0. But by the left-hand side of Lemma 3.2, V(t)W(t)≥d1(V(t))3, so V(t)→0.
(ii) As with part i, it is adequate to show that W is bounded above zero, so suppose it is not, and that there is an unbounded strictly increasing sequence of times, {tn}, n≥1, with W(tn)→0. Therefore, we can find a t0≥τ such that W(t0)≤min{w∗,d1(v∗)2} and, by replacing {tn} with a terminal subsequence if necessary, assume that tn>t0 and W(tn)<e−mTW(t0) for all n≥1. Now define a second sequence {un}, n≥1, by letting un be the largest number less than tn such that W(un)=W(t0). (It certainly may be the case that {un} is not strictly increasing and it could even have a maximal element.) Using (3.2), we derive the string of inequalities
∫tntn−TE≤emT∫tntn−Te−m(tn−u)E(u)du=emTW(tn)≤W(t0)=W(un) |
=∫unun−Te−m(un−u)E(u)du≤∫unun−TE |
by which, looking back at the proof of Lemma 3.4, we can see that ∫tnunF≤0. Now, by the definition of un, on the interval [un,tn], W≤w∗ and, by Lemma 3.3, d1V2≤W≤d1(v∗)2, so V≤v∗. Recall that Lemma 3.10ii states
V'(t)=(c−lv∗)V(t)(w∗−W(t))+lv∗W(t)(v∗−V(t)))−bF(t) |
Integrating this equation from un to tn gives us
V(tn)−V(un)=(c−lv∗)∫tnunV(w∗−W)+lv∗∫tnunW(v∗−V)−b∫tnunF |
Note that by our assumption that c≥a and the second form of (3.3a), cv∗−l≥0. By our definition of {un} then, the right-hand side is nonnegative, so V(un)≤V(tn). But by Lemma 3.2, V(tn)→0, so, V(un)→0 and, by the same lemma, W(t0)d2=W(un)d2≤V(un), a contradiction.
This theorem nails down the gross behavior of V and W except in the case that c<a<c+lq, which doesn't happen if l=0, i.e., there is no pre-lysis leakage. As per the definitions at the beginning of this section, c−a=p(be−mT−1)−(h+k), which looks like net production minus destruction of virus, but discounting the contribution of leakage. Curiously, our proof of the theorem relies on its nonnegativity, not just the positivity of c+lq−a.
Corollary 3.12: Suppose c+lq>a and c≥a. Then
(i) w∗−W∼0 and v∗−V∼0.
(ii) If v∗−V, respectively w∗−W, does not change sign at arbitrarily large values of t, then V(t)→v∗, respectively W(t)→w∗.
Proof: (i) We already know, by Lemma 3.10i, that V(w∗−W)∼0 and W(v∗−V)∼0, but by the theorem, V and W are bounded above zero, and the result follows from Definition 3.3.
(ii) This follows from part i and Corollary 3.6.
Part ii of the corollary says that if a trajectory does not converge, then V(t) passes back and forth through v∗ or W(t) through w∗ forever. This means that a trajectory can only converge, converge in one coordinate and oscillate in the other, or oscillate in both, the last of which by no means proves, but does suggest, that the trajectory circles (v∗,w∗). Corollary 3.15 will narrow these options a bit.
Next, we will investigate further the long-term convergence-type behavior of V and W. We shall be looking at their behavior on (long) intervals of the form I=[u0,t0] where τ≤u0. We'll be manipulating rational functions of V(t),V(t−T),W(t),W(t−T) with denominators bounded away from zero, all of which are bounded Lipschitz functions with bounded Lipschitz derivatives.
Lemma 3.13: Let ϵ>0 be given and suppose that for some C=C(ϵ)>0, f:I=[u0,t0]→R satisfies |f(t0)|<ϵ and f'(t)=−C(f(t)−ε(t)) on I, where |ε(t)|<ϵ. Then |f(t)|<ϵ on I.
Proof: Suppose, to the contrary, that f(t)=ϵ for some t∈I, where t is as small as possible. Then t≠u0 and, since f<ϵ on the interval [u0,t), it must be the case that f'(t)≥0. However, this is contradicted by
f'(t)=−C(f(t)−ε(t))=−C(ϵ−ε(t))<0 |
Similarly, f(t)>−ϵ for t∈I.
Theorem 3.14: If c+lq>a, c≥a, and there are arbitrarily long intervals on which w∗−W doesn't change sign, then there are arbitrarily long intervals on which (V(t),W(t)) is arbitrarily close to (v∗,w∗).
Proof: First, we claim that there are arbitrarily long intervals on which W(t) is arbitrarily close to w∗ and W'(t) is arbitrarily close to zero as well. Assume that there are arbitrarily long intervals on which W≥w∗. The proof is similar for W≤w∗ and one or the other must be true for infinitely many such intervals. Let ϵ>0 be given, let M be a Lipschitz constant for W, and let Λ be the length of an interval on which we wish to establish that W−w∗<ϵ. Generality is not lost if we assume that Λ≥ϵM, so consider such an interval, I. Reasoning as in the proof of Lemma 3.2, for any t∈I, if W(t)−w∗=ϵ, then the right triangle with height ϵ and base either [t,t+ϵM] or [t−ϵM,t], whose area is ϵ22M lies beneath the plot of W(t)−w∗. Now, by Corollary 3.12i, there is a C>0 such that |∫tu(W−w∗)|≤C for all τ≤u<t. Let the integer n be so large that nϵ22M>C and let τ≤u0<t0 such that W(t)−w∗≥0 on [u0,t0] and t0−u0=nΛ. Suppose that W(t)−w∗ took the value ϵ on all of the subintervals [u0+jΛ,u0+(j+1)Λ], j=0,…,n−1. Then we would have the contradiction
nϵ22M≤∑n−1j=0∫u0+(j+1)Λu0+jΛ(W−w∗)=∫t0u0(W−w∗)≤C |
so, it must be the case that W(t)−w∗<ϵ on at least one of these subintervals of length Λ.
Now, let ϵ>0 be given and M be a Lipschitz constant for W'. By what has been done above, we can find an interval I of length at least ϵM and arbitrarily large on which |w∗−W|<min{ϵ,ϵ28M}, immediately implying that |W(t)−W(u)|<ϵ24M for any t,u∈I. Suppose first that W'(t)≥ϵ for some t∈I, Then at least one of t±ϵ2M∈I, say t+ϵ2M∈I. By the definition of M, W'≥ϵ2 on [t,t+ϵ2M], so W(t+ϵ2M)−W(t)≥ϵ24M, a contradiction. Of course, the reasoning is analogous for W'≤−ϵ and/or t−ϵ2M∈I. The proof of this claim is reminiscent of that of Lemma 3.5 (Barbãlat's Lemma). Now, we attend to the rest of the theorem.
First, we consider the case m>0 and let ϵ>0 and Λ>0 be given. Let the integer n be so large that e−nmTv∗<ϵ. The first form of (3.1b) can be written
V(t)−e−mTV(t−T)1−W(t−T)1−W(t)=W'(t)+mW(t)p(1−W(t)) |
Now, by the first form of (3.3d) and the definition of q, mw∗p(1−w∗)=(1−e−mT)v∗ and
1−W(t−T)1−W(t)−1=1−W(t−T)−1+W(t)1−W(t)=W(t)−W(t−T)1−W(t)=(W(t)−w∗)+(w∗−W(t−T))1−W(t) |
W'(t)+mW(t)p(1−W(t))−mw∗p(1−w∗)=(1−w∗)(W'(t)+mW(t))−(1−W(t))mw∗p(1−W(t))(1−w∗)=(1−w∗)W'(t)+m(W(t)−w∗)p(1−W(t))(1−w∗) |
Therefore, looking at the last terms of the two equations above, we see that there is a C>0 which does not depend on ϵ such that if I is an interval on which |w∗−W(t)|<ϵ, |W'(t)|<ϵ, and t−T,t∈I, then
V(t)−e−mTV(t−T)=(1−e−mT)v∗+ε(t) |
where |ε(t)|<Cε. By what we have done above, there is, in fact, such an interval, I, of length Λ+nT. Let t be in the terminal subinterval of I with length Λ. Then, we have a telescoping sum and simplifications
V(t)−e−nmTV(t−nT)=∑n−1j=0e−jmT(V(t−jT)−e−mTV(t−(j+1)T)) |
=∑n−1j=0e−jmT((1−e−mT)v∗+ε(t−jT)) |
=(1−e−nmT)v∗+∑n−1j=0e−jmTε(t−jT) |
|V(t)−v∗|≤e−nmTv∗+Cϵ∑n−1j=0e−jmT<ϵ+Cϵ1−e−mT=(1+C1−e−mT)ϵ |
Next suppose that m=0, ϵ>0 and Λ>0 have been given, and let I be an interval of length Λ on which |w∗−W(t)|<ϵ and |W'(t)|<ϵ. The first form of (3.1b) is W'(t)=E(t)−E(t−T)=F(t), so |F(t)|<ϵ, as well. By (3.2), we have
W(t)=p∫tt−TV(1−W)=p∫tt−TV(1−w∗+w∗−W)=p(1−w∗)∫tt−TV+p∫tt−TA |
W(t)p(1−w∗)−∫tt−TV=∫tt−TA(t)1−w∗;w∗p(1−w∗)−∫tt−TV=∫tt−TA(t)1−w∗−W(t)−w∗p(1−w∗) |
By the first form of (3.3d), v∗T=w∗p(1−w∗) so, using the last equation,
|v∗T−∫tt−TV|=|∫tt−TA(t)1−w∗−W(t)−w∗p(1−w∗)|<(vST1−w∗+1p(1−w∗))ϵ |
Now, suppose further that l=0. By Lemma 3.10ii, V'(t)=cA(t)−bF(t), so that |V'(t)|<(cvS+b)ϵ. Therefore, for any t−T,t∈I, if x∈[t−T,t], then |V(t)−V(x)|<(cvS+b)Tϵ. Integrating over [t−T,t] and then using the last inequality just above, we have
|TV(t)−∫tt−TV|≤∫tt−T|V(t)−V(x)|dx<(cvS+b)T2ϵ |
Then, by this inequality and the one just above,
|V(t)−v∗|<((cvS+b)T+vS1−w∗+1q(1−w∗))ϵ |
Finally, suppose that m=0 but l>0, and I is as above. By Lemma 3.10ii,
(V(t)−v∗)'=V'(t)=(c−lv∗)A(t)+lv∗B(t)−bF(t) |
B(t)=W(t)(v∗−V(t))=w∗(v∗−V(t))+(W(t)−w∗)(v∗−V(t)) |
=w∗(v∗−V(t))+ζ(t) |
where |ζ(t)|<vSϵ, so
(V(t)−v∗)'=lw∗v∗(v∗−V(t))+lv∗ζ(t)+(c−lv∗)A(t)−bF(t) |
(V(t)−v∗)'=−lw∗v∗(V(t)−v∗)+ε(t) |
where |ε(t)|<Cϵ for some C>0. Now, we could have picked I such that |w∗−W(t)|<lw∗v∗Cϵ and |W'(t)|<lw∗v∗Cϵ on I, and would have deduced the last equation above except with |ε(t)|<lw∗v∗ϵ. We now could apply Lemma 3.13 if we knew that |V(u0)−v∗| was sufficiently small, where u0 was the left endpoint of I. But by Corollary 3.12i, for any ϵ>0, V(t) cannot remain at least ϵ distant from v∗ on I if I is long enough. Therefore, taking I twice that long allows us to assume that |V(t0)−v∗| is small enough.
Corollary 3.15: If c+lq>a, c≥a, and W converges, then (V(t),W(t))→(v∗,w∗).
Proof: By Corollary 3.12ii, if W converges, then W(t)→w∗. Then, by Lemma 3.5 W'(t)→0. This allows us to find for any ϵ>0 an interval I=[u0,∞) on which both |w∗−W(t)|<ϵ and |W'(t)|<ϵ. Using the reasoning in the proof of the theorem, we conclude that |V(t)−v∗|<ϵ on [t0,∞) for some t0≥u0.
Corollary 3.16: If c+lq>a, c≥a, and (v∗,w∗) is locally asymptotically stable, then either (V(t),W(t))→(v∗,w∗) or the distance between successive values of t for which W(t)=w∗ is bounded.
Because of this corollary, although it wouldn't guarantee convergence, it would be nice to know if (v∗,w∗) were locally asymptotically stable when c+lq>a. In Section 5, we shall provide some conditions under which that is the case. We have not been able to prove results analogous to Theorem 3.14, respectively Corollary 3.15, whose hypothesis is that v∗−V doesn't change sign over arbitrarily long intervals, respectively converges.
One of the questions in the study of analogous systems of differential equations without delay is the existence of cyclic solutions. A priori, we know of no reason why there could not be cyclic solutions, of any period, to (3.1a) and (3.1b). The existence of such solutions suggests other solution trajectories which spiral into a cycle, or spin away from one. It would be interesting to know if ours was one of those.
For (V,W) to be a cyclic solution of period T it would be necessary, but far from sufficient, that V(−T)=V(0) and W(−T)=W(0). Our solution couldn't be one because V(−T)=0≠v0=V(0) (even though W(−T)=W(0)=0), but it could spiral toward or away from one.
It is possible, however, to specify a delay-differential system, derived from a differential system having cyclic solutions of a period T, which has those same solutions. Simply take initial conditions to be such a cycle with any starting point. Then add functions of the form f(t)−f(t−T) to the right-hand side of the differential system, where f depends only on the solutions.
The extreme case of T=0 provides an "undelayed system" of differential equations,
V'(t)=cV(t)(1−W(t))−aV(t)+lW(t) |
W'(t)=(1−e−mT)E(t)−mW(t) |
obtained from Eqs (3.1) by removing F(t) from the right-hand side (not removing only E(t−T)).
One commonly employed method for ruling out the possibility of cycles in two-dimensional phase space is the Bendixson-Dulac criterion. The simplest form of it would eliminate cycles of any period if, where V'and W'are as in the equations above, ∂V'∂V+∂W'∂W does not take both signs on Σ°. In fact,
∂V'∂V+∂W'∂W=∂∂V(cV(1−W)−aV+lW)+∂∂W(p(1−e−mT)V(1−W)−mW) |
=c(1−W)−a−p(1−e−mT)V−m |
Therefore, there are no cycles if c≤a. If c+lq≤a as well, we know that all solutions converge to zero, but we still don't know much about outcome if c<a<c+lq, except for this result, which rules out only a very specific type of cycle.
In this section, we shall look at the stability of the fixed points (0,0) and, when c+lq>a and c≥a, (v∗,w∗). This translates to determining the characteristic function, G, of the system linearized around its fixed points and finding conditions which guarantee that either some or none of its zeros have positive real parts. If all zeros of G have negative real part, then the fixed point will be locally asymptotically stable. We use the second form of (3.1a) for V' and the last form of (3.1b) for W'.
V'(t)=cV(t)(1−W(t))−aV(t)+lW(t)−be−mT(E(t)−E(t−T)) |
W'(t)=m(qV(t)(1−W(t))−W(t))+e−mT(E(t)−E(t−T)) |
Observe that
∂∂V(cV(1−W)−aV+lW)=c(1−W)−a;∂∂W(cV(1−W)−aV+lW)=l−cV |
∂∂V(qV(1−W)−W)=q(1−W);∂∂W(qV(1−W)−W)=−(1+qV) |
∂∂VE=∂∂VpV(1−W)=p(1−W);∂∂WE=−pV |
Linearizing around (0,0) but leaving out the last terms and subtracting λI, we get the matrix
A=(c−a−λlmq−m−λ) |
|A|=(c−a−λ)(−m−λ)−mlq=λ2+(m−(c−a))λ−m(c+lq−a) |
This determinant is the characteristic function for the undelayed system (T=0) and clearly, if a<c+lq, it has a positive root, making (0,0) unstable for any delay-differential system whose undelayed system is this one. Otherwise, by Theorem 3.11i, all trajectories converge to (0,0).
Next, assuming there is a positive fixed point and linearizing around (v∗,w∗) but leaving out the last terms and subtracting λI, we get the matrix
A=(c(1−w∗)−a−λl−cv∗mq(1−w∗)−(1+qv∗)−λ)=(−lw∗v∗−λl−cv∗mw∗v∗−m1−w∗−λ) |
where the second equality follows from the first form of (3.3a), (3.3b), and the second form of (3.3d). Again, A is the matrix of the linearization of the undelayed system. Therefore,
|A|=(lw∗v∗+λ)(m1−w∗+λ)−mw∗v∗(l−cv∗) |
=λ2+(m1−w∗+lw∗v∗)λ+mlw∗v∗(1−w∗)−mw∗v∗(l−cv∗) |
Looking on the right and using equations (3.3d) and second form of (3.3c),
mlw∗v∗(1−w∗)−mw∗v∗(l−cv∗)=mlw∗v∗(11−w∗−1)+cmw∗=mlqw∗+cmw∗ |
=(c+lq)mw∗=m(c+lq−a)=maqv∗ |
|A|=λ2+(m1−w∗+lw∗v∗)λ+maqv∗ |
We need to add to A the matrix corresponding to the last terms of the defining equations above. We have observed above that ∂∂VE=p(1−W) and ∂∂WE=−pV, so this matrix is pe−mT(1−e−λT)B, where
B=(−b(1−w∗)bv∗1−w∗−v∗) |
Then, G(λ)=|A+pe−mT(1−e−λT)B|. To calculate G(λ), we use the multilinearity on rows of the determinant, noting that |B|=0
G(λ)=|A|+pe−mT(1−e−λT)(−|lw∗v∗+λcv∗−l1−w∗−v∗|−b|1−w∗−v∗mw∗v∗−m1−w∗−λ|) |
−|lw∗v∗+λcv∗−l1−w∗−v∗|−b|1−w∗−v∗mw∗v∗−m1−w∗−λ|=|1−w∗−v∗lw∗v∗+λcv∗−l|+|1−w∗−v∗−mbw∗v∗mb1−w∗+bλ| |
=|1−w∗−v∗(l−mb)w∗v∗+λcv∗−l+mb1−w∗+bλ|=λ|1−w∗−v∗1b|+ |
|1−w∗−v∗(l−mb)w∗v∗cv∗−l+mb1−w∗| |
=(b(1−w∗)+v∗)λ+(1−w∗)(cv∗−l)+mb+(l−mb)w∗ |
=(b(1−w∗)+v∗)λ+mb(1−w∗)+(cv∗−l)(1−w∗)+lw∗ |
Note that the constant term is strictly positive as, either l>0or (second form of (3.3a)) cv∗−l=w∗v∗(c−a)>0. Using the second form of (3.3a) again, (3.3b), and the first form of (3.3d), we manipulate the far right-hand side to get
(cv∗−l)(1−w∗)=c−aq(1−w∗)(1−w∗)=c+lq−aq−l=a(c+lq−a)qa−l=av∗−l |
G(λ)=λ2+(m1−w∗+lw∗v∗)λ+maqv∗ |
+pe−mT(1−e−λT)((b(1−w∗)+v∗)λ+mb(1−w∗)+av∗−l(1−w∗)) |
To make our characteristic function look less unwieldy, let
α=pe−mT(b(1−w∗)+v∗)>0 |
β=mb(1−w∗)+av∗−l(1−w∗)b(1−w∗)+v∗>0 |
δ=m1−w∗+lw∗v∗ |
and note that δ=0 if and only if m=l=0. We now rewrite the characteristic function as
G(λ)=λ2+δλ+maqv∗+α(1−e−λT)(λ+β) | (5.1) |
=λ2+(δ+α)λ+(maqv∗+αβ)−αe−λT(λ+β) |
Initially, we deal with G in an ad hoc fashion. First, consider G restricted to [0,∞). The first form of (5.1) demonstrates that G is positive on (0,∞), and G(0)=maqv∗, so G(0)=0 if and only if m=0. In that case(v∗,w∗) is unstable, so we shall assume from now on that m>0. It's clear that G(−λ)=−G(λ), so we may restrict our search for other roots to those with positive imaginary parts. Suppose that λ is a zero of G, where λ=x+iy, x≥0 and y>0. The following equations represent the real part of G(λ)=0 and the imaginary part divided by y.
x2−y2+δx+maqv∗+α(x+β)(1−e−xTcosyT−yx+βe−xTsinyT)=0 | (5.2a) |
2x+δ+α(1−e−xTcosyT+x+βye−xTsinyT)=0 | (5.2b) |
Only parts ii and iv of the following proposition are directly useful for ruling out zeros of G which have positive real parts.
Proposition 5.1: Suppose that G(λ)=0, where λ=x+iy, x≥0 and y>0. Let f(w∗)=3π(3π−1)(1−w∗)+3π(3π−1)2(1−w∗)2+w∗(1−w∗)3. Then
(i) x(2yT−1)<β
(ii) δ<β<a;m<c(1−w∗)2<α(1−w∗)
(iii) 3π2<yT;x<β3π−1<a3π−1
(iv) T>3π2α(1+f(w∗))
Proof: (i) Since m>0, δ>0 so, for (5.2b) to be true, the last term on its left-hand side must be negative. Therefore,
1−e−xT(cosyT−x+βysinyT)<0;exT<cosyT−x+βysinyT |
It is a straightforward exercise in elementary calculus to show that for C∈R, the function cosz+Csinz is maximized when z=tan−1C and that the maximum is √1+C2. Therefore, it must be the case that
1+xT≤exT<√1+(x+βy)2 |
xT<√1+(x+βy)2−1=x+βy√1+(x+βy)2+1<x+β2y |
2xyT<x+β;x(2yT−1)<β |
(ii) Again, because the last term on the left-hand side of (5.2b) is negative, sinyT<0. Therefore, |cosyT|<1 and yT>π. To save space, let ζ=α(1−e−λT) and ξ=δ+ζ. Observe, since 1−e−λT=1−e−xTcosyT+ie−xTsinyT, that ζ and, therefore, ξ have positive real and negative imaginary parts. Rewrite G as a quadratic function of λ and apply the quadratic formula to get
G(λ)=λ2+ξλ+(maqv∗+ζβ);2λ=−ξ+√ξ2−4(maqv∗+ζβ)=−ξ+σ |
where σ is the correct one of the two possible choices for the square root. Since λ has nonnegative real part and positive imaginary part, 0<Reξ≤Reσ and Imξ<Imσ. Since the real part of ζ is positive, Reσ2=Re(ξ2−4(maqv∗+ζβ))<Reξ2. Therefore,
(Imξ)2=(Reξ)2−Reξ2<(Reσ)2−Reσ2=(Imσ)2;|Imξ|<|Imσ| |
Now, should Imσ<0, along with Imξ<0 and Imσ>Imξ, then |Imσ|<|Imξ|, contradicting the inequality immediately above. Therefore, Imσ>0 and so Imσ>−Imξ, yielding
−Imξ2=2(Reξ)(−Imξ)<2(Reσ)(Imσ)=Imσ2=Imξ2−4βImζ |
2βImζ<Imξ2=2(Reξ)(Imξ)=2(Reξ)(Imζ) |
β>Reξ=δ+Reζ=m1−w∗+lw∗v∗+α(1−e−xTcosyT)>m |
However, from its definition above, β≤mb(1−w∗)+av∗b(1−w∗)+v∗, a weighted average of a and m with strictly positive weights. Therefore, a>β>δ and, by the first form of (3.3a),
c(1−w∗)−m1−w∗=a−lw∗v∗−m1−w∗=a−δ>0 |
and by the definitions of c and α, α>c(1−w∗).
(iii) Now, suppose that 3π2≥yT, so that cosyT≤0. Since Reζ=α(1−e−xTcosyT), it follows that α≤Reζ. Therefore, by the first form of (3.3a),
β>δ+Reζ≥δ+α>m1−w∗+lw∗v∗+c(1−w∗)=m1−w∗+a |
which contradicts part ii. Consequently, cosyT>0, 3π2<yT and, by part i, x<β3π−1.
(iv) Returning to (5.2a) once more, with the added information that cosyT>0, we find
0<x2−y2+δx+maqv∗+α(x+y+β) |
y2−αy−(x2+(δ+α)x+(maqv∗+αβ))<0 |
y2−αy−α2g(x)<0 where g(x)=(xα)2+(1+δα)xα+βα+mα∙aα∙qv∗
y<α(1+√1+4g(x))2<α(1+1+2g(x))2=α(1+g(x)) |
By the second inequality of part ii and the fact that c≥a,
aα<ac(1−w∗)≤11−w∗ |
By the second inequality of part iii, the first inequality of part ii, and (3.3b),
xα<13π−1⋅aα;δα<aα;βα<aα;mα<aα;qv∗=w∗1−w∗ |
g(x)<3π(3π−1)(1−w∗)+3π(3π−1)2(1−w∗)2+w∗(1−w∗)3=f(w∗)≈1.1181−w∗+.132(1−w∗)2+w∗(1−w∗)3 |
Finally, by the first inequality of part iii and what we have just done,
T=yTy>3π2α(1+g(x))>3π2α(1+f(w∗)) |
In [1], there are results which pertain to characteristic functions of the form P(λ)−Q(λ)e−λT where P and Q are real polynomials, P has degree 2, its coefficient of λ is positive, and its constant term is nonnegative, as is the case for our characteristic function but, unfortunately, only when Q is a monomial of degree 2 or less. In [21], a more sophisticated method to rule out roots of transcendental polynomials having positive real parts is explored, including a proof of a generalization of the following result from polynomials to many analytic functions.
Theorem: Suppose P and Q are polynomials with real coefficients, without common pure imaginary root, and that (P−Q)(0)≠0. If f(y)=|P(iy)|2−|Q(iy)|2 has no positive roots and the roots of the real polynomial P−Q all have negative real parts, then for any T>0, all the zeros of P(z)−Q(z)e−Tz have negative real parts as well. Conversely, if P−Q has a root with nonnegative real part, then P(z)−Q(z)e−Tz has zeros with nonnegative real part for any T>0.
We apply this theorem to our situation. Let P(z)=z2+δz+maqv∗+α(z+β) and Q(z)=α(z+β), so that G(λ)=P(λ)−Q(λ)e−λT. Then the root of Q is real and (P−Q)(0)=maqv∗≠0. In addition, since δ>0 is the coefficient of z in P−Q, both roots of that polynomial have negative real parts. Therefore, we are left to deal with the f of the theorem.
P(iy)=−y2+maqv∗+αβ+i(δ+α)y;Q(iy)=α(β+iy) |
|P(iy)|2=y4+(maqv∗)2+α2β2−2y2(maqv∗+αβ)+2αβmaqv∗+y2(δ2+α2+2αδ) |
|Q(iy)|2=α2(β2+y2) |
|P(iy)|2−|Q(iy)|2=y4+(δ(δ+2α)−2(maqv∗+αβ))y2+(maqv∗)2+2αβmaqv∗ |
The positive real roots of this quartic are the square roots of the positive roots of
F(y)=y2+(δ(δ+2α)−2(maqv∗+αβ))y+(maqv∗)2+2αβmaqv∗ |
We need to find conditions under which F has no positive roots, i.e., either that
δ(δ+2α)≥2(maqv∗+αβ) |
orthatδ(δ+2α)<2(maqv∗+αβ)and |
(δ(δ+2α)−2(maqv∗+αβ))2−4maqv∗(maqv∗+2αβ)<0 |
δ2(δ+2α)2−4δ(δ+2α)(maqv∗+αβ)+4(maqv∗+αβ)2−4maqv∗(maqv∗+2αβ)<0 |
δ2(δ+2α)2−4δ(δ+2α)(maqv∗+αβ)+4α2β2<0 |
δ(δ+2α)−2(maqv∗+αβ)−2(maqv∗+αβ)+4α2β2δ(δ+2α)<0 |
Therefore, by the last inequality above, it is sufficient that
maqv∗+αβ≥2α2β2δ(δ+2α);1+maqv∗αβ≥2αβδ(δ+2α) |
δ(δ+2α)≥2αβ1+maqv∗αβ;δ2+2αδ−2αβ1+maqv∗αβ≥0 |
This will be true if and only if δ exceeds the positive root of the quadratic on the right-hand side. For the moment, let A=β1+maqv∗αβ. Then, the positive root of the quadratic is α(√1+2Aα−1). Rationalizing the numerator, we have
√1+2Aα−1=2Aα1+√1+2Aα<Aα |
Consequently, it suffices that δ≥A. This represents a modest improvement over Proposition 5.1ii. It is interesting that stability depends, in part, on the relative magnitudes of m and a, the strengths of clearance of infected cells versus clearance of virus, both of which are defensive features of the host.
The model of viremia proposed here is a relatively simple one, with a two-dimensional system describing viral and infected cell quantities over time, and a single discrete delay representing the time from cell infection until lysis. The parameters of the system include the size of virus inoculum ν0, rate of penetration of virus into cells p, rate of viral leakage from infected cells prior to lysis l, magnitude of the time lag until lysis T, number of virions expelled at lysis b, clearance rate of infected cells m, and two clearance rates of virions h and k.
The long-term behavior of the solution which exhibits dependence on parameter values other than ν0, is investigated with some success. The extinction of virus and infected cells, represented by (0,0) in phase space is always a fixed point of the system and, when it is the only one, extinction is, indeed, ensured. The only other possible situation, occurring when viral production "outweighs" destruction, specifically when p(be−mT−1)+lq>h+k, is the presence of a second, positive, fixed point (v∗,w∗) whose coordinates are fairly simple functions of the parameter values. In this case, our knowledge of the outcome is less precise. As long as a similar parameter inequality, p(be−mT−1)≥h+k, which is independent of leakage, is also satisfied, we can demonstrate persistence of virus and infected cells above a positive threshold, but we have been able to prove convergence to (v∗,w∗) only in special circumstances. However, we can show that (0,0) is unstable in the case of persistence, and also find conditions on parameter values which imply that the signs of the real parts of all the zeros of the characteristic function linearized at (v∗,w∗) are negative. Interestingly, one such condition is m≥p+k, pitting the effectiveness of infected cell removal versus virus removal. To obtain these results, we first employ an ad hoc approach, and then try, alternatively, a very nice theorem, with the latter yielding a somewhat stronger result.
There are some reasons for looking at this rather simple model despite the fact that it affords a less rich description of the biological process than many of the others and, even so, the results obtained are suboptimal. The methods of proof are not standard and may have applications elsewhere. In particular, the "integrally bounded" notion of Definition 3.3 has allowed us to show that the density of virus, respectively, that of infected cells either converges to or oscillates through v∗, respectively, w∗ when the density of virus or infected cells remains bounded above zero. Precisely because of its simplicity, expressions and equations are generally less numerous, shorter, and less opaque than for other models. The inequalities in the last paragraph suggest how greatly persistence might be affected if one could, by treatment or vaccine, modify particular parameter values of the system.
A number of questions have been left undecided. Is there a counterpart to Corollary 3.15 telling us that if V converges, then W does as well? Is there a relation between the "periods" of W−w∗ and V−v∗ when either has arbitrarily large zeros? Can we do better with regard to stability analysis? What happens when c+lq>a>c?
Thanks are given to the academic editor and reviewers for helpful suggestions.
The author declares that there is no conflict of interest.
[1] | G. I. Marchuk, Mathematical Models in Immunology, Optimization Software, Inc., Publications Division, Springer-Verlag, New York, 1983. |
[2] |
V. Bernhauerová, B. Lisowski, V. Rezelj, M. Vignuzzi, Mathematical modelling of SARS-CoV-2 infection of human and animal host cells reveals differences in the infection rates and delays in viral particle production by infected cell, J. Theor. Biol. , 531 (2021), 110895. https://doi.org/10.1016/j.jtbi.2021.110895 doi: 10.1016/j.jtbi.2021.110895
![]() |
[3] |
F. A. Rihan, V. Gandhi, Dynamics and sensitivity of fractional-order delay differential model for Coronavirus (COVID-19) infection, Prog. Fract. Diff. Appl., 7 (2021), 43–61. https://doi.org/10.18576/pfda/070105 doi: 10.18576/pfda/070105
![]() |
[4] |
L. N. Cooper, Theory of an immune system retrovirus, PNAS, 83 (1986), 9159–9163. https://doi.org/10.1073/pnas.83.23.9159 doi: 10.1073/pnas.83.23.9159
![]() |
[5] | N. Intrator, G. P. Deocampo, L. N. Cooper, Analysis of immune system retrovirus equations, in Theoretical Immunology, Part 2, A. S. Perelson, Ed., Addison-Wesley, Redwood City, Calif., (1988), 85–100. |
[6] | A. McLean, HIV infection from an ecological viewpoint, in Theoretical Immunology, Part 2, A. S. Perelson, Ed., Addison-Wesley, Redwood City, Calif., (1988), 77–84. |
[7] | S. Merrill, AIDS: Background and the dynamics of the decline of immunocompetence, in Theoretical Immunology, Part 2, A. S. Perelson, Ed., Addison-Wesley, Redwood City, Calif., (1988), 59–75. |
[8] | A. S. Perelson, Modeling the interaction of the immune system with HIV, in Mathematical and Statistical Approaches to AIDS Epidemiology (Lecture Notes Biomath., Vol. 83), Ed., C. Castillo-Chavez, Springer-Verlag, New York, (1989), 350–370. https://doi.org/10.1007/978-3-642-93454-4_17 |
[9] |
A. S. Perelson, D. E. Kirschner, R. De Boer, Dynamics of HIV infection of CD4+ T cells, Math. Biosci. , 114 (1993), 81–125. https://doi.org/10.1016/0025-5564(93)90043-A doi: 10.1016/0025-5564(93)90043-A
![]() |
[10] |
A. S. Perelson, R. M. Ribiero, Modeling the within host dynamics of HIV infection, BMC Biol. , 11 (2013), 96. https://doi.org/10.1186/1741-7007-11-96 doi: 10.1186/1741-7007-11-96
![]() |
[11] |
C. Rajivganthi, F. A. Rihan, Global dynamics of a stochastic viral infection model with latently infected cells, Appl. Sci. , 11 (2021), 10484. https://doi.org/10.3390/app112110484 doi: 10.3390/app112110484
![]() |
[12] |
X. Zhou, X. Song, X. Shi, A differential equation model of HIV infection of CD4+ T-cells with cure rate, J. Math. Anal. App. , 342 (2008), 1342–1355. https://doi.org/10.1016/j.jmaa.2008.01.008 doi: 10.1016/j.jmaa.2008.01.008
![]() |
[13] |
A. V. M. Herz, S. Bonhoeffer, R. M. Anderson, R. M. May, M. A. Novak, Viral dynamics in vivo: limitations on estimates of intracellular delay and virus decay, PNAS, 93 (1996), 7247–7251. https://doi.org/10.1073/pnas.93.14.7247 doi: 10.1073/pnas.93.14.7247
![]() |
[14] |
D. Li, W. Ma, Asymptotic properties of a HIV-1 infection model with time delay, J. Math. Ana. App., 335 (2007), 683–691. https://doi.org/10.1016/j.jmaa.2007.02.006 doi: 10.1016/j.jmaa.2007.02.006
![]() |
[15] |
P. W. Nelson, A. S. Perelson, Mathematical analysis of delay differential equation models of HIV-1 infection, Math. Biosci. , 179 (2002), 73–94. https://doi.org/10.1016/S0025-5564(02)00099-8 doi: 10.1016/S0025-5564(02)00099-8
![]() |
[16] |
J. Yang, X. Wang, F. Zhang, A differential equation model of HIV infection of CD4+ T-cells with delay, Discr. Dyn. Nat. Soc. , 2008 (2008). https://doi.org/10.1155/2008/903678 doi: 10.1155/2008/903678
![]() |
[17] |
K. Hattaf, N. Yousfi, Qualitative analysis of a generalized virus dynamics model with both modes of transmission and distributed delays, Int. J. Differ. Equations, 2018 (2018). https://doi.org/10.1155/2018/9818372 doi: 10.1155/2018/9818372
![]() |
[18] |
M. Maziane, K. Hattaf, N. Yousfi, Spatiotemporal dynamics of an HIV infection model with delay in immune response activation, Int. J. Differ. Equations, 2018 (2018). https://doi.org/10.1155/2018/3294268 doi: 10.1155/2018/3294268
![]() |
[19] | D. Kirschner, Using mathematics to understand HIV immune dynamics, Not. AMS, 43 (1996), 191–202. |
[20] | R. Bellman, K. L. Cooke, Differential-Difference Equations, Academic Press, New York, 1963. https://doi.org/10.1063/1.3050672 |
[21] | K. L. Cooke, P. van den Driessche, On zeroes of some transcendental equations, Funkcialaj Ekvacioj, 29 (1986), 77–90. |