Viral Infection Dynamics Model Based on a Markov Process with Time Delay between Cell Infection and Progeny Production

Abstract

Many human virus infections including those with the human immunodeficiency virus type 1 (HIV) are initiated by low numbers of founder viruses. Therefore, random effects have a strong influence on the initial infection dynamics, e.g., extinction versus spread. In this study, we considered the simplest (so-called, ‘consensus’) virus dynamics model and incorporated a delay between infection of a cell and virus progeny release from the infected cell. We then developed an equivalent stochastic virus dynamics model that accounts for this delay in the description of the random interactions between the model components. The new model is used to study the statistical characteristics of virus and target cell populations. It predicts the probability of infection spread as a function of the number of transmitted viruses. A hybrid algorithm is suggested to compute efficiently the system dynamics in state space domain characterized by the mix of small and large species densities.

1. Introduction

Since the mid 1990s, mathematical modelling of human virus infections has been intensively developed for various infections including those with the human immunodeficiency virus type 1 (HIV) and the hepatitis B virus. Reviews on existing approaches are presented, i.e., in [1,2,3,4]. Mathematical modelling of virus-target cell interactions based on the understanding of the underlying biological processes can provide deeper insights into the mechanisms of the infection dynamics, see [3,5,6,7]. The majority of these studies are based on deterministic descriptions with continuous populations of virus particles (virions) and target cells. This approximation is accurate when the number of virions, infected cells and other interacting components are large enough that is typical for the later stages of an infection process. However, this deterministic framework has its drawbacks because an infection process is intrinsically stochastic. Populations of interacting species undergo stochastic fluctuations that are more profound at the early stage of an infection when the number of virions and infected cells are still small. Furthermore, random fluctuations may lead to extinction of an infection process. Such effects are well suited for being analyzed with methods of stochastic processes theory.

The stochasticity of a virus infection may be accounted for in models using stochastic differential equations (SDEs) governed by Brownian motions (BMs), i.e., [8,9,10,11]. This modelling approach seems more realistic and can explain phenomena like virus extinction (cf. [12]). However, the applicability of such an approach is restricted by the daunting task of parameter identification for processes like volatility or the probability density function (PDF). For example, the volatility of the Brownian process varies with the population size of interacting components and that is not easy to account for.

Another approach to model stochastic infection dynamics is the use of a discrete or continuous Markov Chain (MC) in the framework of the Monte-Carlo method. This approach and algorithms for numerical simulations have been developed first for chemical kinetics (cf. [13,14]). Pearson et al. [15] considered a two-component model for early infections with the use of a MC. A similar simulation of the ’consensus’ three-component virus dynamics model (cf. [2]) was proposed and studied in [16]. In such an approach, there is no necessity to work out parameters of Brownian motion or other processes because the transition rates (propensities) for the Markov chain are determined directly from parameters of the deterministic equations. A comparison between the SDE approach and the discrete and continuous MC approach for a simple population dynamics model can be found in [17].

Insightful models of virus infections may take the intracellular phase of the virus life-cycle into account and thus, the time delay between infection of a cell and the production of new virus progeny (cf. [18]). A deterministic model with a fixed time delay has been proposed by Herz et al. [19] (see also [20,21]). A more sophisticated model accounting for two distributed time delays was developed in [22]. Stochastic modelling of infection dynamics with fixed and distributed time delays based on the SDE approach can be found in [11,23,24].

In our present work we generalized the ideas of [16] on a model with time delay and used a Markov Process (MP) with a fixed delay instead of a Markov Chain (MC) considered in [16]. This complicates computation but makes the model more realistic. We also propose a hybrid model for effective computation of statistical parameters of the stochastic virus dynamics.

In Section 2, we describe the well-known deterministic model of virus infection dynamics with delay. The related stochastic model is presented in Section 3. In Section 4, the hybrid modelling algorithm is developed. We summarize the results in Section 5.

2. Deterministic Model

The standard and classic ODE system for the three species virus dynamics has been introduced in [25,26] (see also [1,2,19]):

$$\begin{array}{cc}\hfill \dot{x}& =\lambda -dx-\beta xv\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill \dot{y}& =\beta xv-ay\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill \dot{v}& =ky-uv\hfill \end{array}$$

(1c)

where x, y and v are, respectively, the numbers of uninfected cells, infected cells and the number of free virus particles (virions) in a fixed volume compartment. The authors suppose that uninfected cells are produced at a constant rate, $\lambda$, and die at a rate $dx$; virions infect uninfected cells at a rate proportional to the product of their numbers, $\beta xv$; infected cells produce free virus at a rate $ky$; infected cells die at a rate $ay$; free virus particles are removed from the system at a rate $uv$.

Herz et al. [19] have modified ODEs (1) by including the delay accounting for a latent period between the time when target cells are contacted with virions and the time when proviral DNA integrates into the cellular genome and proposed the following DDE system

$$\begin{array}{cc}\hfill \dot{x}& =\lambda -dx-\beta xv\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill \dot{y}& =\beta x(t-\tau )v(t-\tau )e^{-a\tau }-ay\hfill \end{array}$$

(2b)

$$\begin{array}{cc}\hfill \dot{v}& =ky-uv\hfill \end{array}$$

(2c)

where $\tau$ is the time lag between the fusion of the virion with the cell membrane and integration of proviral HIV DNA in the cell genome. The term $\beta x(t-\tau )v(t-\tau )e^{-a\tau }$ indicates that the secretion rate for infected cells is proportional to the number of uninfected cells and virions at the time $t-\tau$ decreased by factor $e^{-a\tau }$ because of natural and immune-mediated death infected cells with the rate a in time $\tau$. Variables in all other terms are considered at time t.

Ciupe et al. [27] have proposed a four species model with account for the second type delay: the time lag between virion’s DNA penetrates the cell and new virions are produced and released. We present here a similar, three species DDE version:

$$\begin{array}{cc}\hfill \dot{x}& =\lambda -dx-\beta xv\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill \dot{y}& =\beta xv-ay\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill \dot{v}& =ky(t-\tau )e^{-a\tau }-uv.\hfill \end{array}$$

(3c)

Here, $\tau$ is the time delay between penetration of a virion into a cell and release of new virions. Term $ky(t-\tau )e^{-a\tau }$ indicates that the growth rate for virions is proportional to the number of cells infected by virus at the time $t-\tau$ decreased by a factor $e^{-a\tau }$ because of natural and immune-mediated death of infected cells with the rate a in time $\tau$.

A more general model accounting of four species interaction for the both types of delay has been proposed by Pawelek et al. [22]. Its three species version is studied in [28].

Applying results obtained in [28] to DDEs (3), we can show that the basic reproduction number of the model is

$$R_0=\frac{\beta \lambda k}{adu}{e}^{-a\tau }.$$

(4)

The model has two equilibrium states: the infection-free equilibrium $\{x=x_0,y=0,v=0\}$ where $x_0=\lambda /d$ and the infection equilibrium $\{x=x_{\infty },y=y_{\infty },v=v_{\infty }\}$ where

$$x_{\infty }=\frac{x_0}{R_0},y_{\infty }=\frac{du}{\beta k}\left(R_0-1\right)e^{a\tau },v_{\infty }=\frac{d}{\beta }\left(R_0-1\right).$$

(5)

The infection equilibrium exists and is asymptotically stable if $R_0>1$.

Let $v_0$ virions arrive at instant $t=0$. It is natural to take $x_0=\lambda /d$ as an initial number of uninfected cells in the virus dynamics and to assume that there were no infected cells before first virions arrived. Then we have the following initial conditions

$$x\left(0\right)=x_0,y\left(0\right)=0,v\left(0\right)=v_0.$$

(6)

For DDEs it is necessary to set also the history function defined in the interval $[-\tau ,0)$:

$$x\left(t\right)=x_0,y\left(t\right)=0,v_0\left(t\right)=0,t\in [-\tau ,0).$$

(7)

Equation (3) with initial conditions (6) and history (7) form a complete initial value problem (sometimes called the initial data problem) for the DDEs and uniquely define the virus dynamics for $t>0$. They can be numerically integrated, for example, with the use of Matlab function dde23().

Nowak and May [2] have elaborated typical parameters for a typical HIV virus infection process:

$$\lambda ={10}^5,d=0.1,a=0.5,\beta =2\times {10}^{-7},k=100,u=5.$$

(8)

These parameters have units of inverse days: ${\mathrm{d}}^{-1}$. These values have been used in a number of works including our previous work [16].

Examples of computations of the DDEs with parameters (8) and different types of delay are shown in Figure 1. The same value $\tau =1$ d is used in both cases: for the first (2) and the second (3) types of delay. We denote the solution for the first case as $x_1\left(t\right),y_1\left(t\right),v_1\left(t\right)$ and for the second case as $x_2\left(t\right),y_2\left(t\right),v_2\left(t\right)$.

Observe that for the same delay, the dynamics in both cases is identical for the number of uninfected cells $x_1\left(t\right)\equiv x_2\left(t\right)$ and the number of virions $v_1\left(t\right)\equiv v_2\left(t\right)$ for all t, whereas the dynamics of infected cells differs by time shifting and scaling: $y_1\left(t\right)=y_2(t+\tau )e^{-a\tau }$. This means that contribution of both delays is similar into the virus dynamics.

This fact motivates a detailed study of the model with the second type delay (3) to build an equivalent stochastic model. Building of a stochastic model with the first type delay (2) is more involved and will be considered in the subsequent works.

3. Stochastic Model

3.1. Markov Process Approach

Now we account for a discrete nature of the populations and random interactions between different species. We consider an approach based on Markov process formulation extending directly DDEs (3).

Let $X,Y,V\in {\mathbb{Z}}^{+}$ be non-negative integer numbers of uninfected, infected cells and virions, respectively. Then we assume that their dynamics obeys the following time continuous Markov process (MP) with delay (

Table 1. Markov Process with Delay.

#	Transition	Propensity	Description
1	$X\left(t\right)\to X\left(t\right)+1$	$\lambda$	$\mathrm{Uninfected}\mathrm{cell}\mathrm{birth}$
2	$X\left(t\right)\to X\left(t\right)-1$	$dX\left(t\right)$	$\mathrm{Uninfected}\mathrm{cell}\mathrm{death}$
3	$\begin{array}{c}X\left(t\right)\to X\left(t\right)-1\hfill \\ Y\left(t\right)\to Y\left(t\right)+1\hfill \end{array}$	$\beta X\left(t\right)V\left(t\right)$	$\mathrm{Cell}\mathrm{infection}$
4	$Y\left(t\right)\to Y\left(t\right)-1$	$aY\left(t\right)$	$\mathrm{Infected}\mathrm{cell}\mathrm{death}$
5	$V\left(t\right)\to V\left(t\right)+1$	$k_{\tau }Y(t-\tau )$	$\mathrm{Virion}\mathrm{birth}$
6	$V\left(t\right)\to V\left(t\right)-1$	u$V\left(t\right)$	$\mathrm{Virion}\mathrm{death}$

).

Where

$$k_{\tau }=k{e}^{-a\tau }$$

(9)

is an effective (cell-death adjusted) free virus production rate (per cell).

The initial conditions for MP (Table 1) are

$$X\left(0\right)=X_0,Y\left(0\right)=0,V\left(0\right)=V_0,X_0,V_0\in {\mathbb{Z}}^{+}$$

(10)

where $X_0$ has the Poisson distribution with the rate $x_0=\lambda /d$

$$P(X_0=n)=\frac{x_0^n}{n!}{e}^{-x_0}.$$

(11)

in line with [16]. For large $x_0$ it can be approximated by the Gaussian distribution with mean and variance equal to $x_0$. Therefore the coefficient of variation (CV), defined as the standard deviation over the mean, equals $x_0^{-1/2}$ that gives ${10}^{-3}$ for $x_0={10}^6$. Thus the probability density function (PDF) is very narrow for realistic numbers of uninfected cells x, so the effect of the stochasticity is negligible and $X_0$ can be approximated by the integer closest to $x_0=\lambda /d$.

A Markov process with delay also should include history at time $t\in [-\tau ,0)$. Obviously for the virus dynamics we have absence of infected cells before the infection started:

$$Y\left(t\right)=0,t\in [-\tau ,0).$$

(12)

Relation (Table 1) implies that variables $X,Y,V$ are non-negative as the propensity of transitions 2, 3, 4, 6, in which the number of one of the species decreases, vanishes as soon as $X=0$ (transitions 2 and 3) or $Y=0$ (transition 4) or $V=0$ (transition 6).

Process (Table 1) can be re-written in terms of Poisson processes [29]:

$$X\left(t\right)=X\left(0\right)+{\mathrm{Po}}_1\left(\lambda t\right)-{\mathrm{Po}}_2\left(d{\int }_0^{t}X\left(s\right)\mathrm{d}s\right)-{\mathrm{Po}}_3\left(\beta {\int }_0^{t}X\left(s\right)V\left(s\right)\mathrm{d}s\right)$$

(13)

$$Y\left(t\right)=Y\left(0\right)+{\mathrm{Po}}_3\left(\beta {\int }_0^{t}X\left(s\right)V\left(s\right)\mathrm{d}s\right)-{\mathrm{Po}}_4\left(a{\int }_0^{t}Y\left(s\right)\mathrm{d}s\right)$$

(14)

$$V\left(t\right)=V\left(0\right)+{\mathrm{Po}}_5\left(k{e}^{-a\tau }{\int }_0^{(t-\tau )\vee 0}Y\left(s\right)\mathrm{d}s\right)-{\mathrm{Po}}_6\left(u{\int }_0^{t}V\left(s\right)\mathrm{d}s\right)$$

(15)

where ${\mathrm{Po}}_i(·),i=1,\dots ,6$ are independent Poisson processes.

Consider the scaled MP $(\widehat{X}\left(t\right),\widehat{Y}\left(t\right),\widehat{V}\left(t\right))$

$$\widehat{X}\left(t\right)={\mathsf{\Lambda }}^{-1}X\left(\mathsf{\Lambda }t\right),\widehat{Y}\left(t\right)={\mathsf{\Lambda }}^{-1}Y\left(\mathsf{\Lambda }t\right),\widehat{V}\left(t\right)={\mathsf{\Lambda }}^{-1}V\left(\mathsf{\Lambda }t\right)$$

(16)

described by the following equations

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \widehat{X}\left(t\right)& =\widehat{X}\left(0\right)+{\mathsf{\Lambda }}^{-1}{\mathrm{Po}}_1\left(\mathsf{\Lambda }\lambda t\right)-{\mathsf{\Lambda }}^{-1}{\mathrm{Po}}_2\left(\mathsf{\Lambda }d{\int }_0^t\widehat{X}\left(s\right)\mathrm{d}s\right)\hfill \\ \hfill & -{\mathsf{\Lambda }}^{-1}{\mathrm{Po}}_3\left(\mathsf{\Lambda }\beta {\int }_0^t\widehat{X}\left(s\right)\widehat{V}\left(s\right)\mathrm{d}s\right)\hfill \end{array}\end{array}$$

(17a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \widehat{Y}\left(t\right)& =\widehat{Y}\left(0\right)+{\mathsf{\Lambda }}^{-1}{\mathrm{Po}}_3\left(\mathsf{\Lambda }\beta {\int }_0^t\widehat{X}\left(s\right)\widehat{V}\left(s\right)\mathrm{d}s\right)\hfill \\ \hfill & -{\mathsf{\Lambda }}^{-1}{\mathrm{Po}}_4\left(\mathsf{\Lambda }a{\int }_0^t\widehat{Y}\left(s\right)\mathrm{d}s\right)\hfill \end{array}\end{array}$$

(17b)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \widehat{V}\left(t\right)& =\widehat{V}\left(0\right)+{\mathsf{\Lambda }}^{-1}{\mathrm{Po}}_5\left(\mathsf{\Lambda }k{e}^{-a\tau }{\int }_0^{(t-\tau )\vee 0}\widehat{Y}\left(s\right)\mathrm{d}s\right)\hfill \\ \hfill & -{\mathsf{\Lambda }}^{-1}{\mathrm{Po}}_6\left(\mathsf{\Lambda }u{\int }_0^t\widehat{V}\left(s\right)\mathrm{d}s\right)\hfill \end{array}\end{array}$$

(17c)

were $\widehat{X}\left(0\right)=X\left(0\right),\widehat{Y}\left(0\right)=Y\left(0\right),\widehat{V}\left(0\right)=V\left(0\right)$ and parameters $\lambda$, d, $\beta$, a, k, u and $\tau$ are properly scaled.

Proposition 1.

If $(\widehat{X}\left(0\right),\widehat{Y}\left(0\right),\widehat{V}\left(0\right))\to (x_0,y_0,v_0)$ then for any $T>0$ the scaled process (17) weakly converges as $\mathsf{\Lambda }\to \infty$ to the solution of the system (3) with the initial condition $(x_0,y_0,v_0)$

$$(\widehat{X}\left(t\right),\widehat{Y}\left(t\right),\widehat{V}\left(t\right))\to \left(x\left(t\right),y\left(t\right),v\left(t\right)\right)$$

(18)

on the interval $0\le t\le T$.

The proof is presented in Appendix A.

Algorithm for numerical simulation of this model is proposed by Anderson [30] (namely, Algorithm 7 in this work) and based on the next reaction method due to Gibson and Bruck [31]. It is also described and studied by Banks et al. [29] (namely, Algorithm 2 in this work).

3.2. Direct Numerical Simulation

We will study the virus dynamics described by MP (Table 1) numerically. The simulation algorithm is a slightly modified Algorithm 2 described in [29].

To describe the numerical algorithm modelling MP (Table 1) we introduce the state vector $\mathit{X}={[X\left(t\right),Y\left(t\right),V\left(t\right)]}^T$ and the $N\times M$ matrix of transitions

$$\mathit{T}=\left[\begin{array}{cccccc}\hfill +1& \hfill -1& \hfill -1& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill +1& \hfill -1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill +1& \hfill -1\end{array}\right]$$

(19)

where $N=3$ is the number of variables, $M=6$ is the number of transitions in the process. If the mth transition occurs at a certain instant, then the mth column of matrix $\mathit{T}$ should be added to the state vector $\mathit{X}$.

For the MP with delay we have to arrange two arrays: ${\mathit{t}}^{\mathrm{Y}}$ to store the instants when Y has been changed and ${\mathit{Y}}^{\mathrm{Y}}$ to store the number of infected cells Y at those instances. These arrays contain a history for the Y component of the MP.

To handle arrays ${\mathit{t}}^{\mathrm{Y}}$ and ${\mathit{Y}}^{\mathrm{Y}}$, we use two counters $i_1$ and $i_2$ where $i_1$ indicates the position of the the earliest event to be accounted in ${\mathit{t}}^{\mathrm{Y}},{\mathit{Y}}^{\mathrm{Y}}$ while $i_2$ indicates the position of the last event stored in those arrays. The set of events to be accounted is not empty if $i_2\ge i_1$.

The process dies out at instant t if $V\left(t\right)=0$ and $Y\left(t\right)=0$ and the set of events to be accounted is empty, i.e., $i_2<i_1$. This is a difference from the process without delay for which the condition of extinction is simply $Y\left(t\right)=V\left(t\right)=0$.

To model the process, we introduce a propensity vector $\mathit{\nu }=[{\nu }_1,\dots ,{\nu }_M]$, a vector of integrals $\mathit{J}=[J_1,\dots ,J_M]$, and also vector $\mathit{R}=[R_1,\dots ,R_M]$ where $R_m=-ln\left(r_m\right)$ with $r_m$ be random numbers uniformly distributed on $[0,1]$.

In view of Table 1 the calculation of the propensity vector $\mathit{\nu }$ depends on four components: $X\left(t\right),Y\left(t\right),V\left(t\right),Y(t-\tau )$, therefore it is convenient to introduce the augmented state vector $\check{\mathit{X}}=[X\left(t\right),Y\left(t\right),V\left(t\right),Y(t-\tau )]$.

We also introduce an augmented vector of expected timesteps for all the transitions $\mathit{t}=[t_1,\dots ,t_{M+1}]$ which contains an additional component $t_{M+1}$. Time $t_{M+1}-\tau$ indicates an instant when the number of infected cells Y has been changed, thus $t_{M+1}$ is the instant when this change affects the dynamics.

3.3. Full Stochastic Numerical Simulations Results

Algorithm 1 for the direct numerical simulations of MP (Table 1) has been implemented in the C language with the use of the PCG library [32] for the random number generator. The script calling the program in loop to obtain $N_r={10}^3$ non-extinct realizations (see below) for various $V_0$ is written in the GNU Bash. It utilizes the GNU parallel tool [33] to parallelize executions of different realizations on a maximum available number of threads.

Algorithm 1 Implementation of the Markov process with delay.

Set the current time $t=0$ and the final time of the process $t_f$; Initialise the augmented state vector $\check{\mathit{X}}={[X_0,0,V_0,0]}^T$ and set $\mathit{J}$: $J_m=0,m=1,\dots ,M$; Compute the propensity vector $\mathit{\nu }\left(\check{\mathit{X}}\right)$ using values of propensity indicated in Table 1; Generate M random numbers $r_m$ and compute vector $\mathit{R}$: ${\displaystyle R_m=-ln\left(r_m\right),m=1,\dots ,M}$; Compute the time steps to next event in all the transitions: ${\displaystyle \Delta t_m=R_m/{\nu }_m,m=1,\dots ,M}$ and set $\Delta t_{M+1}=\infty$; Allocate arrays ${\mathit{t}}^{\mathrm{Y}}$ and ${\mathit{y}}^{\mathrm{Y}}$, set $t_1^{\mathrm{Y}}=0$, $y_1^{\mathrm{Y}}=0$; set counters $i_1=1$, $i_2=0$ for these arrays; Find transition p with the minimal time step: $\Delta t_p=min\{\Delta t_1,\dots ,\Delta t_{M+1}\}$; Update vector $\mathit{J}$: $J_m=J_m+{\nu }_m\Delta t_p$, $m=1,\dots ,M$; Update the current time $t=t+\Delta t_p$; If$p<M+1$then (a) update the state vector $\check{\mathit{X}}$: ${\check{X}}_n={\check{X}}_n+T_{np},n=1,2,3$ (b) set $J_p=0$ (c) if$p=3$ or $p=4$ (transition in which Y has been changed) then set $i_2=i_2+1$, $Y^{\mathrm{Y}}\left(i_2\right)=Y$, $t_{i_2}^{\mathrm{Y}}=t^{\prime }$ If$p=M+1$then (a) update the state vector $\check{\mathit{X}}$: ${\check{X}}_{M+1}=Y_{i_1}$ (b) set $i_1=i_1+1$ (c) if$i_2<i_1$ and ${\check{X}}_2=\check{X_3}=0$ then flag the realisation as extinct; Update the propensity vector $\mathit{\nu }=\mathit{\nu }(\check{\mathit{X}})$; If$p\le M$then generate random number r and compute $R_p=-ln\left(r\right)$; Update the time step vector $\Delta t_m=(R_m-J_m)/{\nu }_m$, $m=1,\dots ,M$; If$i_2\ge i_1$then set $\Delta t_{M+1}=t^Y\left(i_1\right)-t$ else set $\Delta t_{M+1}=\infty$; Store the current state and time; If$t<t_f$then go to 7 otherwise terminate the computation.

Two examples of the MP numerical simulation are presented in Figure 2. The left plot represents an established viral infection process. A realization in which numbers of virions and infected cells can reach their peak values is a non-degenerate realization and corresponds to a developed infection disease. We observe relatively high stochastic fluctuations in the number of virions and infected cells at the very earlier stage of infection when those numbers are not large. Those fluctuations become relatively low at later stages with increase of the numbers of all species. To make the fluctuations more visible parameter $\lambda$ is decreased and parameter $\beta$ is increased by a factor of 10. Such scaling does not change the deterministic dynamics for the same ratio $v_0/x_0$.

The plot in Figure 2b represents the degenerated (extinct) realization. In contrast to the deterministic case, if the infection is modelled via Markov process the infection dynamics can extinct in some realizations, i.e., it can reach the state $V\left(t_e\right)=0$, $Y(t\in [t_e-\tau ,t_e])\equiv 0$ at a certain time $t_e$ after which the viral dynamics is terminated. As usual this occurs at the eclipse phase when numbers of virions and uninfected cells are still not too large.

The probability of infection establishment, $P_i$, and the extinction probability, $P_e=1-P_i$, strongly depend on the initial number of virions $V_0$ but also on the delay time $\tau$. The described above numerical simulation enables computation of these probabilities and their dependance on parameters of the model by counting the number of non-degenerated (infection) realization $N_i$ and degenerated (extinct) realizations:

$$P_i=\underset{N_i,N_e\to \infty }{lim}\left(N_i/N\right),P_e=\underset{N_i,N_e\to \infty }{lim}\left(N_e/N\right)$$

(20)

where $N=N_i+N_e$ is the number of all realizations.

For finite N, the $N_i/N$ and $N_e/N$ ratios are distributed as the sum ${\sum }_{j=1}^N{\xi }_j$ where ${\xi }_j$ are independent and identically distributed (iid) Bernoulli random variables with the success probability $P_i$ and $P_e$, respectively. We can show that, if $P_i$ and $P_e$ are not close to zero or unity then the probability density function (PDF) of $N_i/N$ and $N_e/N$ for $N\to \infty$ tend to the Gaussian law with the standard deviation

$${\sigma }_{P_i}={\sigma }_{P_e}=\sqrt{N{P}_i{P}_e}\approx \sqrt{N_i{N}_e/N}$$

(21)

Equation (21) enables estimation of the accuracy of computation of $P_{i,e}$ for given number of realizations. The accuracy within two-three significant digits can be guaranteed for $N={10}^6$ taken in our modelling study.

The computed dependence of $P_i$ on $\tau$ for different $V_0$ and fixed other parameters (8) is shown in Figure 3a.

As it is seen in the figure, the infection probability $P_i$ decays fast when the time delay $\tau$ approaches four days. Remind that in virtue of Equation (4), the basic reproduction number exponentially depends on $\tau$ decreasing from $R_0=8$ at $\tau =0$ down to unity at $\tau =4.16$ d. For this reason one can expect that the infection probability $P_i$ vanishes at $\tau \ge 4.16$ d. However, the plots in Figure 3a show that $P_i$ remains finite even $\tau >4.16$ d.

To explain this phenomenon we recall that a formally non-degenerated process does not represent a developed infection dynamics when $R_0$ is close to unity. The infection peak is developing very slow: hundreds of days (that is not typical for HIV infection), the number of virions is smaller than the number of uninfected cells, the number of uninfected cells varies insignificantly. Therefore, when the time delay $\tau$ exceeds four days, the process with the fixed parameters (8) does not reproduce the observed HIV viral infection dynamics, and the parameters (8) should be adjusted.

Note that if we keep the parameter k to be fixed then the effective free virions production rate $k_{\tau }=k{e}^{-a\tau }$ decays exponentially with the growth of time delay $\tau$. Remember that the free virions production rate $k=100$ d${}^{-1}$ has been elaborated by Nowak and May [2] from the available HIV viral dynamics data (although without delay). Therefore the essential reduction in the actual virus producing rate with the growth of time delay looks unnatural. It seems more reasonable to consider another scenario to study the effect of the time delay variation on the system dynamics: to keep the actual free virus producing rate $k_{\tau }$ to be constant while varying the time delay $\tau$. In the model studied here we can set $k_{\tau }=100$ d${}^{-1}$ as it is found in [2]. Thus, instead of relation (9), we employ

$$k_{\tau }=k_{\tau =0}=const.$$

(22)

This is equivalent to that we increase parameter k multiplying it by $e^{a\tau }$ in DDE (3) and in MP (Table 1). Note that keeping $k_{\tau }$ constant we also keep constant the basic reproduction number $R_0\left(\tau \right)=\left(\beta \lambda k_{\tau }\right)/\left(adu\right)$ defined by (4).

The scenario corresponding to the constant parameter k reflects the following extreme situation when a productively infected cell does not die until a certain quantity of the virions is produced and released. For example, the productively infected cell entering the cell cycle could be considered as escaping the cell death before the division cycle is completed. In addition, some cells can keep producing the virions in oscillatory ways without dying for time periods larger then the considered delays [34]. Hence, the case of a constant k is instructive. We would prefer to keep it in our consideration.

This is an alternative approach to study dependence of the viral infection process on the time delay which is more appropriate for large $\tau$. Dependance of the extinction probability on $\tau$ for different $V_0$ is plotted in Figure 3b. It demonstrates that $P_i$ is independent of $\tau$ in this approach at least within the accuracy of computation.

From here an important and non-obvious result follows: the probability of infection development $P_i$ depends on parameters of a stochastic system combined in $R_0$ (4) which rigorously speaking is the basic reproduction number of the equivalent deterministic system.

To study statistical characteristics: the average trajectories and deviation caused by random fluctuations we should note that the current PDF of species populations can be significantly asymmetric and have more than one peak, i.e., strongly non-Gaussian. This is clearly seen from Figure 4 where examples of histograms of uninfected cells numbers are shown for selected instances.

For example in the second scenario $k_{\tau }=const$ the standard deviation exceeds the mean value of uninfected cells in the interval $18.18<t<18.74$ d. In this case it is more correct to deal with quantiles and median than with mean value and standard deviation. First of all it has sense to consider statistics of extinct and non-extinct realizations separetely. In the extinct realizations, the number of uninfected cells is close to the initial number of cells $x_0={10}^6$ that gives a high peak visible in both plots in Figure 4.

The plots of the current median and interdecimal range (IDR), for all the species in non-extinct realizations are shown in Figure 5 and Figure 6 for $\tau =1$ d and different number of initial number of virions $V_0$. The coloured patches indicate the regions where 80% of all realizations are located. The mean values, $\overline{X},\overline{Y},\overline{V}$, are also plotted by dot-dashed curves. They are close to the median lines in the most parts. The maximum discrepancy between mean, median and deterministic trajectories takes place in the time interval of the fast decrease of the uninfected cells.

The plots show that the IDR regions width remains approximately the same for various number of $V_0$ getting slightly thinner with the growth of $V_0$. Observe also that the smaller number of virions the greater the mean/median number of virions in non-extinct realizations exceeds the deterministic values. Furthermore, the peak of mean/median number of virions in non-extinct realizations appears slightly earlier than peak in the deterministic solution: the smaller number of $V_0$ the greater is the shift.

4. The Hybrid Stochastic Model

The direct numerical simulation is rather time consuming because of a large number of species to be accounted: $X=O\left(X_0\right)\sim {10}^6,V_{max}>X_0$. This causes the time steps between events to be very small. To compute statistical characteristics of the infection process, it is required to run a huge number of realizations especially for evaluation of the infection/extinction probability. Therefore one needs to consider methods to accelerate the computation process.

Plots in Figure 2a show that the randomness in realizations appears mainly due to the relatively high fluctuations at the beginning of viral dynamics while the numbers of virions and infected cells are not large. The curves become rather smooth at a later time, the numbers of all dynamic participants are large (several orders), and the fluid dynamics limit works rather well. This dynamics is typical for an infection process with small initial number of virions ($V_0\ll X_0$): see [16] for viral dynamics and [35,36,37,38] for an epidemic outbreak. Therefore, we can employ the hybrid modelling approach developed in [16,35,37] in which the stochastic dynamics is to be split into two main phases: a genuine stochastic and quasi-deterministic ones.

Phase 1 is the time interval, at which the numbers of virions and infected cells are not large, and hence, the stochasticity of interaction between species has to be taken into account. The number of uninfected cells is large enough and its randomness is not essential in all the phases. Its global variation remains also small for long time after infection starts. The last fact enabled the authors of [16] (similar to what is done in [35,36,37,38]) to simplify description of the earlier phase model by excluding the uninfected cells dynamics from the total process approximating the number of uninfected cell by its initial value: $X\approx x_0$.

In a developed infection process, populations $Y,V$ reach large numbers at a certain instant $t_{\star }$. Thus all the populations become large enough that the total process can be switched to a deterministic behaviour described by DDE (3) (phase 2).

Here we propose a modified method for phase 1 in which the number of uninfected cells $x\left(t\right)$ is not approximated by the constant $x_0$ but can vary obeying a differential equation. This makes the algorithm more flexible as at the switching point the number of uninfected cells can be differ significantly from it initial value $x_0$ (remaining $X\gg 1$).

This coupled deterministic-stochastic process can be described as follows

$$\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}=\lambda -dx\left(t\right)-\beta V\left(t\right)x\left(t\right)\phantom{x-\beta V\left(t\right)x\left(t\right)}$$

(23)

Process (

Table 2. Hybrid Markov Process with Delay.

#	Transition	Rate	Description
1	$Y\left(t\right)\to Y\left(t\right)+1$	$\beta x\left(t\right)V\left(t\right)$	$\mathrm{Cell}\mathrm{infection}$
2	$Y\left(t\right)\to Y\left(t\right)-1$	$aY\left(t\right)$	$\mathrm{Infected}\mathrm{cell}\mathrm{death}$
3	$V\left(t\right)\to V\left(t\right)+1$	$k_{\tau }Y(t-\tau )$	$\mathrm{Virion}\mathrm{birth}$
4	$V\left(t\right)\to V\left(t\right)-1$	$uV\left(t\right)$	$\mathrm{Virion}\mathrm{death}$

) is a reduced MP with delay. Differential Equation (23) has a stochastic parameter $V\left(t\right)$.

In this approach switching from simplified stochastic process (23) to deterministic process (3) occures at the time $t=t_{\star }$ at which the following weaker condition should be fulfilled:

$$min\{Y,V\}\ge X_{\star }\gg 1,$$

(24)

and, therefore, $x\left(t_{\star }\right)$ can be essentially different from $x_0$. Here $X_{\star }$ is a threshold—the minimum number of a component at which its dynamics can be accurately approximated by a deterministic model. For example, in the computations below, we use $X_{\star }={10}^3$.

The initial conditions for process (23) are similar to (10):

$$x\left(0\right)=x_0=\lambda /d,Y\left(0\right)=0,V\left(0\right)=V_0,V_0\in {\mathbb{Z}}^{+}.$$

(25)

The history is also described by (12).

The convergence of MP (23)–(25)–(12) to the deterministic process described by (3)–(6)–(7) is similar to that for MP (Table 1)–(10)–(12).

Stochastic process (Table 2) has been implemented in the same way as MP (Table 1) and described in Algorithm 2. The state vector becomes $\mathit{X}={[Y\left(t\right),V\left(t\right)]}^T$ and the matrix of events has the size $N\times M=2\times 4$:

$$\mathit{T}=\left[\begin{array}{cccc}+1& -1& \phantom{+}0& \phantom{+}0\\ \phantom{+}0& \phantom{+}0& +1& -1\end{array}\right].$$

(26)

The augmented vector for this process is $\check{\mathit{X}}=[Y\left(t\right),V\left(t\right),Y(t-\tau )]$. The propensity vector $\mathit{\nu }={[{\nu }_1,\dots ,{\nu }_M]}^T$ now depends on x and $\check{\mathit{X}}$. The dependence is indicated in Table 2. All other vectors $\mathit{J}={[J_1,\dots ,J_M]}^T$, $\mathit{R}={[R_1,\dots ,R_M]}^T$, $\mathit{t}={[t_1,\dots ,t_{M+1}]}^T$ used in Algorithm 1 are the same with the new $M=4$.

The two-step predictor-corrector method is implemented to integrate ODE (23)

$$\begin{array}{cc}\hfill \tilde{x}& =x_i+\Delta t_i\left(\lambda -d{x}_i-\beta V_i{x}_i\right)\hfill \end{array}$$

(27a)

$$\begin{array}{cc}\hfill x_{i+1}& =x_i+\frac{1}{2}\Delta t_i\left(2\lambda -d{x}_i-d\tilde{x}-\beta V_i{x}_i-\beta V_{i+1}\tilde{x}\right)\hfill \end{array}$$

(27b)

where $\Delta t_i=t_{i+1}-t_i$ is the time interval between two subsequent events in process (Table 2); $x_i=x\left(t_i\right),V_i=V\left(t_i\right)$.

Algorithm 2 Implementation of the hybrid Markov process with delay.

Set the current time $t=0$ and the final time of the process $t_f$; Set $x\left(0\right)=x_0$ and initialise the augmented state vector $\check{\mathit{X}}={[0,V_0,0]}^T$ and vector $\mathit{J}$: $J_m=0,m=1,\dots ,M$; Compute the propensity vector $\mathit{\nu }(x,\check{\mathit{X}})$ using values of propensity indicated in Table 2; Generate M random numbers $r_m$ and compute vector $\mathit{R}$: ${\displaystyle R_m=-ln\left(r_m\right),m=1,\dots ,M}$; Compute the time steps to the next event in all the transitions: ${\displaystyle \Delta t_m=R_m/{\nu }_m,m=1,\dots ,M}$ and set $\Delta t_{M+1}=\infty$; Allocate arrays ${\mathit{t}}^{\mathrm{Y}}$ and ${\mathit{y}}^{\mathrm{Y}}$, set $t_1^{\mathrm{Y}}=0$, $y_1^{\mathrm{Y}}=0$; set counters $i_1=1$, $i_2=0$ for these arrays; Find transition p with the minimal time step: $\Delta t_p=min\{\Delta t_1,\dots ,\Delta t_{M+1}\}$; Update vector $\mathit{J}$: $J_m=J_m+{\nu }_m\Delta t_p$, $m=1,\dots ,M$; Update the current time $t=t+\Delta t_p$; If$p<M+1$then (a) set $\Delta x=\lambda -dx-\beta {\check{X}}_{N}x$ (recall ${\check{X}}_N=V,N=2$) (b) set $\tilde{x}=x+\Delta t_p\Delta x$ (c) update the state vector $\check{\mathit{X}}$: ${\check{X}}_n={\check{X}}_n+T_{np},n=1,2$ (d) set $x=x+\frac{1}{2}\Delta t_p\left(\Delta x+(\lambda -d\tilde{x}-\beta {\check{X}}_N\tilde{x})\right)$ (e) set $J_p=0$ (f) if$p=1$ or $p=2$ (transition in which Y has been changed) then set $i_2=i_2+1$, $Y^{\mathrm{Y}}\left(i_2\right)={\check{X}}_1$, $t^{\mathrm{Y}}\left(i_2\right)=t+\tau$; If$p=M+1$then (a) update the state vector $\check{\mathit{X}}$: ${\check{X}}_{M+1}=Y_{i_1}$ (b) set $i_1=i_1+1$; If$i_2<i_1$ and ${\check{X}}_1={\check{X}}_2=0$ then flag the realisation as extinct and terminate the computation. Update the propensity vector $\mathit{\nu }=\mathit{\nu }(x,\check{\mathit{X}})$; If$p\le M$then generate random number r and compute $R_p=-ln\left(r\right)$; Update the time step vector $\Delta t_m=(R_m-J_m)/{\nu }_m$, $m=1,\dots ,M$; If$i_2\ge i_1$then set $\Delta t_{M+1}=t^Y\left(i_1\right)-t$ else set $\Delta t_{M+1}=\infty$; Store the current state and time; If$t\ge t_f$then terminate the computation. If$min\{{\check{X}}_1,{\check{X}}_2\}\ge X_{\star }$then set $t_{\star }=t$ and switch to phase 2 otherwise go to 7.

If the process has not died out, the computation is continued till time $t_{\star }$ when condition (24) is satisfied. For $t\ge t_{\star }$ the process can be switched to a deterministic one described by DDE (3) with the initial conditions

$$x=x\left(t_{\star }\right),y=Y\left(t_{\star }\right),v=V\left(t_{\star }\right)$$

(28)

at time $t=t_{\star }$ where $x\left(t_{\star }\right),Y\left(t_{\star }\right),V\left(t_{\star }\right)$ are the final results of computation (23).

For DDE (3) we have to define also the history $y\left(t\right),t\in [t_{\star }-\tau ,t_{\star })$. It is complicated to employ the computed stochastic piece-wise function $Y\left(t\right)$ as it contains too many parameters and its use requires elaborate search within intervals between its discontinuities. Instead we employ an approximate solution of DDEs (3) valid at early stage of viral dynamics with $x\left(t\right)=const=x_0$.

Substituting $x=x_0$ into DDEs (3) we reduce it down to two equations

$$\begin{array}{cc}\hfill \dot{y}& =\beta x_{0}v-ay\hfill \end{array}$$

(29a)

$$\begin{array}{cc}\hfill \dot{v}& =k_{\tau }y(t-\tau )-uv\hfill \end{array}$$

(29b)

In contrast to the case $\tau =0$ considered in [16], Equation (29) do not admit a solution in the closed form for all times and should be integrated numerically. Nevertheless, in the interval $[0,\tau ]$ the solution is explicit: here the number of virions decreases exponentially: $v(t\le \tau )=v_{0}exp\{-ut\}$, $y(t\le \tau )\equiv 0$.

DDEs (29) admit an exponential (self-similar) solution

$$v=C{e}^{\alpha t},y=C\frac{\beta x_0}{\alpha +a}{e}^{\alpha t}$$

(30)

where C is a constant determined by $v_0$, $\alpha$ is a positive solution to the equation

$$(\alpha +a)(\alpha +u)=\beta x_0{k}_{\tau }{e}^{-\alpha \tau }.$$

(31)

Dependance of $\alpha$ on $\tau$ for parameters (8) is shown in Figure 7.

Solution (30)–(31) describes an intermediate asymptotic behaviour in the interval between the initial decreasing of the virions number and approaching the infection peak. This asymptotics is clearly seen in Figure 1. Almost linear parts of the curves for $y\left(t\right)$ and $v\left(t\right)$ plotted in logarithmic scale indicate an exponential growing function. An almost exponential infection growth is also seen in the plot of a non-extinct realization in Figure 2a.

Therefore we approximate it by an exponential function utilising the self-similar solution (30) at stage 2 by selecting time $t_{\star }$ such that the instant of the history beginning, $t_{\star }-\tau$, belongs to stage 2. To this aim we solve Equation (31) with respect to parameter $\alpha$ and approximate function $y\left(t\right)$, necessary for the history, as

$$y\left(t\right)=Y\left(t_{\star }\right)e^{\alpha (t-t_{\star })},t\in (t_{\star }-\tau ,t_{\star }).$$

(32)

Computations show that the use of the hybrid model accelerates the process by several orders. For example, the full model with $x_0={10}^6,v_0=8,\tau =0.4,t_{final}=20$ d needs 32.2 s of the CPU time on Intel© Core™ i7-7500U CPU 2.70GHz×2 whereas the hybrid version is computed in 3 ms on average.

An example of several non-extinct realizations computed by the hybrid method for $V_0=100$, $\tau =1$ d and $X_{\star }={10}^3$ is shown in Figure 8. Deterministic part of the curves computed by integration DDEs (3) (phase 2) are drawn by lighter colours. Switching from phase 1 to phase 2 is indicated by circles on the curves for $Y\left(t\right)$ and $V\left(t\right)$. In these realizations, the switching occurs when $Y=min\{Y,V\}=X_{\star }={10}^3$. Value of V in these realizations is greater: $V\left(t_{\ast }\right)\approx 5\times {10}^3$. Thus, $t_{\star }$ varies for different realizations.

Figure 8 illustrates that the stochastic dynamics in phase 1 causes the scattering of the deterministic curves depicting the infection dynamics in phase 2.

The statistical characteristics, median and interdecimal range, for non-extinct realizations obtained by the hybrid method are shown in Figure 9 for all populations for $V_0=20$ and $\tau =1$ d. The notations are the same as in Figure 5. Here for comparison, the same statistical characteristic obtained using the direct simulations (see Figure 5b) are plotted by dash-dot (mean) and thin solid (IDR) lines by a darker colour. Observe that both models give the same statistical characteristics with a high accuracy.

The infection probability $P_i$ computed by the hybrid method is shown in Figure 3 by open black circles and black dashed lines. Here it is also seen that the hybrid method provides high enough accuracy for computation of this important characteristic of infection process.

Another important characteristics of infection process is the time of infection development that can be defined as the time lag between first virions arrived and peak of the infection. It is a challenging problem to determine peak of a stochastic process as it can have high local fluctuations shifting the global maximum from expected position. The paper [35] employs the Gaussian processes to fit the model. The hybrid model allows computation the maximum of the number of infected cells, $Y_{max}$, or the number of virions, $V_{max}$, easily as the peak is located on the smooth deterministic part of the curve (phase 2). The mean value and standard deviation (SD) of $Y_{max}$ and $V_{max}$ versus time delay for different initial number of virions $V_0$ are shown in Figure 10. The averaging is performed over all non-extinct realizations. For comparison, analogous curves obtained by deterministic modelling are indicated as well. The second scenario: $k_{\tau },R_0=const$ is considered in these simulations. In the first scenario: $k_{\tau },R_0\propto e^{-a\tau }$ (9) we obtain very large peak time up to ${10}^3$ d.

Observe that the mean peak is shifted to earlier time for small number of initial virions (1 and 10) and not noticeable for larger numbers (${10}^2$ and ${10}^3$). Especially this shift is large for $V_0=1$, so that the corresponding curve almost coincides with that for $V_0=10$, therefore the later one is drawn by the dashed line. The coefficient of variation (SD over mean) is about 5%–6% for $V_0=1,10,{10}^2$ and about 1% for $V={10}^3$. Therefore the SD patch is very thin and not visible for $V_0={10}^3$ in this scale.

5. Conclusions

In this work, a stochastic viral dynamics model with the time lag between virions production and infection of cells is developed on the base of a Markov process with a time delay. The model has been computationally implemented and studied numerically. The model provides a useful tool to calculate statistical characteristics of virus infection dynamics.

The key statistical characteristics of a virus infection process has been computed: variation of the mean, median and interdecimal range (IDR) for all variables in time, the infection and extinction probabilities and time of the infection development. The dependence of the infection dynamics on the time delay between cell infection and virus progeny production has been studied. Two approaches to incorporate time delays have been tested: (a) the fixed parameters approach in which the reproduction number exponentially decays with the growth of the time delay, and (b) the constant reproduction number approach in which the parameters of the model are adjusted to preserve a constant reproduction number. It is shown that with the second approach, the extinction probability is independent of time delay but the infection peak-time grows non-linearly with the increase of the time delay and differs from the peak-time predicted by the deterministic model for a small number of initial virions. It is also shown that for a small number of initial virions, the number of virions during the infection process, averaged over the non-extinct realizations, exceeds that number calculated via the deterministic model.

A novel and fast computational algorithm to simulate the viral dynamics based on a Markov process with time delay has been proposed, implemented and compared with the full stochastic MP model. In this hybrid model, the dynamics of the components with large numbers of state variables is computed by integration of the ODE/DDE. This essentially accelerates the simulation and computation of the process statistics parameters. It is shown numerically that this hybrid modelling algorithm provides appropriate accuracy in computation of the statistical parameters.

In subsequent work, the full and the hybrid scheme has to be generalized by accounting for a greater number of interacting components (like in the model described in [39]), for a distributed delay, and for more than one delay. This will enable to develop more realistic virus infection models that shall help to better understand regulation and sensitivity of the underlying biological processes.