A two-dimensional discrete delay-differential system model of viremia

1.   Introduction

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.

2.   The model and bounds on its solution

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 ${\nu }_{0}$ is introduced into the formerly virus-free blood.

We denote the rate at which virions penetrate attacked cells by $p$, where $0 < p\le 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\ge 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 $\eta \left(t\right)$ and another, from circulating antibodies and other immune mechanisms, with relative rate $\kappa \left(t\right)$. For the infected cells, we assume a simple clearance independent of host cells at a relative rate $\mu \left(t\right)$. These functions, $\eta$, $\kappa$, and $\mu$ defined on ($-T, \infty$), 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\left(t\right)$ be the density of free circulating virions at time $t$ (equivalently, the number of virions per host cell) and $W\left(t\right)$ the density of those host cells which have been infected (warped?). We define $E$ by $E\left(t\right)=pV\left(t\right)(1-W\left(t\right))$, the rate at which virus infects previously uninfected host cells. Virions are cleared by intact host cells at the rate $\eta V\left(1-W\right)=\frac{\eta }{p}E$, and by host cell independent mechanisms at the rate $\kappa 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.

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 $\; {\mathrm{\Sigma }}\;$ be the closed strip in two-dimensional phase space defined by $V\left(t\right)\ge 0$ and $0\le W\left(t\right)\le 1$. The following proposition establishes bounds on the region in which the trajectory of $(V\left(t\right), W\left(t\right))$ lies.

Proposition 2.1: Let $\gamma =lT+b{e}^{T}$. Then, for all $t > 0$,

(i) $(V\left(t\right), W\left(t\right))\in {\; {\mathrm{\Sigma }}\; }^{°}$, the interior of $\; {\mathrm{\Sigma }}\;$.

(ii) $V\left(t\right)\le {v}_{0}+\gamma +\; {\mathrm{m}}\; \; {\mathrm{a}}\; \; {\mathrm{x}}\; \{\gamma, \frac{\gamma }{pT}\}.$

(iii) $W\left(t\right)\le 1-\frac{{e}^{-p{v}_{S}T}}{1+p{v}_{S}T}$, where ${v}_{S}=\; {\mathrm{sup}}\; V\left(t\right)$ (which exists by ii).

Proof: (i) The initial conditions (2.1c) state that $V\left(0\right) > 0$ and $W\left(0\right)=0$. By (2.1b), we also have ${W}^{\text{'}}\left(0\right)=E\left(0\right)=pV\left(0\right) > 0$. Therefore, by continuity and these inequalities, $(V\left(t\right), W\left(t\right))\in {\; {\mathrm{\Sigma }}\; }^{°}$ on an interval $(0, u)$, for some $u > 0$. Suppose $u$ is maximal with this property implying, by continuity and the fact that ${\; {\mathrm{\Sigma }}\; }^{°}$ is open, that $(V\left(u\right), W\left(u\right))$ lies on the boundary of $\; {\mathrm{\Sigma }}\;$. Let ${\nu }_{u}=\underset{t\in [0, u]}{\; {\mathrm{max}}\; }V\left(t\right)$. Then, by (2.1b), for $t\in [0, u)$,

Integrating from $0$ to $u$, using continuity of $W$ for positive $t$, and then exponentiating, we get

Next, for $t\in (0, u)$, by (2.1a), the definition of $E\left(t\right)$, and what we have just shown, and letting $H$ and $K$ be upper bounds for $\eta$ and $\kappa$ respectively,

Integrating from $0$ to $u$, and then exponentiating as above, we get

Finally, we take a derivative and use (2.1b) again, yielding for $t > 0$,

If we let $g\left(t\right)$ be the first term on the right-hand side of the last equality (so the second is $g(t-T)$), note that $g\left(t\right)\equiv 0$ on $(-T, 0]$, and integrate from $0$ to $t$, using simple substitution on the right, we get

and, since $E > 0$ on $(0, u)$, $W\left(u\right) > 0$. Therefore, $(V\left(u\right), W\left(u\right))\in {\; {\mathrm{\Sigma }}\; }^{°}$, a contradiction, so $(V\left(t\right), W\left(t\right))\in {\; {\mathrm{\Sigma }}\; }^{°}$ for all $t > 0$.

(ii) Next, we construct a special sequence of times, $0={t}_{0} < {t}_{1} < \dots < {t}_{n} < \dots$ which increases without bound. Assuming we have determined ${t}_{n}$, if $V\left({t}_{n}\right)\le V\left(t\right)$ for $t\in [{t}_{n}, {t}_{n}+T]$, take ${t}_{n+1}={t}_{n}+T$ and, otherwise, let ${t}_{n+1}$ be any value in $[{t}_{n}, {t}_{n}+T]$ for which $V\left(t\right)$ attains its minimum on that closed interval. In either case ${t}_{n+1}-{t}_{n}\le T$ and ${t}_{n+2}\ge {t}_{n}+T$ so, indeed, ${t}_{n}\to \infty$. We first wish to show by induction that for all $n$, $V\left({t}_{n}\right) < {v}_{0}+\; {\mathrm{m}}\; \; {\mathrm{a}}\; \; {\mathrm{x}}\; \{\gamma, \frac{\gamma }{pT}\}$. Since $V\left({t}_{0}\right)=V\left(0\right)={v}_{0}$, this is true for $n=0$, so suppose it true for some $n\ge 0$. If ${t}_{n+1}$ was determined as a minimum for $V\left(t\right)$ on [${t}_{n}, {t}_{n}+T]$, then $V\left({t}_{n+1}\right) < V\left({t}_{n}\right)$, so the statement is true for $n+1$. Otherwise, ${t}_{n+1}={t}_{n}+T$ and ${t}_{n}$ is a minimum for $V\left(t\right)$ on [${t}_{n}, {t}_{n}+T]$. The following inequalities are consequences of (2.1a) and (2.2) and the facts that $W\left(t\right) < 1$ and $\mu \left(t\right) < 1$ for all $t$.

Integrating (2.3) from ${t}_{n}$ to ${t}_{n+1}$ and using (2.4), we get

If $pT\ge 1$, $V\left({t}_{n+1}\right)\le \gamma$, and the statement is true for $n+1$. If $pT < 1$, then, $\gamma < \frac{\gamma }{pT}$, and by the induction assumption,

and, again, the statement is true. Now any positive $t\in$ [${t}_{n}, {t}_{n+1}$] for some $n\ge 0$. Since ${t}_{n+1}-{t}_{n}\le T$, weakening (2.3) above to ${V}^{\text{'}}\left(t\right)\le l+bE(t-T)$, integrating from ${t}_{n}$ to $t$, and then using (2.4), we get

which finishes the proof of ii.

(iii) First, we claim that for any positive $t$, there is an ${t}_{0}\in [t-T, t]$ such that $W\left({t}_{0}\right)\le \frac{p{v}_{S}T}{1+p{v}_{S}T}$. Suppose this were not true. Then, we would have for all $x\in [t-T, t]$,

so, by (2.2) and the fact that $t\in [t-T, t]$, implying that $W\left(t\right) > \frac{p{v}_{S}T}{1+p{v}_{S}T}$,

a contradiction, so there is, indeed, an ${t}_{0}\in [t-T, t]$ such that $1-W\left({t}_{0}\right)\ge \frac{1}{1+p{v}_{S}T}$. Now, integrating the inequality $\frac{-{W}^{\text{'}}\left(t\right)}{1-W\left(t\right)}\ge -p{\nu }_{S}$, which follows from (2.1b) as in the proof of part i, this time from ${t}_{0}$ to $t$, and exponentiating, we have

Proposition 2.1 places $V$ and $W$ and, because they satisfy Eqs (2.1), their derivatives as well, in the real algebra of bounded C¹ 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 $\left|f\left(t\right)-f\left(u\right)\right|\le M\left|t-u\right|$).

3.   Ultimate behavior of the solution

Let $\tau$ be so large that for $t\ge \tau -T$, $\mu \left(t\right)\equiv m$, $\eta \left(t\right)\equiv h$, and $\kappa \left(t\right)\equiv k$, making the system autonomous with constant coefficients. We also assume that $\tau > 2T$ to guarantee that $V$ and $W$ are of class C² for $t\ge \tau$. It will be convenient to define $a=p+k > 0$, and $c=pb{e}^{-mT}-h$. Also, let $F\left(t\right)={e}^{-mT}(E\left(t\right)-E\left(t-T\right))$, $q=\frac{p(1-{e}^{-mT})}{m}$ for $m > 0$, and $q=pT$ if $m=0$. (Interestingly, $\underset{m\to 0}{\; {\mathrm{lim}}\; }\frac{1-{e}^{-mT}}{m}=T$, by L'Hôspital's Rule and, as a function of $m$, $q$ is decreasing, so $q\le pT$). Let $J\left(t\right)=qV\left(t\right)(1-W\left(t\right))-W\left(t\right)$. Then, for $t\ge \tau$, the system becomes

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}^{\text{'}}\left(t\right)=E\left(t\right)-E\left(t-T\right)=F\left(t\right)$. The purpose of the term $J=qV\left(1-W\right)-W$ in the final forms of defining Eqs (3.1) will become apparent shortly. We can also rewrite (2.2) as

By a fixed point, $({v}^{*}, {w}^{*})$, of the system, we mean a constant solution $V\left(t\right)\equiv {v}^{*}$, $W\left(t\right)\equiv {w}^{*}$ to Eqs (3.1) on $(-T, \infty)$. We are interested in them because any point in $\; {\mathrm{\Sigma }}\;$ to which $(V\left(t\right), W\left(t\right))$ might converge must be a fixed point.

Proposition 3.1: Fixed points for the system include $\left(0, 0\right)$ and, if and only if $c+lq > a$, $\left({v}^{*}, {w}^{*}\right)=(\frac{c+lq-a}{qa}, \frac{c+lq-a}{c+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)\in \; {\mathrm{\Sigma }}\;$ is a fixed point of the autonomous system if and only if

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}^{*})\in \; {\mathrm{\Sigma }}\;$. Then by the first, respectively second, equation,

Substituting (3.3b) into the first form of (3.3a), we get $\left(c+lq\right)\left(1-{w}^{*}\right)=a$ so, since ${w}^{*} > 0$, a second fixed point exists if and only if $c+lq > a$ and, in that case,

From now on, whenever $c+lq > a$, $\left({v}^{*}, {w}^{*}\right)$ will denote the nonzero fixed point given by this proposition.

Lemma 3.2: There are positive constants ${d}_{1}$ and ${d}_{2}$ such that for all $t\ge \tau$,

Proof: By the first form of (3.1a), and the fact that $W\le 1$, ${V}^{\text{'}}\left(t\right)\ge -(h+a)V\left(t\right)$. For $x\in [t-T, t]$, integrating from $x$ to $t$ and exponentiating, we get

Therefore, by (3.2), we have

so we can take ${d}_{2}=pT{e}^{(h+a)T}$. To prove the left-hand inequality, let $d$ be the larger of $\frac{{v}_{S}}{T}$ and a Lipschitz constant for $V$, and observe that for all $t\ge \tau$, the line passing through $(t, V\left(t\right))$ with slope $d$ must cross the $t$-axis at some $u\ge t-\frac{V\left(t\right)}{{v}_{S}}T$ and also lie below the plot of the function $V$ on $[t-T, t]$. Therefore, ${\int }_{t-T}^{t}V$ exceeds the area of the triangle with height $V\left(t\right)$ and base $\frac{V\left(t\right)}{d}$, i.e., $\frac{{\left(V\left(t\right)\right)}^{2}}{2d}$. Again, employing (3.2) along with Proposition 2.1iii, we have

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 [$\tau, \infty)$. We shall say such a function, $f$, is integrally bounded if there is a ${C}_{f} > 0$ such that $\left|{\int }_{u}^{t}f\right|\le {C}_{f}$ for all $\tau \le u < t$ and, in that case, we write $f\sim 0$, a definition relevant only for this paper. For example, $f\left(t\right)=\; {\mathrm{sin}}\; t\sim 0$ and its terminal behavior is cyclic, while $f\left(t\right)=\frac{1}{t}$ even converges to zero, but is not integrally bounded. It is straightforward to demonstrate that if $f\sim 0$ and $g\sim 0$, then $\alpha f+\beta g\sim 0$ for any $\alpha, \beta \in \mathbb{R}$. Even further, if $f\sim 0$ and $g$ is any bounded function, then $gf\sim 0$, and if $g$ is bounded away from zero, then $\frac{f}{g}\sim 0$. The derivative of any bounded function is integrally bounded, as we have ${\int }_{u}^{t}{f}^{\text{'}}=f\left(t\right)-f\left(u\right)$. We shall also say that $f \sim g$ if $f-g\sim 0$, so that if $f \sim g$ and $g\sim h$, then $f\sim h$.

Lemma 3.4: If $f: [\tau, \infty)\to \mathbb{R}$ is bounded, then $f\left(t\right) \sim f\left(t-T\right)$. In particular, $F\left(t\right)={e}^{-mT}\left(E\left(t\right)-E\left(t-T\right)\right)\sim 0$.

Proof: For any $\tau \le u < t$, we have

Lemma 3.5: (Barbãlat's Lemma): Let ${t}_{0}\ge \tau$, $g: [{t}_{0}, \infty)\to \text{R}$ be differentiable with uniformly continuous derivative. Then, if $\underset{t\to \infty }{\; {\mathrm{lim}}\; }g\left(t\right)$ exists (and is finite), $\underset{t\to \infty }{\; {\mathrm{lim}}\; }{g}^{\text{'}}\left(t\right)=0$.

Proof (included for completeness): Suppose ${g}^{\text{'}}\left(t\right)$ does not converge to $0$ as $t\to \infty$, so that there is an $\epsilon > 0$ such that $\left|{g}^{\text{'}}\left(t\right)\right|\ge \epsilon$ for arbitrarily large values of $t$. Since ${g}^{\text{'}}$ is uniformly continuous, there is a $\delta > 0$ such that $\left|{g}^{\text{'}}\right|\ge \frac{ \epsilon }{2}$ on the interval $[t, t+\delta ]$ for all such $t$ and, consequently, $\left|g\left(t\right)-g(t+\delta)\right|\ge \frac{ \epsilon \delta }{2}$, which precludes the convergence of $g\left(t\right)$.

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 $\sim$.

Corollary 3.6: Let ${t}_{0}\ge \tau$ and $g: \left[\tau, \infty \right)\to \mathbb{R}$ be Lipschitz and $g \sim 0$. Then either $\underset{t\to \infty }{\; {\mathrm{lim}}\; }g\left(t\right)=0$ or $g$ changes sign for arbitrarily large values of $t$.

Proof: Suppose that ${t}_{0}\ge \tau$ is such that $g$ takes only nonnegative or nonpositive values on $\left[{t}_{0}, \infty \right)$ and let $G: \left[{t}_{0}, \infty \right)\to \mathbb{R}$ be defined by $G\left(t\right)={\int }_{{t}_{0}}^{t}g$. Then $G$ is monotone and, since $g \sim 0$, it is bounded, so $G$ converges. But, ${G}^{\text{'}}\left(t\right)=g\left(t\right)$ and, since $g$ is Lipschitz, it is uniformly continuous. Lemma 3.5 applies.

Lemma 3.7: Let $f: [\tau, \mathbb{\infty })\to \mathbb{R}$ be continuous and bounded. Then, ${\int }_{x-T}^{x}f \sim 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\ge \tau$ and consider ${\int }_{u}^{t}{\int }_{x-T}^{x}f\left(y\right)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 $\left(u-T, u\right)$ 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

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 \sim 0$

But, if $m=0$, by (3.2) and Lemma 3.7, $W=p{\int }_{t-T}^{t}V\left(1-W\right) \sim pTV\left(1-W\right)=qV\left(1-W\right)$, 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\left)W\right) \sim 0$

Definition 3.9: When $c+lq > a$, define, for all $t\ge \tau,$ $A\left(t\right)=V\left(t\right)({w}^{*}-W\left(t\right))$, $B\left(t\right)=W\left(t\right)({v}^{*}-V(t\left)\right)$, and $D\left(t\right)=\frac{{w}^{*}}{{v}^{*}}V\left(t\right)-W\left(t\right)$. Notice that $A-B={v}^{*}D$.

Lemma 3.10: If $c+lq > a$, then

(i) $A \sim 0$ and $B \sim 0$

(ii) ${V}^{\text{'}}\left(t\right)=\frac{1}{{v}^{*}}{\left(\right(cv}^{*}-l)A+lB)-bF\left(t\right)$

Proof: (i) Remark 3.8ii and the second form of (3.3c) tell us that $A \sim 0$. By Remark 3.8i and (3.3b), we have

so $D \sim 0$. But $B=A-{v}^{*}D$, so $B \sim 0$.

(ii) We have, by the second form of (3.1a), and the first form of (3.3a),

But, we also have

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\le a$, then $(V\left(t\right), W\left(t\right))\to \left(\; {\mathrm{0, 0}}\; \right)$ as $t\to \infty$.

(ii) If $c+lq > a$ and $c\ge a$ (either $c > a$ or $c=a$ and $l > 0)$, then $V\left(t\right)$ and $W\left(t\right)$ are bounded above zero.

Proof: (i) Suppose that $0 < c+lq\le a$. By the right-hand side of Lemma 3.2, it is adequate to show that $V\left(t\right)\to 0$. By Remark 3.8ii, $V(\frac{c+lq-a}{c+lq}-W) \sim 0$. If $c+lq < a$, then $\frac{c+lq-a}{c+lq}-W$ is bounded below zero which, by Definition 3.3, implies that $V \sim 0$ and, by Corollary 3.6, $V\left(t\right)\to 0$. If $c+lq=a$, then $VW \sim 0$ and, this time, $V\left(t\right)W\left(t\right)\to 0$. But by the left-hand side of Lemma 3.2, $V\left(t\right)W\left(t\right)\ge {d}_{1}{\left(V\left(t\right)\right)}^{3}$, so $V\left(t\right)\to 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, {${t}_{n}\}$, $n\ge 1$, with $W\left({t}_{n}\right)\to 0$. Therefore, we can find a ${t}_{0}\ge \tau$ such that $W\left({t}_{0}\right)\le \; {\mathrm{min}}\; \{{w}^{*}, {d}_{1}{\left({v}^{*}\right)}^{2}\}$ and, by replacing {${t}_{n}\}$ with a terminal subsequence if necessary, assume that ${t}_{n} > {t}_{0}$ and $W\left({t}_{n}\right) < {e}^{-mT}W\left({t}_{0}\right)$ for all $n\ge 1$. Now define a second sequence {${u}_{n}\}$, $n\ge 1$, by letting ${u}_{n}$ be the largest number less than ${t}_{n}$ such that $W\left({u}_{n}\right)=W\left({t}_{0}\right)$. (It certainly may be the case that {${u}_{n}\}$ is not strictly increasing and it could even have a maximal element.) Using (3.2), we derive the string of inequalities

by which, looking back at the proof of Lemma 3.4, we can see that ${\int }_{{u}_{n}}^{{t}_{n}}F\le 0$. Now, by the definition of ${u}_{n}$, on the interval [${u}_{n}, {t}_{n}$], $W\le {w}^{*}$ and, by Lemma 3.3, ${d}_{1}{V}^{2}\le W\le {d}_{1}{\left({v}^{*}\right)}^{2}$, so $V\le {v}^{*}$. Recall that Lemma 3.10ii states

Integrating this equation from ${u}_{n}$ to ${t}_{n}$ gives us

Note that by our assumption that $c\ge a$ and the second form of (3.3a), ${cv}^{*}-l\ge 0$. By our definition of {${u}_{n}\}$ then, the right-hand side is nonnegative, so $V\left({u}_{n}\right)\le V\left({t}_{n}\right)$. But by Lemma 3.2, $V\left({t}_{n}\right)\to 0$, so, $V\left({u}_{n}\right)\to 0$ and, by the same lemma, $\frac{W\left({t}_{0}\right)}{{d}_{2}}=\frac{W\left({u}_{n}\right)}{{d}_{2}}\le V\left({u}_{n}\right)$, 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\left(b{e}^{-mT}-1\right)-(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\ge a$. Then

(i) ${w}^{*}-W\sim 0$ and ${v}^{*}-V\sim 0$.

(ii) If ${v}^{*}-V$, respectively ${w}^{*}-W$, does not change sign at arbitrarily large values of $t$, then $V\left(t\right)\to {v}^{*}$, respectively $W\left(t\right)\to {w}^{*}$.

Proof: (i) We already know, by Lemma 3.10i, that ${V(w}^{*}-W) \sim 0$ and $W({v}^{*}-V) \sim 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\left(t\right)$ passes back and forth through ${v}^{*}$ or $W\left(t\right)$ 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 $\left({v}^{*}, {w}^{*}\right)$. 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=[u}_{0}, {t}_{0}]$ where $\tau \le {u}_{0}$. We'll be manipulating rational functions of $V\left(t\right), V\left(t-T\right), W\left(t\right), W\left(t-T\right)$ with denominators bounded away from zero, all of which are bounded Lipschitz functions with bounded Lipschitz derivatives.

Lemma 3.13: Let $\epsilon > 0$ be given and suppose that for some $C=C\left(\epsilon \right) > 0$, $f: {I=[u}_{0}, {t}_{0}]\to \mathbb{R}$ satisfies $\left|f\left({t}_{0}\right)\right| < \epsilon$ and ${f}^{\text{'}}\left(t\right)=-C(f\left(t\right)-\varepsilon (t\left)\right)$ on $I$, where $\left|\varepsilon \left(t\right)\right| < \epsilon$. Then $\left|f\left(t\right)\right| < \epsilon$ on $I$.

Proof: Suppose, to the contrary, that $f\left(t\right)= \epsilon$ for some $t\in I$, where $t$ is as small as possible. Then $t\ne {u}_{0}$ and, since $f < \epsilon$ on the interval [${u}_{0}, t$), it must be the case that ${f}^{\text{'}}\left(t\right)\ge 0$. However, this is contradicted by

Similarly, $f\left(t\right) > - \epsilon$ for $t\in I$.

Theorem 3.14: If $c+lq > a$, $c\ge a$, and there are arbitrarily long intervals on which ${w}^{*}-W$ doesn't change sign, then there are arbitrarily long intervals on which $(V\left(t\right), W\left(t\right))$ is arbitrarily close to $({v}^{*}, {w}^{*})$.

Proof: First, we claim that there are arbitrarily long intervals on which $W\left(t\right)$ is arbitrarily close to ${w}^{*}$ and $W\text{'}\left(t\right)$ is arbitrarily close to zero as well. Assume that there are arbitrarily long intervals on which $W\ge {w}^{*}$. The proof is similar for $W\le {w}^{*}$ and one or the other must be true for infinitely many such intervals. Let $\epsilon > 0$ be given, let $M$ be a Lipschitz constant for $W$, and let $\; {\mathrm{\Lambda }}\;$ be the length of an interval on which we wish to establish that $W-{w}^{*} < \epsilon$. Generality is not lost if we assume that $\; {\mathrm{\Lambda }}\; \ge \frac{ \epsilon }{M}$, so consider such an interval, $I$. Reasoning as in the proof of Lemma 3.2, for any $t\in I$, if $W\left(t\right)-{w}^{*}= \epsilon$, then the right triangle with height $\epsilon$ and base either [$t, t+\frac{ \epsilon }{M}$] or [$t-\frac{ \epsilon }{M}, t$], whose area is $\frac{{ \epsilon }^{2}}{2M}$ lies beneath the plot of $W\left(t\right)-{w}^{*}$. Now, by Corollary 3.12i, there is a $C > 0$ such that $\left|{\int }_{u}^{t}(W-{w}^{*})\right|\le C$ for all $\tau \le u < t$. Let the integer $n$ be so large that $\frac{n{ \epsilon }^{2}}{2M} > C$ and let $\tau \le {u}_{0} < {t}_{0}$ such that $W\left(t\right)-{w}^{*}\ge 0$ on [${u}_{0}, {t}_{0}$] and ${t}_{0}-{u}_{0}=n\; {\mathrm{\Lambda }}\;$. Suppose that $W\left(t\right)-{w}^{*}$ took the value $\epsilon$ on all of the subintervals [${u}_{0}+j\; {\mathrm{\Lambda }}\; ,{u}_{0}+\left(j+1\right)\; {\mathrm{\Lambda }}\;$], $j=0, \dots, n-1$. Then we would have the contradiction

so, it must be the case that $W\left(t\right)-{w}^{*} < \epsilon$ on at least one of these subintervals of length $\; {\mathrm{\Lambda }}\;$.

Now, let $\epsilon > 0$ be given and $M$ be a Lipschitz constant for $W\text{'}$. By what has been done above, we can find an interval $I$ of length at least $\frac{ \epsilon }{M}$ and arbitrarily large on which $\left|{w}^{*}-W\right| < \; {\mathrm{min}}\; \{ \epsilon, \frac{{ \epsilon }^{2}}{8M}\}$, immediately implying that $\left|W\left(t\right)-W\left(u\right)\right| < \frac{{ \epsilon }^{2}}{4M}$ for any $t, u\in I$. Suppose first that $W\text{'}\left({t}\right)\ge \epsilon$ for some ${t} \in I$, Then at least one of ${t} \pm \frac{ \epsilon }{2M}\in I$, say ${t}+\frac{ \epsilon }{2M}\in I$. By the definition of $M$, $W\text{'} \ge \frac{ \epsilon }{2}$ on [${t}, {t}+\frac{ \epsilon }{2M}]$, so $W({t}+\frac{ \epsilon }{2M})-W\left({t}\right)\ge \frac{{ \epsilon }^{2}}{4M}$, a contradiction. Of course, the reasoning is analogous for ${W}^{\text{'}}\le - \epsilon$ and/or ${t}-\frac{ \epsilon }{2M}\in 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 $\epsilon > 0$ and $\; {\mathrm{\Lambda }}\; > 0$ be given. Let the integer $n$ be so large that ${e}^{-nmT}{v}^{*} < \epsilon$. The first form of (3.1b) can be written

Now, by the first form of (3.3d) and the definition of $q$, $\frac{m{w}^{*}}{p\left(1-{w}^{*}\right)}=(1-{e}^{-mT}){v}^{*}$ and

Therefore, looking at the last terms of the two equations above, we see that there is a $C > 0$ which does not depend on $\epsilon$ such that if $I$ is an interval on which $\left|{w}^{*}-W\left(t\right)\right| < \epsilon$, $\left|{W}^{\text{'}}\left(t\right)\right| < \epsilon$, and $t-T, t\in I$, then

where $\left|\varepsilon \left(t\right)\right| < C \varepsilon$. By what we have done above, there is, in fact, such an interval, $\text{I}$, of length $\; {\mathrm{\Lambda }}\; +nT$. Let $t$ be in the terminal subinterval of $\text{I}$ with length $\; {\mathrm{\Lambda }}\;$. Then, we have a telescoping sum and simplifications

Next suppose that $m=0$, $\epsilon > 0$ and $\; {\mathrm{\Lambda }}\; > 0$ have been given, and let $I$ be an interval of length $\; {\mathrm{\Lambda }}\;$ on which $|{w}^{*}-W\left(t\right)| < \epsilon$ and $\left|{W}^{\text{'}}\left(t\right)\right| < \epsilon$. The first form of (3.1b) is ${W}^{\text{'}}\left(t\right)=E\left(t\right)-E\left(t-T\right)=F\left(t\right)$, so $\left|F\left(t\right)\right| < \epsilon$, as well. By (3.2), we have

By the first form of (3.3d), ${v}^{*}T=\frac{{w}^{*}}{p\left(1-{w}^{*}\right)}$ so, using the last equation,

Now, suppose further that $l=0$. By Lemma 3.10ii, ${V}^{\text{'}}\left(t\right)=cA\left(t\right)-bF\left(t\right)$, so that $\left|{V}^{\text{'}}\left(t\right)\right| < (c{v}_{S}+b) \epsilon$. Therefore, for any $t-T, t\in I$, if $x\in [t-T, t]$, then $\left|V\left(t\right)-V\left(x\right)\right| < (c{v}_{S}+b)T \epsilon$. Integrating over $[t-T, t]$ and then using the last inequality just above, we have

Then, by this inequality and the one just above,

Finally, suppose that $m=0$ but $l > 0$, and $I$ is as above. By Lemma 3.10ii,

where $\left|\zeta \left(t\right)\right| < {v}_{S} \epsilon$, so

where $\left|\varepsilon \left(t\right)\right| < C \epsilon$ for some $C > 0$. Now, we could have picked $I$ such that $\left|{w}^{*}-W\left(t\right)\right| < \frac{l{w}^{*}}{{v}^{*}C} \epsilon$ and $\left|{W}^{\text{'}}\left(t\right)\right| < \frac{l{w}^{*}}{{v}^{*}C} \epsilon$ on $I$, and would have deduced the last equation above except with $\left|\varepsilon \left(t\right)\right| < \frac{l{w}^{*}}{{v}^{*}} \epsilon$. We now could apply Lemma 3.13 if we knew that $\left|V\left({u}_{0}\right)-{v}^{*}\right|$ was sufficiently small, where ${u}_{0}$ was the left endpoint of $I$. But by Corollary 3.12i, for any $\epsilon > 0$, $V\left(t\right)$ cannot remain at least $\epsilon$ distant from ${v}^{*}$ on $I$ if $I$ is long enough. Therefore, taking $I$ twice that long allows us to assume that $\left|V\left({t}_{0}\right)-{v}^{*}\right|$ is small enough.

Corollary 3.15: If $c+lq > a$, $c\ge a$, and $W$ converges, then $(V\left(t\right), W\left(t\right))\to ({v}^{*}, {w}^{*})$.

Proof: By Corollary 3.12ii, if $W$ converges, then $W\left(t\right)\to {w}^{*}$. Then, by Lemma 3.5 ${W}^{\text{'}}\left(t\right)\to 0$. This allows us to find for any $\epsilon > 0$ an interval $I=[{u}_{0}, \infty)$ on which both $\left|{w}^{*}-W\left(t\right)\right| < \epsilon$ and $\left|{W}^{\text{'}}\left(t\right)\right| < \epsilon$. Using the reasoning in the proof of the theorem, we conclude that $\left|V\left(t\right)-{v}^{*}\right| < \epsilon$ on $[{t}_{0}, \infty)$ for some ${t}_{0}\ge {u}_{0}$.

Corollary 3.16: If $c+lq > a$, $c\ge a$, and $({v}^{*}, {w}^{*})$ is locally asymptotically stable, then either $(V\left(t\right), W\left(t\right))\to ({v}^{*}, {w}^{*})$ or the distance between successive values of $t$ for which $W\left(t\right)={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.

4.   Cycles and the undelayed system

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\left(-T\right)=V\left(0\right)$ and $W\left(-T\right)=W\left(0\right)$. Our solution couldn't be one because $V\left(-T\right)=0\ne {v}_{0}=V\left(0\right)$ (even though $W\left(-T\right)=W\left(0\right)=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\left(t\right)-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,

obtained from Eqs (3.1) by removing $F\left(t\right)$ from the right-hand side (not removing only $E\left(t-T\right)$).

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}^{\; {\mathrm{\text{'}}\; }}$and ${W}^{\text{'}}$are as in the equations above, $\frac{\partial {V}^{\text{'}}}{\partial V}+\frac{\partial {W}^{\text{'}}}{\partial W}$ does not take both signs on ${\; {\mathrm{\Sigma }}\; }^{°}$. In fact,

Therefore, there are no cycles if $c\le a$. If $c+lq\le 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.

5.   Local considerations

In this section, we shall look at the stability of the fixed points $(0, 0)$ and, when $c+lq > a$ and $c\ge 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}^{\text{'}}$ and the last form of (3.1b) for ${W}^{\text{'}}$.

Observe that

Linearizing around $(0, 0)$ but leaving out the last terms and subtracting $\lambda {\bf{I}}$, we get the matrix

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 $\lambda {\bf{I}}$, we get the matrix

where the second equality follows from the first form of (3.3a), (3.3b), and the second form of (3.3d). Again, ${\bf{A}}$ is the matrix of the linearization of the undelayed system. Therefore,

Looking on the right and using equations (3.3d) and second form of (3.3c),

We need to add to ${\bf{A}}$ the matrix corresponding to the last terms of the defining equations above. We have observed above that $\frac{\partial }{\partial V}E=p(1-W)$ and $\frac{\partial }{\partial W}E=-pV$, so this matrix is $p{e}^{-mT}(1-{e}^{-\lambda T}){\bf{B}}$, where

Then, $G\left(\lambda \right)=\left|{\bf{A}}+p{e}^{-mT}(1-{e}^{-\lambda T}){\bf{B}}\right|$. To calculate $G\left(\lambda \right)$, we use the multilinearity on rows of the determinant, noting that $\left|{\bf{B}}\right|=0$

Note that the constant term is strictly positive as, either $\text{l} > 0$or (second form of (3.3a)) $c{v}^{*}-l=\frac{{w}^{*}}{{v}^{*}}\left(c-a\right) > 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

To make our characteristic function look less unwieldy, let

and note that $\delta =0$ if and only if $m=l=0$. We now rewrite the characteristic function as

Initially, we deal with $G$ in an ad hoc fashion. First, consider $G$ restricted to $[0, \infty)$. The first form of (5.1) demonstrates that $G$ is positive on $(0, \infty)$, and $G\left(0\right)=maq{v}^{*}$, so $G\left(0\right)=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\left(\stackrel{-}{\lambda }\right)=\stackrel{-}{G\left(\lambda \right)}$, so we may restrict our search for other roots to those with positive imaginary parts. Suppose that $\lambda$ is a zero of $G$, where $\lambda =x+iy$, $x\ge 0$ and $y > 0$. The following equations represent the real part of $G\left(\lambda \right)=0$ and the imaginary part divided by $y$.

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\left(\lambda \right)=0$, where $\lambda =x+iy$, $x\ge 0$ and $y > 0$. Let $f\left({w}^{*}\right)=\frac{3\pi }{(3\pi -1)(1-{w}^{*})}+\frac{3\pi }{{(3\pi -1)}^{2}{\left(1-{w}^{*}\right)}^{2}}+\frac{{w}^{*}}{{\left(1-{w}^{*}\right)}^{3}}$. Then

(i) $x\left(2yT-1\right) < \beta$

(ii) $\delta < \beta < a; \; \; \; \; \; \; \; \; \; \; m < c{\left(1-{w}^{*}\right)}^{2} < \alpha (1-{w}^{*})$

(iii) $\frac{3\pi }{2} < yT; \; \; \; \; \; \; \; \; \; \; x < \frac{\beta }{3\pi -1} < \frac{a}{3\pi -1}$

(iv) $T > \frac{3\pi }{2\alpha (1+f({w}^{*}\left)\right)}$

Proof: (i) Since $\text{m} > 0$, $\delta > 0$ so, for (5.2b) to be true, the last term on its left-hand side must be negative. Therefore,

It is a straightforward exercise in elementary calculus to show that for $C\in \mathbb{R}$, the function $\; {\mathrm{cos}}\; z+C\; {\mathrm{sin}}\; z$ is maximized when $z={\; {\mathrm{tan}}\; }^{-1}C$ and that the maximum is $\sqrt{1+{C}^{2}}$. Therefore, it must be the case that

(ii) Again, because the last term on the left-hand side of (5.2b) is negative, $\; {\mathrm{sin}}\; yT < 0$. Therefore, $\left|\; {\mathrm{cos}}\; yT\right| < 1$ and $yT > \pi$. To save space, let $\zeta =\alpha (1-{e}^{-\lambda T})$ and $\xi =\delta +\zeta$. Observe, since $1-{e}^{-\lambda T}=1-{e}^{-xT}\; {\mathrm{cos}}\; yT+i{e}^{-xT}\; {\mathrm{sin}}\; yT$, that $\zeta$ and, therefore, $\xi$ have positive real and negative imaginary parts. Rewrite $G$ as a quadratic function of $\lambda$ and apply the quadratic formula to get

where $\sigma$ is the correct one of the two possible choices for the square root. Since $\lambda$ has nonnegative real part and positive imaginary part, $0 < {\mathrm{Re}}\xi \le {\mathrm{Re}}\sigma$ and ${\mathrm{Im}}\xi < {\mathrm{Im}}\sigma$. Since the real part of $\zeta$ is positive, ${\mathrm{Re}}{\sigma }^{2}={{\mathrm{Re}}(\xi }^{2}-4\left(maq{v}^{*}+\zeta \beta \right)) < {\mathrm{Re}}{\xi }^{2}$. Therefore,

Now, should ${\mathrm{Im}}\sigma < 0$, along with ${\mathrm{Im}}\xi < 0$ and ${\mathrm{Im}}\sigma > {\mathrm{Im}}\xi$, then $\left|{\mathrm{Im}}\sigma \right| < \left|{\mathrm{Im}}\xi \right|$, contradicting the inequality immediately above. Therefore, ${\mathrm{Im}}\sigma > 0$ and so ${\mathrm{Im}}\sigma > -{\mathrm{Im}}\xi$, yielding

However, from its definition above, $\beta \le \frac{mb\left(1-{w}^{*}\right)+a{v}^{*}}{b\left(1-{w}^{*}\right)+{v}^{*}}$, a weighted average of $a$ and $m$ with strictly positive weights. Therefore, $a > \beta > \delta$ and, by the first form of (3.3a),

and by the definitions of $c$ and $\alpha$, $\alpha > c(1-{w}^{*})$.

(iii) Now, suppose that $\frac{3\pi }{2}\ge yT$, so that $\; {\mathrm{cos}}\; yT\le 0$. Since ${\mathrm{Re}}\zeta =\alpha (1-{e}^{-xT}\; {\mathrm{cos}}\; yT)$, it follows that $\alpha \le {\mathrm{Re}}\zeta$. Therefore, by the first form of (3.3a),

which contradicts part ii. Consequently, $\; {\mathrm{cos}}\; yT > 0$, $\frac{3\pi }{2} < yT$ and, by part i, $x < \frac{\beta }{3\pi -1}$.

(iv) Returning to (5.2a) once more, with the added information that $\; {\mathrm{cos}}\; yT > 0$, we find

${y}^{2}-\alpha y-{\alpha }^{2}g\left(x\right) < 0$ where $g\left(x\right)=({\frac{x}{\alpha })}^{2}+\left(1+\frac{\delta }{\alpha }\right)\frac{x}{\alpha }+\frac{\beta }{\alpha }+\frac{m}{\alpha }\bullet \frac{a}{\alpha }\bullet q{v}^{*}$

By the second inequality of part ii and the fact that $c\ge a$,

By the second inequality of part iii, the first inequality of part ii, and (3.3b),

Finally, by the first inequality of part iii and what we have just done,

In ^[1], there are results which pertain to characteristic functions of the form $P\left(\lambda \right)-Q\left(\lambda \right){e}^{-\lambda T}$ where $P$ and $Q$ are real polynomials, $P$ has degree 2, its coefficient of $\lambda$ 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)\left(0\right)\ne 0$. If $f\left(y\right)={\left|P\left(iy\right)\right|}^{2}-{\left|Q\left(iy\right)\right|}^{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\left(z\right)-Q\left(z\right){e}^{-Tz}$ have negative real parts as well. Conversely, if $P-Q$ has a root with nonnegative real part, then $P\left(z\right)-Q\left(z\right){e}^{-Tz}$ has zeros with nonnegative real part for any $T > 0$.

We apply this theorem to our situation. Let $P\left(z\right)={z}^{2}+\delta z+maq{v}^{*}+\alpha \left(z+\beta \right)$ and $Q\left(z\right)=\alpha (z+\beta)$, so that $G\left(\lambda \right)=P\left(\lambda \right)-Q\left(\lambda \right){e}^{-\lambda T}$. Then the root of $Q$ is real and $\left(P-Q\right)\left(0\right)=maq{v}^{*}\ne 0$. In addition, since $\delta > 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.

The positive real roots of this quartic are the square roots of the positive roots of

We need to find conditions under which $F$ has no positive roots, i.e., either that

Therefore, by the last inequality above, it is sufficient that

This will be true if and only if $\delta$ exceeds the positive root of the quadratic on the right-hand side. For the moment, let $A=\frac{\beta }{1+\frac{maq{v}^{*}}{\alpha \beta }}$. Then, the positive root of the quadratic is $\alpha \left(\sqrt{1+2\frac{A}{\alpha }}-1\right)$. Rationalizing the numerator, we have

Consequently, it suffices that $\delta \ge 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.

6.   Discussion

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 ${\nu }_{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 ${\nu }_{0}$, is investigated with some success. The extinction of virus and infected cells, represented by $\left(0, 0\right)$ 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\left(b{e}^{-mT}-1\right)+lq > h+k$, is the presence of a second, positive, fixed point $\left({v}^{*}, {w}^{*}\right)$ 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\left(b{e}^{-mT}-1\right)\ge 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 $\left({v}^{*}, {w}^{*}\right)$ only in special circumstances. However, we can show that $\left(0, 0\right)$ 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 $\left({v}^{*}, {w}^{*}\right)$ are negative. Interestingly, one such condition is $m\ge 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$?

Acknowledgments

Thanks are given to the academic editor and reviewers for helpful suggestions.

Conflict of interest

The author declares that there is no conflict of interest.