# Understanding a dynamical system for virus populations

I am reading the book: M. A. Novak, R. M. May: Virus Dynamics: Mathematical Principles of Immunology and Virology. Oxford University Press, 2000.

I came across few places that I don’t understand. In the model, uninfected cells react with free virus to give rise to infected cells; the rate constant is $\beta$. Infected cells produce free virions at rate $k$. Uninfected cells, free virus and infected cells die at rates $d$, $u$ and $a$, respectively. Uninfected cells are replenished at rate $\lambda$.

The schematic of the model is:

The model equations are [eq. (3.1) p. 18 in the link]:

\begin{alignat}{1} \dot{x} &= λ - dx - βxv,\\ \dot{y} &= βxv - ay, \\ \dot{v} &= ky - uv, \end{alignat}

where $x$, $y$, and $v$ are the populations of uninfected cells, infected cells and virions.

The basic reproductive ratio, $R_0$ is the number of newly infected cells that arise from any infected cell when almost all cells are uninfected (and the system is near its equilibrium):

$$R_0= \frac{\beta \lambda k}{adu}$$

What I don’t understand is:

1. It is mentioned that if $R_0<1$ the virus will not spread since every infected cell will produce on average less than one infected cell. If we start with $N$ infected cells, then on average, we expect roughly ${\ln N}\over {\ln(1/R_0)}$ rounds of replication before the virus population dies out. How is this ${\ln N}\over {\ln(1/R_0)}$ found?

2. If $R_0>1$ then virus will initially grow exponentially, and $r_0$ is the exponential growth rate of the population. Then it is said $r_0$ is given by the larger root of the equation $r_0^2+(a+u)r_0+au (1-R_0)=0$. How did they come up with this equation?

3. Then they say in the equation $r_0^2+(a+u)r_0+au (1-R_0)=0$ if $u \gg a+r_0$ we find the approximation $r_0=a(R_0-1)$. How is this $r_0=a(R_0-1)$ obtained?

• Have you at least found a solution to the three differential equations? – marshal craft Aug 31 '17 at 15:14
• Basically I think you should be able to solve for $x,y,v$ then eliminate them yielding equation involving only the above rates etc. Or put them in form equal to $R_0$. – marshal craft Aug 31 '17 at 15:20
• Wrzlprmft, is that a medical peculiarity, can one not still consider the question as asked hypothetically? Or somehow is it mathematically over constrained? – marshal craft Aug 31 '17 at 15:37
• To me it seems necessary, you have flow of uninfected cell in, flow out, virus flow in per infected cell, virus flow out, infected cell flow in and flow out. Three things uninfected cells, virus, infected cells. – marshal craft Aug 31 '17 at 15:40
• @marshalcraft "you should be able to solve for x,y,v" Typically this is the one task one cannot achieve, hence the need to turn to single parameters $R_0$ and $r_0$ describing (qualitatively) (some aspects of) the solutions. – Did Aug 31 '17 at 18:47

The equations are related to discrete approximations of the infection process of the virus.

1. Let $F_n$ be the number of infected cells at replication step $n$. Then, we have that $F_{n+1} \approx R_0 F_n$, $F_0=N$. Then we can solve this recurrence relation to get $F_n = R_0^nN$. We want to know for what $n$ we get $F_n <1$. So we solve $R_0^n N <1$ for $n$. Taking logarithms we get that $$n\ln(R_0)+\ln(N) < 1$$
subtracting, and dividing by $\ln(R_0)$ (and remember that $ln(R_0) < 0$ since $R_0 <1$, so we have to switch the direction of the inequality) we get $$n > \frac{-\ln(N)}{\ln(R_0)} = \frac{\ln(N)}{\ln \frac{1}{R_0}}$$

2. Here is my best guess for this one. Suppose that we insert a single virion into the population when it is at equilibrium. We then have that as time progresses $v(t) = e^{r_0t}$. Plugging this into the $\dot{v}$ equation gives that $ky=(r_0+u)e^{r_0t}$. Differentiating the $\dot{v}$ equation again gives $$\ddot{v} + u\dot{v} = k\beta x v - a(ky).$$ Using the above expression for $ky$ and the fact that $x \approx \frac{\lambda}{d}$ we get $$\ddot{v} + u\dot{v} - \frac{k \beta \lambda}{d}v = -a(r_0+u)e^{r_0t}.$$ Plug in $v = e^{r_0t}$, which then allows you to cancel the $e^{r_0t}$ from both sides and this gives $$r_0^2+ur_0- \frac{k \beta \lambda}{d} = -ar_0-au.$$ Thus $$r_0^2 + (u+a)r_0 + au-\frac{k\beta\lambda}{d}=0$$ Then, finally, we factor $au$ out of the constant terms and use the definition of $R_0$ to get $$r_0^2 + (u+a)r_0 + au (1-R_0) = 0.$$

3. Expand the equation to get $r_0^2 + a r_0 + ur_0 + au(1-R_0) = 0$. Dividing both sides by $u$ gives $$r_0 \left( \frac{r_0+a}{u} \right) + r_0+a(1-R_0)=0$$ Since $u \gg r_0+a$ we get that $\frac{r_0+a}{u} \ll 1$. since the leading term in the above equation is so small, it can be neglected. Tossing it out gives $$r_0 + a(1-R_0) = 0 \Rightarrow r_0 = a(R_0-1).$$

• Thank you very much for the answer. Can you tell me how you decided to differentiate $\dot{v}$ again in the first place. Also, have they decided to chose $u \gg r_0+a$ just to make the math simple or is it based on some biological meaning? – clarkson Sep 1 '17 at 1:13
• So, If you look at Wrzlprmft's answer below, it will give a less ad hoc type of way at arriving at the answer. However, the reason that I chose to differentiate it again is because it seemed like the quadratic equation was the characteristic equation for a second order ODE. Since it related to $v$, it made sense to differentiate the $v$ equation. – kareemmatheson Sep 1 '17 at 17:41
• As for the size condition, it says that the rate that the free virus "dies" is much higher than both the rate at which it is reproducing (inside the cells) and the rate at which infected cells are dying. This suggests to me that the condition is saying that the environment outside the cell is extremely hostile to the virus. I haven't done much in-host viral modeling, but that would be my best guess. – kareemmatheson Sep 1 '17 at 17:45
• @ kareemmatheson @Wrzlprmft Would it be possible to find $r_o$ like term in a much complexed system with around 6 non linear ODEs like the system that is found in this article in equations 6-11. ncbi.nlm.nih.gov/pmc/articles/PMC3066295 – clarkson Sep 4 '17 at 11:33
• @clarkson: You can always linearise the system and apply what I described in my answer. However, you should be aware that you make some approximations here (linearity, virus populations are not discrete, largest eigenvector dominates dynamics) that have a limited range of validity. You have to check whether these apply. – Wrzlprmft Sep 4 '17 at 11:43

An alternative way to solve problem 2 using some knowledge on dynamical systems:

We first insert our knowledge of being near the uninfected equilibrium (with $x=\frac{λ}{d}$) into the differential equations:

\begin{alignat}{4} \dot{y} &=~& - ay &~+~& \frac{λβ}{d} v, \\ \dot{v} &=& ky &~-~& uv, \end{alignat}

This is a linear system of differential equations, which can be rewritten as a matrix–vector multiplication (whose matrix I denote as $A$):

$$\frac{\mathrm{d}}{\mathrm{d}t} \begin{pmatrix}y\\v\end{pmatrix} = \begin{pmatrix} -a & \frac{λβ}{d} \\ k & -u \end{pmatrix} \begin{pmatrix}y\\v\end{pmatrix} =: A \begin{pmatrix}y\\v\end{pmatrix}$$

The solutions to such a differential equation are known to be of the form:

$$α_1 e^{ρ_1 t} w_1 + α_2 e^{ρ_2 t} w_2,$$

where the $w_i$ are the eigenvectors of $A$, the $ρ_i$ are the corresponding eigenvalues and the $α_i$ are constants determined by the initial conditions. Now, the component corresponding to the largest eigenvalue will dominate the other one and thus we get exponential growth with this eigenvalue as a growth rate. (As $v$ and $y$ depend on each other, everything grows with the same rate.)

So, all that is left to do is to determine the characteristic polynomial for $A$ (whose roots are the eigenvalues):

$$p_A(r_0) = \det(r_0𝟙-A) = \det\pmatrix{r_0+a & \frac{-λβ}{d} \\ -k & r_0+u} = (r_0+a)(r_0+u) - \frac{kλβ}{d} \\ = r_0^2 +(a+u)r_0 + au - auR_0 = r_0^2 +(a+u)r_0 + au (1 - R_0)$$