Stochastic investigation of HIV infection and the emergence of drug resistance

Damilola Olabode; Libin Rong; Xueying Wang; Damilola Olabode; Libin Rong; Xueying Wang

doi:10.3934/mbe.2022054

Mathematical Biosciences and Engineering

2022, Volume 19, Issue 2: 1174-1194. doi: 10.3934/mbe.2022054

Previous Article Next Article

Research article Special Issues

Stochastic investigation of HIV infection and the emergence of drug resistance

1.
Department of Mathematics and Statistics, Washington State University, Pullman, WA 99164, USA
2.
Department of Mathematics, University of Florida, Gainesville, FL 32611, USA

Academic Editor: Yang Kuang

Received: 12 January 2021 Accepted: 22 October 2021 Published: 30 November 2021

Drug-resistant HIV-1 has caused a growing concern in clinic and public health. Although combination antiretroviral therapy can contribute massively to the suppression of viral loads in patients with HIV-1, it cannot lead to viral eradication. Continuing viral replication during sub-optimal therapy (due to poor adherence or other reasons) may lead to the accumulation of drug resistance mutations, resulting in an increased risk of disease progression. Many studies also suggest that events occurring during the early stage of HIV-1 infection (i.e., the first few hours to days following HIV exposure) may determine whether the infection can be successfully established. However, the numbers of infected cells and viruses during the early stage are extremely low and stochasticity may play a critical role in dictating the fate of infection. In this paper, we use stochastic models to investigate viral infection and the emergence of drug resistance of HIV-1. The stochastic model is formulated by a continuous-time Markov chain (CTMC), which is derived based on an ordinary differential equation model proposed by Kitayimbwa et al. that includes both forward and backward mutations. An analytic estimate of the probability of the clearance of HIV infection of the CTMC model near the infection-free equilibrium is obtained by a multitype branching process approximation. The analytical predictions are validated by numerical simulations. Unlike the deterministic dynamics where the basic reproduction number $\mathcal{R}_0$ serves as a sharp threshold parameter (i.e., the disease dies out if $\mathcal{R}_0 < 1$ and persists if $\mathcal{R}_0 > 1$ ), the stochastic models indicate that there is always a positive probability for HIV infection to be eradicated in patients. In the presence of antiretroviral therapy, our results show that the chance of clearance of the infection tends to increase although drug resistance is likely to emerge.

Keywords:

Citation: Damilola Olabode, Libin Rong, Xueying Wang. Stochastic investigation of HIV infection and the emergence of drug resistance[J]. Mathematical Biosciences and Engineering, 2022, 19(2): 1174-1194. doi: 10.3934/mbe.2022054

Related Papers:

[1]	Wenshuang Li, Shaojian Cai, Xuanpei Zhai, Jianming Ou, Kuicheng Zheng, Fengying Wei, Xuerong Mao . Transmission dynamics of symptom-dependent HIV/AIDS models. Mathematical Biosciences and Engineering, 2024, 21(2): 1819-1843. doi: 10.3934/mbe.2024079
[2]	Damilola Olabode, Jordan Culp, Allison Fisher, Angela Tower, Dylan Hull-Nye, Xueying Wang . Deterministic and stochastic models for the epidemic dynamics of COVID-19 in Wuhan, China. Mathematical Biosciences and Engineering, 2021, 18(1): 950-967. doi: 10.3934/mbe.2021050
[3]	Shengqiang Liu, Lin Wang . Global stability of an HIV-1 model with distributed intracellular delays and a combination therapy. Mathematical Biosciences and Engineering, 2010, 7(3): 675-685. doi: 10.3934/mbe.2010.7.675
[4]	Gesham Magombedze, Winston Garira, Eddie Mwenje . Modelling the immunopathogenesis of HIV-1 infection and the effect of multidrug therapy: The role of fusion inhibitors in HAART. Mathematical Biosciences and Engineering, 2008, 5(3): 485-504. doi: 10.3934/mbe.2008.5.485
[5]	Churni Gupta, Necibe Tuncer, Maia Martcheva . A network immuno-epidemiological model of HIV and opioid epidemics. Mathematical Biosciences and Engineering, 2023, 20(2): 4040-4068. doi: 10.3934/mbe.2023189
[6]	Linda J. S. Allen, Vrushali A. Bokil . Stochastic models for competing species with a shared pathogen. Mathematical Biosciences and Engineering, 2012, 9(3): 461-485. doi: 10.3934/mbe.2012.9.461
[7]	Yan Wang, Tingting Zhao, Jun Liu . Viral dynamics of an HIV stochastic model with cell-to-cell infection, CTL immune response and distributed delays. Mathematical Biosciences and Engineering, 2019, 16(6): 7126-7154. doi: 10.3934/mbe.2019358
[8]	Divine Wanduku . A nonlinear multi-population behavioral model to assess the roles of education campaigns, random supply of aids, and delayed ART treatment in HIV/AIDS epidemics. Mathematical Biosciences and Engineering, 2020, 17(6): 6791-6837. doi: 10.3934/mbe.2020354
[9]	Zhuolin Qu, Benjamin H. McMahon, Douglas J. Perkins, James M. Hyman . Staged progression epidemic models for the transmission of invasive nontyphoidal Salmonella (iNTS) with treatment. Mathematical Biosciences and Engineering, 2021, 18(2): 1529-1549. doi: 10.3934/mbe.2021079
[10]	Xiaodan Sun, Yanni Xiao, Zhihang Peng . Modelling HIV superinfection among men who have sex with men. Mathematical Biosciences and Engineering, 2016, 13(1): 171-191. doi: 10.3934/mbe.2016.13.171

Abstract

1. Introduction

Drug resistance (antimicrobial and anti-neoplastic) is a global concern, threatening the treatment of infectious diseases and cancer based diseases, which results in prolonged illness, disability and death [1]. The World Health Organization (WHO) defines antimicrobial resistance as a microorganism's resistance to an antimicrobial drug that was once able to treat an infection by that microorganism [2]. Microbes include protozoa, fungi, virus and bacteria. Drug resistance can be associated with cases such as clinical misuse, natural occurrence, self medication. Antibiotics (treatment for bacterial infection) drug resistance has been the most common and long standing form of drug resistance as far back as in the 1950s'. Anti-neoplastic resistance (i.e., synonymous with chemotherapy resistance) is the ability of cancer cells to survive and grow despite different anti-cancer therapies [3]. Lots of works have been dedicated to drug resistance in antimicrobial [4-6] and in anti-neoplastic [3,7-9]. Beyond the studies above are some mathematical insights into drug resistance. For instance, in 1998, Blower et al. [10] used mathematical model to suggest control measures for herpes epidemics that would prevent the emergence of substantial levels of antiviral drug resistance. In 2009, Cohen et al. [11] used mathematical models to offer insight into the acquisition and amplification of drug resistance in patients with tuberculosis.

One of the leading causes of death in most sub saharan African countries is the acquired immune deficiency syndrome (AIDS), a chronic, potentially life-threatening condition caused by HIV. Currently, seven classes of drug^* are used to inhibit the life cycle of HIV through the antiretroviral therapy (ART). In the treatment of HIV, ART uses the combination of three or more drugs from two or more of the above classes [13]. The commonly used classes of drug are the nucleoside reverse transcriptase inhibitors (NRTIs) and the protease inhibitors (PIs). Drug combination has proved to suppress the plasma viral loads to below the clinical detection limit. However, a complete suppression of the viral detection is not seen due to either host or viral factors, which in essence decreases the survival rate of the host and increases HIV infection. These factors include imperfect adherence to drug regimen, drug resistance, poor drug absorption (e.g., see [13-17]).

^* The nucleoside reverse transcriptase inhibitors (NRTIs), non-nucleoside reverse transcriptase inhibitors (NNRTIs), protease inhibitors (PIs), fusion inhibitors, CCR5 antagonists, post-attachment inhibitors, and integrase strand transfer inhibitors (INSTIs) [12]

Drug resistance in HIV has been a major factor for treatment failure [18]. Larder et al. [14] defined drug resistance as the ability for the virus to replicate even in the presence of drugs. HIV drug resistance has been found to be as a result of incomplete suppression of viral loads rather than the infection of resistant strain from other people.

Several works have been dedicated to getting insight into the dynamics of HIV in the presence of drug resistance. In 1992, McLean and Nowak [] investigated the competitive interactions between zidovudine-sensitive (an NRTI HIV drug class) and resistant strains of HIV within the context of interaction between CD4 $^+$ T cells and HIV, using a deterministic approach. In 1997, Kirschner and Webb [] studied strategies in a monotherapy treatment of HIV infection in the presence of both sensitive and resistance strains at different levels of CD4 $^+$ T cells counts. In 1997, Nowak et al. [21] provided an analytical approximation for the rate of the emergence of resistant strain in different compartments of the virus population. In 1998, Ribeiro et al. [22] calculated the expected frequency of the drug-resistant strain in patients who have not received treatments. Their result indicated that it is possible that in the onset of treatment in these patients, the resistant strain could be selected. In 1998, Kepler and Perelson [23] suggested that spatial heterogeneity can be a factor in the evolution of drug resistance. In 2005, Murray and Perelson [24] showed the influence of quasi species on the resistance to zidovudine monotherapy. In 2007, Rong et al. [13] analytically derived expressions that specify the condition in which the resistant strain is selected. In 2010, Vaidya et al. [25] studied patients with resistance to enfuvirtide (HIV-1 fusion inhibitor). Their result shows that the rapid replacement of resistant virus during an ART interruption was attributed to the associated fitness cost to the resistance strain. In 2012, Kitayimbwa et al. [13] used a deterministic model with both forward and backward mutations to show the coexistence of both the sensitive and resistant strains, and they derived the conditions for dominance of the viral strain in the presence and absence of ART. Moreno-Gamez et al. [26] employed a simple birth-death process model to study multidrug-resistance. Their results show that imperfect drug penetration can lead to monotherapy at sites of bacterial residence despite treatment with several drugs.

Most of the above works are built upon deterministic models. Only a very few of them has been done via a stochastic approach to study the emergence of HIV-1 drug resistance [26,27]. It is worthy mentioning that there is a body of mathematical models using stochastic formulations to investigate within-host viral HIV dynamics (e.g., [27-36] and the references therein). For instance, Perelson et al. [37] developed a stochastic model of the early HIV-1 infection to probe the earliest virus-host events following virus inoculation. Hill et al. [31] focused on the latent infected cells to predict delays in viral rebound via a stochastic formulation. Conway and her collaborators studied how prophylaxis with antiretroviral drugs can reduce the risk of infection [29] and why some HIV-infected patients experienced post-treatment control [30]. A recent paper combined stochastic and deterministic within-host models to estimate the time of infection [38]. But these models did not consider the emergence of drug resistance during therapy.

The primary goal of this paper is to investigate HIV-1 infection and the emergence of drug resistance. The early stage of the infection (i.e., the first few hours to days following HIV exposure) may determine whether the infection can be successfully established. However, the numbers of infected cells and viruses during the early stage are extremely low and stochasticity may play a critical role in dictating the fate of infection. To that end, we use stochastic modeling to study the probability of the clearance of HIV infection, time to the clearance of HIV infection, and time to the establishment of infection.

The remainder of the article is organized as follows. Section 2 describes the deterministic model proposed in [13,15] and summarizes their results. In Section 3, we develop a continuous-time Markov chain (CTMC) model and employ a multitype branching process theory to approximate the dynamics of the CTMC model near the infection-free equilibrium. Numerical simulations are performed in Section 4. Finally, the paper is concluded with some discussion.

2. Deterministic model and analysis

We consider an ordinary differential equation model that has been described and analyzed in [15]. An extreme case of this model without backward mutation was studied in an earlier work [13]. The system is given as follows:

$\begin{equation} \begin{split} \dfrac{dT}{dt} & = \lambda -dT-{(1-\epsilon_{RT}^{s})} k_sV_sT - {(1-\epsilon_{RT}^{r})} k_rV_rT \\ \dfrac{dT_s}{dt} & = (1-u){(1-\epsilon_{RT}^{s})}k_sV_sT +z{(1-\epsilon_{RT}^{r})}k_rV_rT -\delta T_s \\ \dfrac{dV_s}{dt} & = {(1-\epsilon_{PI}^{s})}N_s \delta T_s - cV_s\\ \dfrac{dT_r}{dt} & = u{(1-\epsilon_{RT}^{s})}k_sV_sT +(1-z){(1-\epsilon_{RT}^{r})}k_rV_rT - \delta T_r \\ \dfrac{dV_r}{dt} & = {(1-\epsilon_{PI}^{r})}N_r \delta T_r - cV_r \end{split} \end{equation}$

(2.1)

where $t$ is the time variable, $T(t)$ is the concentration of susceptible T-cells at time $t$ , $T_s(t)$ is the concentration of productively infected cells by the drug-sensitive virus. $V_s(t)$ denotes the concentration of drug-sensitive virus. $k_s$ is the constant rate at which the uninfected cells are infected by the drug-sensitive strain. $T_r(t)$ , $V_r(t)$ and $k_r$ are the the concentration of productively infected cells, the concentration of virus and the infection rate associated with resistant strain, respectively. The proportions $u$ and $z$ depict the forward mutation (sensitive strain to resistant strain) and the backward mutation (resistant strain to sensitive strain), respectively, with $0 < u\le 1, 0\le z\le 1$ . The $\epsilon_j^i$ 's ( $i = s, r$ and $j = RT, PI$ ) are constant drug efficacies corresponding to the reverse transcriptase (RT) and the protease inhibitor (PI) for the sensitive or resistant strain with $0 \leq \epsilon_{RT}^{s}, \epsilon_{RT}^{r}, \epsilon_{PI}^{s}, \epsilon_{PI}^{r} < 1$ . The definition and the baseline values of the rest of the parameters are given in Table 1. In the absence of treatment (i.e. all the drug efficacies are 0), the model can be used to study whether the infection can be successfully established and whether the mutant virus, due to mutations, preexists before therapy. During treatment (with positive drug efficacies), the model investigates how the drug-resistant virus develops and competes with the wild-type virus.

Table 1. Parameter values and definitions.

Parameter	Definition	Value	Reference
$T_0$	Initial target cell count	$10^6 \ mL^{-1}$	[39]
$d$	Death rate of uninfected cells	$0.01 \ day^{-1}$	[40]
$\lambda$	Recruitment rate of uninfected cells	$10^4 \ mL^{-1} day^{-1}$	[41]
$k_s$	Infection rate of target cells by wild-type virus	$2.4 \times 10^{-8} \ mL \ day^{-1}$	[41]
$k_r$	Infection rate of target cells by drug-resistance virus	$2.0 \times 10^{-8} \ mL \ day^{-1}$	[13]
$u$	Forward mutation rate	$2.24 \times 10^{-5}$	[25]
$z$	Backward mutation rate	$1.73 \times 10^{-5}$	[25]
$\delta$	Death rate of infected cells	$1 \ day^{-1}$	[42]
$N_s$	Burst size of drug-sensitive strain	$3000$	[13]
$N_r$	Burst size of drug resistance strain	$2000$	[13]
$c$	Clearance rate of free virus	$23 \ day^{-1}$	[43]

| Show Table

DownLoad: CSV

2.1. Basic reproduction number

The dynamics of the system were studied by Kitayimbwa et al. in []. The basic reproduction number $\mathcal{R}_0$ is defined as the expected number of secondary cases produced by one primary case in an otherwise susceptible population. Using the next generation method [], it follows from the analysis in [] that $\mathcal{R}_0$ associated with system (2.1) is

$\begin{equation} \begin{split} \mathcal{R}_0 & = \dfrac{\lambda}{2cd} \bigg[ \big((1-u)k_s N_s{(1-\epsilon^s)} + (1-z)k_r N_r{(1-\epsilon^r)} \\ &+ \sqrt{ \bigg((u-1)k_s N_s {(1-\epsilon^s)} + (z-1)k_r N_r {(1-\epsilon^r)}\bigg)^2 - 4 (1-u-z)k_r k_s N_r N_s {(1-\epsilon^s)(1-\epsilon^r)}}\bigg] \end{split} \end{equation}$

(2.2)

Here $(1-\epsilon^s) = (1-\epsilon_{RT}^s)(1-\epsilon_{PI}^s) \text{ and } (1-\epsilon^r) = (1-\epsilon_{RT}^r)(1-\epsilon_{PI}^r)$ . For completeness, the derivation of $\mathcal{R}_0$ is provided in the Appendix.

2.1.1. Case 1: Forward mutation ( $z = 0 \text{ and }0 < u < 1$ )

Suppose $z = 0$ , $0 < u < 1$ , and no implementation of the ART (i.e., all $\epsilon_j^i = 0$ ) for $j = RT$ , $PI$ and $i = s$ , $r$ in model (2.1). In this case, only the forward mutation is allowed and the associated deterministic dynamics have been analyzed in []. In what follows, we summarize their results. There are at most three biologically feasible steady states of the form $\bar{E} = \big(\bar{T}, \bar{T}_s, \bar{V}_s, \bar{T}_r, \bar{V}_r \big)$ . The first one is the infection-free equilibrium (IFE) that always exists, and it is given by

$E_0 = \Big(\frac{ \lambda}{d}, 0, 0, 0, 0\Big).$

The second one is a boundary steady state

$E_r = \bigg(\frac{c}{k_r N_r}, 0, 0, (\mathcal{R}_r -1) \frac{dc}{k_r N_r \delta}, (\mathcal{R}_r -1) \frac{d}{k_r} \bigg)$

which captures that the survival strain is the drug-resistant strain. The third one is an interior steady state, which represents the coexistence of both viral strains and is given by

$E_c = \bigg(\frac{c}{(1-u)k_s N_s}, \tilde{T}_s, \frac{N_s \delta}{c} \tilde{T}_s, \tilde{T}_r, \frac{N_r \delta}{c} \tilde{T}_r \bigg)$

where

$\begin{equation*} \tilde{T}_s = \frac{[(1-u)\sigma - 1][(1-u)\mathcal{R}_s -1]\lambda}{(\sigma -1)(1-u)\mathcal{R}_s \delta}, \quad \tilde{T}_r = \frac{[(1-u)\mathcal{R}_s -1] u \sigma \lambda}{(\sigma -1)(1-u)\mathcal{R}_s \delta}, \quad \sigma = \frac{\mathcal{R}_s}{\mathcal{R}_r}. \end{equation*}$

Here

$\begin{equation*} \mathcal{R}_s = \frac{k_s N_s \lambda}{dc}, \quad \mathcal{R}_r = \frac{k_r N_r \lambda}{dc} \end{equation*}$

are the reproduction numbers of the wild-type and the drug-resistant strain, respectively, by assuming that the system is decoupled.

It is obvious that if $z = 0$ in model (2.2) (i.e., the case of the forward mutation), then the basic reproduction number $\mathcal{R}_0$ is given by

$\mathcal{R}_0 = \max \Big\{ (1-u)\, \mathcal{R}_s, \mathcal{R}_r\Big\}.$

The stability result of the forward mutation obtained in [13] is summarized as follows.

Theorem 2.1 ([13], Propositions 1 and 2]).

$a.$ The IFE $E_0$ is locally asymptotically stable if $\mathcal{R}_0 < 1$ , and it is unstable if $\mathcal{R}_0 > 1$ .

$b.$ The boundary steady state with only drug-resistant virus $E_r$ exists if and only if $\mathcal{R}_r > 1$ . It is locally asymptotically stable if $\mathcal{R}_r > (1-u)\mathcal{R}_s$ and unstable if $\mathcal{R}_r < (1-u)\mathcal{R}_s.$

$c.$ The coexistence steady state $E_c$ exists and is locally asymptotically stable if and only if

$\mathcal{R}_s > \max\left\{\frac{1}{1-u}, \frac{\mathcal{R}_r}{1-u}\right\}.$

2.1.2. Case 2: Forward and backward mutation ( $0 < u, z < 1$ )

Suppose that $0 < u < 1, \, 0 < z < 1, \text{ and }{\epsilon_j^i = 0}$ in model (2.1), then both forward and backward mutations are allowed. The dynamical result obtained in [] is presented below. In this case, we have up to two steady states. One is the IFE $E_0$ and another one is the endemic equilibrium (EE), $E_1 = (T^*, T_s^*, V_s^*, T_r^*, V_r^*)$ , where

$\begin{equation} \begin{split} T^* & = \frac{\lambda}{d(1-u-z)\mathcal{R}_s \mathcal{R}_r} \big[(1-u)\mathcal{R}_s+(1-z)\mathcal{R}_r-\mathcal{R}_0\big], \\ T_s^* & = \frac{\lambda}{\delta \tau} \bigg[ \mathcal{R}_0-\mathcal{R}_r + \frac{z \tau - \mathcal{R}_0}{\mathcal{R}_s (1-u-z)}+\frac{1}{\sigma (1-u-z)}\bigg], \\ V_s^* & = \frac{N_s \delta}{c} T_s^*, \\ T_r^* & = \frac{\lambda}{\tau \delta} \bigg[\frac{\mathcal{R}_0}{\mathcal{R}_r (1-u-z)}+(\mathcal{R}_s - \mathcal{R}_0) - \frac{\sigma(1-u)+u}{(1-u-z)}\bigg], \\ V_r^* & = \frac{N_r \delta}{c} T_r^*, \end{split} \end{equation}$

(2.3)

with

$\sigma = \frac{\mathcal{R}_s}{\mathcal{R}_r}, \text{ and } \tau = \mathcal{R}_s - \mathcal{R}_r.$

The stability result [15], Theorem 1 & Theorem 3] for system (2.1) with the forward and backward mutations is summarized as follows.

Theorem 2.2 ([15], Theorems 1, 3]).

$(1)$ (Global stability) If $\mathcal{R}_0 < 1$ in model (2.2), then the IFE $E_0$ is the unique equilibrium solution and it is globally asymptotically stable in $\mathbb{R}^5_{+}$ .

$(2)$ (Local stability) If $\mathcal{R}_0 > 1$ , system (2.1) has equilibrium solutions: the IFE $E_0$ and the EE $E_1$ . Moreover, the IFE $E_0$ is unstable and the EE $E_1$ is locally asymptotically stable.

This result indicates that $\mathcal{R}_0$ serves as an infection threshold; i.e., if $\mathcal{R}_0 < 1$ , the infection will be cleared, whereas if $\mathcal{R}_0 > 1$ , the infection will become established. It is worthy mentioning that it is straightforward to extend the analysis of the deterministic dynamics to the case where drug treatment is included, and a similar result remains to be valid; i.e., the corresponding basic reproduction number is a threshold parameter of the HIV infection by using the deterministic model (2.1).

2.2. Sensitivity Analysis

To determine the effect that model parameters on $\mathcal{R}_0$ , we perform a sensitivity analysis by computing the normalized sensitivity index $\mathcal{S}$ which is introduced in [45] and defined as follows:

$\begin{equation*} \mathcal{S} = \left. \frac{d \mathcal{R}_0}{dp} \cdot \frac{p}{\mathcal{R}_0} \right|_{p = p_0} \end{equation*}$

where $p$ is the parameter of interest and $p_0$ is the base value of this parameter.

Our results of the sensitivity analysis for the case without and with drug treatment, ART, are summarized in and , respectively. The obtained normalized sensitivity indices are ranked from the largest to the smallest by using the corresponding magnitudes. The sign of the sensitive index indicates the direction of the change in $\mathcal{R}_0$ . More specifically, the positive (negative) sign represents that $\mathcal{R}_0$ will increase (decrease) as the parameter increases. The sensitivity analysis displayed in shows that the parameters $\lambda$ , $k_s$ , $N_s$ , $c$ , and $d$ are the most influential parameters in terms of the relative change in $\mathcal{R}_0$ . This suggests that the basic reproduction number is most sensitive to the wild-type strain parameters rather than the mutant strain parameters in the absence of ART. Thus, the contribution to the basic reproduction number from de novo mutation is minor. This is reasonable as the mutation rate is very small and drug treatment has not been considered in this basic reproduction number. In contrast, with the presence of drug treatment, shows that $\lambda$ , $k_r$ , $N_r$ , $c$ , and $d$ are the most influential parameters for $\mathcal{R}_0$ . This indicates that the drug-resistant strain is likely to emerge and dominate the virus population when it is selected by the drug pressure (i.e. drug treatment is more effective in inhibiting the wild-type strain). The dynamics of the resistant strain will be illustrated later in a numerical simulation.

Table 2. Relative sensitivity analysis of parameters in the absence of ART with the wild-type strain dominating the viral dynamics. Here

$\epsilon_{RT}^s = \epsilon_{RT}^r = \epsilon_{PI}^s = \epsilon_{PI}^r = 0$ and the rest of parameter values are provided in Table 1.

Parameter	$\mathcal{S}$
$\lambda$	$1$
$k_s$	$1$
$N_s$	$1$
$d$	$-1$
$c$	$-1$
$u$	$-2.24 \times 10^{-5}$
$N_r$	$1.088 \times 10^{-9}$
$k_r$	$1.076 \times 10^{-9}$
$z$	$4.844 \times 10^{-10}$

| Show Table

DownLoad: CSV

Table 3. Relative sensitivity analysis of parameters under ART with the resistant strain dominating the viral dynamics. Here

$\epsilon_{RT}^s = 0.5, \epsilon_{RT}^r = 0.02, \epsilon_{PI}^s = 0.1, \text{and }\epsilon_{PI}^r = 0.01$ and the rest of parameter values are the same as in Table 1.

Parameter	$\mathcal{S}$
$\lambda$	$1$
$k_r$	$1$
$N_r$	$1$
$d$	$-1$
$c$	$-1$
$z$	$-1.73 \times 10^{-5}$
$N_s$	$1.203 \times 10^{-8}$
$k_s$	$1.173 \times 10^{-8}$
$u$	$1.959 \times 10^{-9}$

| Show Table

DownLoad: CSV

3. Stochastic models

ODE models typically assume that the sizes of the compartments are large enough that the members are homogeneously mixed, or at least that there is homogeneous mixing in each subgroup if the population is stratified by activity levels. However, at the beginning of an infection, there is a very small number of infected cells or virions and the infection is indeed a stochastic event depending on heterogeneity among cells/virions and patterns of contacts between them. Hence it suggests that the homogeneous mixing at the beginning of an infection may not be a good assumption and the ODE models may be inappropriate when the initial infection is small. In particular, it has been showed that stochasticity may play a significant role when the viral load is low [28,46].

Therefore, we use a continuous-time Markov chain (CTMC) model to study the variability of the disease dynamic, and then we apply the theory of the multitype branching process (e.g., see [47-52] and references therein) to approximate the dynamics of the CTMC model near the IFE.

3.1. Continuous-time Markov Chain (CTMC)

The continuous-time Markov chain (CTMC) is a stochastic process that is continuous in time and discrete in state space. The ODE model (2.1) serves as a framework for formulation of the CTMC model. For simplicity, the same notation is used for the variables as in the ODE model. Let ${X}(t) = (T(t), T_s(t), V_s(t), T_r(t), V_r(t))$ denote the discrete random vector for the states, where for each $t\in [0, \infty)$ ,

$\begin{equation*} T(t), T_s(t), V_s(t), T_r(t), V_r(t)\in \big\{0, 1, 2, 3, \cdots\big\}. \end{equation*}$

Let $\Delta {X}(t) = {X}(t+\Delta t) - {X}(t)$ . In the CTMC model, there are twelve independent Poisson processes and the corresponding infinitesimal transition probabilities are defined in Table 4.

Table 4. State transitions and associated rates for the CTMC model.

Event	Description	Transition $\Delta {X}(t)$	Transition rate, $r_i$
1	Recruitment of $T$	$(1, 0, 0, 0, 0)$	$\lambda$
2	Death of $T$	$(-1, 0, 0, 0, 0)$	$dT$
3	Infection due to sensitive strain	$(-1, 1, 0, 0, 0)$	$(1-u){(1-\epsilon_{RT}^s)}k_s V_s T$
4	Forward mutation	$(-1, 0, 0, 1, 0)$	$u {(1-\epsilon_{RT}^s)}k_s V_s T$
5	Backward mutation	$(-1, 1, 0, 0, 0)$	$z{(1-\epsilon_{RT}^r)}k_rV_r T$
6	Infection caused by resistant strain	$(-1, 0, 0, 1, 0)$	$(1-z){(1-\epsilon_{RT}^r)}k_rV_r T$
7	Death of $T_s$	$(0, -1, 0, 0, 0)$	$\delta T_s$
8	Production of $V_s$	$(0, 0, 1, 0, 0)$	${(1-\epsilon_{PI}^s)}N_s \delta T_s$
9	Clearance of $V_s$	$(0, 0, -1, 0, 0)$	$c V_s$
10	Death of $T_s$	$(0, 0, 0, -1, 0)$	$\delta T_r$
11	Production of $V_r$	$(0, 0, 0, 0, 1)$	${(1-\epsilon_{PI}^r)}N_r \delta T_r$
12	Clearance of $V_r$	$(0, 0, 0, 0, -1)$	$c V_r$

| Show Table

DownLoad: CSV

3.2. Branching process approximation

A multi-type branching process approximation is applied to the four infectious classes, i.e., $T_s, V_s, T_r, V_r$ , in the CTMC model at the IFE, where $T = \lambda/d$ . The multi-type branching process is linear (in terms of infectious classes) and time homogeneous with independent births and deaths. Thus, we are able to define the probability generating function (pgfs) for the birth and death of each infected class. As a result, we can use this to calculate the probability of the establishment of infection.

In general, let $\{Z(t) : t \in [0, \infty)\}$ be a collection of discrete valued vector of random variables where $Z(t) = (Z_1(t), Z_2(t), \cdots, Z_n(t))^T$ represents the disease classes. We assume that the birth of an offspring and the number of the offspring that is produced by a certain type are independent. We also assume that the individuals of type $j$ follows an identical offspring pgf, i.e., they are independent and identically distributed. Let $X_j(0) = \delta_{jk}$ , where $\delta_{jk} = 1 \text{ if } i = j$ and zero otherwise, for $1\leq j, k \leq n$ . The offspring pgf for the infectious class of type $j$ is a function $f_j:[0, 1]^n \rightarrow [0, 1]$ with

$\begin{align*} f_j(s_1, \cdots, s_n) = \sum\limits_{k_n = 0}^{\infty} \cdots \sum\limits_{k_1 = 0}^{\infty} p_j(k_1, \cdots, k_n) s_1^{k_1} \cdots s_n^{k_n} \end{align*}$

where $p_j(k_1, k_2, \cdots, k_n)$ is the probability that parent of type $j$ has $k_i$ offsprings of type $i$ for $i = 1, 2, \cdots, n$ .

In the case of our model, there are four disease components $n = 4$ and

$Z = (Z_1, Z_2, Z_3, Z_4) = (T_S, V_s, T_r, V_r).$

The associated multitype branching process consists of ten independent events and the corresponding transition rates are obtained directly from the rates in the linearized system of the ODE model (2.1) at the IFE, which is summarized in . Using this, we derive the offspring pgf in terms of $T_S, V_s, T_r, V_r$ and they are given by

$\begin{align*} f_1(s_1, s_2, s_3, s_4)& = \frac{{(1-\epsilon_{PI}^s)}N_s \delta s_1s_2 +\delta }{{(1-\epsilon_{PI}^s)}N_s \delta + \delta}, \\ f_2(s_1, s_2, s_3, s_4)& = \dfrac{\dfrac{(1-u){(1-\epsilon_{RT}^s)}k_s \lambda}{d} s_1 s_2 +c + u{(1-\epsilon_{RT}^s)}k_s \frac{\lambda}{d} s_2 s_3}{{(1-\epsilon_{RT}^s)}k_s \frac{\lambda}{d} + c}, \\ f_3(s_1, s_2, s_3, s_4)& = \frac{{(1-\epsilon_{PI}^r)}N_r \delta s_3s_4 +\delta }{{(1-\epsilon_{PI}^r)}N_r \delta + \delta}, \\ f_4(s_1, s_2, s_3, s_4)& = \frac{\dfrac{(1-z){(1-\epsilon_{RT}^r)}k_r \lambda}{d} s_3 s_4 +c + z{(1-\epsilon_{RT}^r)}k_r \frac{\lambda}{d} s_1 s_4}{{(1-\epsilon_{RT}^r)}k_r \frac{\lambda}{d} + c}. \end{align*}$

Table 5. Transition probabilities

$\hat{r}_i\Delta t+o(\Delta t)$ for branching process (BP) approximation.

Event	Description	$\Delta Z(t)$	BP rate $\hat{r}_i$
1	Infection caused by sensitive strain	$(1, 0, 0, 0)$	$(1-u){(1-\epsilon_{RT}^s)}k_s V_s(t) \frac{\lambda}{d}$
2	Forward mutation	$(0, 0, 1, 0)$	$u{(1-\epsilon_{RT}^s)}k_sV_s(t) \frac{\lambda}{d}$
3	Backward mutation	$(1, 0, 0, 0)$	$z {(1-\epsilon_{RT}^r)}k_r V_r(t) \frac{\lambda}{d}$
4	Infection caused by resistant strain	$(0, 0, 1, 0)$	$(1-z){(1-\epsilon_{RT}^r)}k_rV_r(t)\frac{\lambda}{d}$
5	Death of $T_s$	$(-1, 0, 0, 0)$	$\delta T_s(t)$
6	Production of $V_s$	$(0, 1, 0, 0)$	${(1-\epsilon_{PI}^s)}N_s \delta T_s(t)$
7	Clearance of $V_s$	$(0, -1, 0, 0)$	$c V_s(t)$
8	Death of $T_r$	$(0, 0, -1, 0)$	$\delta T_r(t)$
9	Production of $V_r$	$(0, 0, 0, 1)$	${(1-\epsilon_{PI}^r)}N_r \delta T_r(t)$
10	Clearance of $V_r$	$(0, 0, 0, -1)$	$c V_r(t)$

| Show Table

DownLoad: CSV

By the theory of branching process [], Theorem 1.2], the probability of the clearance of HIV infection of the above system can be determined by solving the fixed points $f_i(q_1, q_2, q_3, q_4) = q_i\ (1\le i\le 4)$ on $[0, 1]^4$ . We know that there is always going to be a fixed point at $(1, 1, 1, 1)$ . We are interested in the minimal fixed point $(q_1, q_2, q_3, q_4)$ in which $q_i\in [0, 1]$ for $1\le i\le 4$ . In general, there is no closed-form expression for the minimal fixed point. So we have to numerically solve $q_i$ ( $1\le i\le 4$ ).

Stochastic threshold:

Let $M = (m_{kj})$ be the expectation matrix of the offspring pgfs. $M$ is a nonnegative $4 \times 4$ matrix whose entry $m_{kj}$ gives the expected number of offsprings in type $k$ produced by an infectious individual in type $j$ ; that is,

$m_{kj} = \left. \frac{\partial f_k}{\partial s_j} \right|_{(s_1, s_2, s_3, s_4) = (1, 1, 1, 1)}.$

Then

$M = (m_{kj}) = \begin{bmatrix} \frac{{(1-\epsilon_{PI}^s)}N_s \delta}{{(1-\epsilon_{PI}^s)}N_s \delta + \delta} & \frac{{(1-\epsilon_{PI}^s)}N_s \delta}{{(1-\epsilon_{PI}^s)}N_s \delta + \delta} &0 &0\\ \frac{\frac{(1-u){(1-\epsilon_{RT}^s)}k_s \lambda}{d}}{\frac{{(1-\epsilon_{RT}^s)}k_s \lambda}{d}+c} & \frac{{(1-\epsilon_{RT}^s)}k_s \frac{\lambda}{d}}{\frac{{(1-\epsilon_{RT}^s)}k_s \lambda}{d}+c}&\frac{\frac{u {(1-\epsilon_{RT}^s)}k_s \lambda}{d}}{\frac{{(1-\epsilon_{RT}^s)}k_s \lambda}{d}+c}&0\\ 0&0& \frac{{(1-\epsilon_{PI}^r)}N_r \delta}{{(1-\epsilon_{PI}^r)}N_r \delta + \delta} & \frac{{(1-\epsilon_{PI}^r)}N_r \delta}{{(1-\epsilon_{PI}^r)}N_r \delta + \delta}\\ \frac{z{(1-\epsilon_{RT}^r)}k_r\frac{\lambda}{d}}{{(1-\epsilon_{RT}^r)}k_r \frac{\lambda}{d}+c} & 0& \frac{(1-z){(1-\epsilon_{RT}^r)}k_r\frac{\lambda}{d}}{{(1-\epsilon_{RT}^r)}k_r \frac{\lambda}{d}+c}& \frac{{(1-\epsilon_{RT}^r)}k_r\frac{\lambda}{d}}{{(1-\epsilon_{RT}^r)}k_r \frac{\lambda}{d}+c} \end{bmatrix}.$

The expectation matrix $M$ is irreducible since the corresponding digraph is strongly connected. It is easy to verify that $F$ and $V$ defined in the appendix are non-negative and non-singular $M$ -matrices (see [53]). By the threshold theorem [54], we obtain the following relationship between the deterministic threshold and the stochastic threshold.

$\begin{equation} \mathcal{R}_0 < 1 \, ( = 1, > 1) \ \ \text{ if and only if }\ \ \rho(M) < 1 \, ( = 1, > 1). \end{equation}$

(3.1)

In the subcritical and critical cases where $\rho(M) \le 1$ (i.e., $\mathcal{R}_0 \le 1$ ), the probability of the ultimate clearance of HIV infection is 1; i.e.,

$P_{e x t} = \lim\limits_{t \rightarrow \infty} {\rm{Prob}}\{\underline{\mathrm{X}}(t) = \overrightarrow{0}\} = 1$

In the supercritical case where $\rho(M) > 1$ (i.e., $\mathcal{R}_0 > 1$ ), the corresponding clearance of HIV infection probability is determined by the minimal fixed point of $f_i(q_1, q_2, q_3, q_4) = q_i$ ( $1\le i \le 4$ ), and

$P_{e x t} = \lim\limits_{t \rightarrow \infty} \operatorname{\rm{Prob}}\{\underline{\mathrm{X}}(t) = \overrightarrow{0}\} = q_{1}^{k_{1}} q_{2}^{k_{2}} q_{3}^{k_{3}} q_{4}^{k_{4}}$

where $k_i = X_i(0)$ for $1\le i \le 4$ .

Consequently, the probability of the establishment of infection is

$P_{out} = 1-P_{ext} = \begin{cases} 0, &\mbox{if $\mathcal{R}_0 \le 1$}, \\ 1-q_1^{k_1}q_2^{k_2}q_3^{k_3}q_4^{k_4}, & \mbox{if $\mathcal{R}_0 > 1$}. \end{cases}$

In Section 4, we estimate the clearance of HIV infection probability by computing $q_1, q_2, q_3$ and $q_4$ numerically. Our result shows that the linear branching process provides a very good approximation for the nonlinear CTMC near the IFE. In particular, the estimate for the clearance of HIV infection agrees closely with simulations of the CTMC model.

4. Numerical simulations

In this section, we numerically solve our CTMC model (defined in Table 4) with the based values of model parameter given in Table 1.

4.1. Viral dynamics without the presence of drug

First, we study the clearance of HIV infection probability in the absence of drugs (i.e., $\epsilon_{RT}^s = \epsilon_{RT}^r = \epsilon_{PI}^s = \epsilon_{PI}^r = 0$ ). Here the clearance of HIV infection is attained when $T_s = V_s = T_r = V_r = 0$ . In the simulation for the CTMC model, we estimate the probability of the clearance of HIV infection using the frequency of the clearance of infection from a total of $10,000$ sample paths, where each sample path is generated by using the Gillespie algorithm [55]. We also compare this estimate from the CTMC model to the corresponding analytical approximation computed from the multitype branching process theory. , displayed the clearance of HIV infection probability as a function of the viral clearance rate $c$ when the infection is initiated by a productively infected $T$ cell in the drug-sensitive and resistant strain (i.e., $T_s(0) = 1$ and $T_r(0) = 1$ ), respectively. This result shows that there is a strong agreement between the result of the CTMC model (obtained from Gillespie algorithm) and theoretical estimate (obtained from the multitype branching process approximation), and indeed the difference of these results is in the third decimal place. In , we plot the deterministic threshold parameter $\mathcal{R}_0$ and the probability of the establishment of infection $P_{out}$ as $c$ is varied, where the criterion for the establishment of infection is when the cumulative sum of $T_s, V_s, T_r$ and $V_r$ reaches $5,000$ . (resp. ) displays the associated result for $T_s(0) = 1$ (resp. $T_r(0) = 1$ ). One can see from this figure that $P_{out}$ is proportional to $\mathcal{R}_0$ . The smaller the $\mathcal{R}_0$ value, the lower the probability of the establishment of infection. More specifically, when $\mathcal{R}_0 \le 1$ , $P_{out} = 0$ , which implies that $P_{ext} = 1$ as we have seen in Eq (3.1). Furthermore, even in the severity of the disease (i.e., $\mathcal{R}_0 \gg 1$ ), $P_{out}$ appears to be still strictly below the unity. For instance, one can see from that, in our baseline case where $c = 23$ , $\mathcal{R}_0 = 3.13$ and $P_{out} = 0.69$ , and hence $P_{ext} = 1-P_{out} = 0.31$ ; namely, there is about 31% of the chance for the clearance of HIV infection when $\mathcal{R}_0 = 3.13$ . This shows that there is a positive chance for the clearance of HIV infection even if $\mathcal{R}_0 > 1$ . Additionally, shows that there is a higher chance for the establishment of infection with infection initiated with one infected $T$ cell in the sensitive strain as compared that in the resistant strain.

Figure 1. Probability of the clearance of HIV infection in the absence of drugs when the disease is initiated by (a) one

$T_s$ cell (b) one

$T_r$ cell. The green line corresponds to the analytic estimate from branching process approximation, and the red stars represent the result from the CTMC model. The estimates of the clearance probability are obtained based on 10,000 sample paths of the CTMC model. Here

$\epsilon_{RT}^s = \epsilon_{RT}^r = \epsilon_{PI}^s = \epsilon_{PI}^r = 0$ and the rest of parameter values are provided in Table 1.

DownLoad: Full-Size Img PowerPoint

Figure 2. Probability of a chronic infection in the absence of drugs when the infection is initiated by (a) one

$T_s$ cell (b) one

$T_r$ cell. The blue (resp., red) curve illustrates

$\mathcal{R}_0$ (resp.

$P_{out}$ ) as a function of

$c$ . The estimates of the clearance probability are obtained based on 10,000 sample paths of the CTMC model. The values of parameters are the same as in Figure 1.

DownLoad: Full-Size Img PowerPoint

Secondly, we study the time to the clearance of HIV infection and the establishment of infection using the CTMC model over $[0,365]$ days. In the numerical results shown, each histogram is generated from $10,000$ simulations of the CTMC model. displays the conditional probability distribution of time to the clearance of HIV infection, given that the clearance occurs within $[0,365]$ days. The left (resp. right) panel shows the result when the infection is initiated by $T_s(0) = 1$ (resp. $T_r(0) = 1$ ). The expected time to the clearance of HIV infection in the case of $T_s(0) = 1$ and $T_r(0) = 1$ are $0.62$ and $1.07$ days, respectively. shows the conditional probability distribution of the establishment of infection given that the establishment of infection occurs during $[0,365]$ days. We see that the establishment happens sooner than when the infection is initiated by $T_s(0) = 1$ as compared that by $T_r(0) = 1$ . The associated mean time to the establishment of infection are $2.08$ and $5.57$ days respectively. The corresponding mean and standard derivation are summarized in Table 6.

Figure 3. Conditional probability distribution of time to the clearance of HIV infection given that the eradication of infection occurs during [0,365] days. This histogram is generated based on

$10,000$ sample paths of the CTMC model. The parameter values are the same as in Figure 1.

DownLoad: Full-Size Img PowerPoint

Figure 4. Conditional probability distribution of time to a chronic infection given that the establishment of infection happens during [0,365] days. This histogram is generated based on

$10,000$ sample paths of the CTMC model. The parameter values are the same as in Figure 1.

DownLoad: Full-Size Img PowerPoint

Table 6. Statistics on time to the clearance of HIV infection and time to the establishment of infection given that the occurrence of an the clearance of HIV infection and the establishment of infection during

$[0,365]$ days. The numerical estimates of these statistics are based on

$10,000$ sample paths of the CTMC model. The parameter values are the same as in Figure 1.

Initiation	$P_{out}$	$E[T_{ext}\|ext]$	$\sigma[T_{ext}\|ext]$	$E[T_{out}\|out]$	$\sigma[T_{out}\|out]$
$T_s=1, V_s=T_r=V_r=0$	0.69	0.62	0.44	2.08	0.61
$T_r=1, T_s=V_s=V_r=0$	0.44	1.07	1.07	5.57	1.78

| Show Table

DownLoad: CSV

4.2. Dynamics of drug resistance under treatment

In this subsection, we investigate the case where the treatment is on (i.e., the corresponding efficacy parameters $\epsilon_{tr}^{st} > 0$ for $tr = RT$ or $PI$ , and $st = s$ or $r$ ). More specifically, for the example presented below, we assume that $\epsilon_{RT}^s = 0.5, \epsilon_{RT}^r = 0.02, \epsilon_{PI}^s = 0.1, \text{and }\epsilon_{PI}^r = 0.01$ . In this case the basic reproduction number is $\mathcal{R}_0 = 1.6873$ , with the reproduction number associated with the sensitive strain is $\mathcal{R}_s = 1.4087$ , and that of the resistant strain is $\mathcal{R}_r = 1.6873$ . This shows that the resistant strain is dominating under this specific the antiretroviral therapy (ART). Since $\mathcal{R}_0 > 1$ , the deterministic model (2.1) predicts that the ODE solution converges to a nontrivial endemic state for every biologically feasible initial condition that is other than the IFE (see the blue dashed curve in Figure 5(b) for an example). Prior to convergence, a high peak of the infection is visible in the ODE solution. However, this is not the case for the stochastic model. In the stochastic model, there is a positive probability for the clearance of the infection to take place. Two sample paths of the CTMC model are graphed in Figure 5, where one represents the clearance of the infection and the other represents a successful establishment of the infection.

Figure 5. Comparison the dynamics between the ODE model and the CTMC model in the presence of the ART. The blue dashed curve displays the ODE solution of

$V_r(t)$ and the black solid curve is a sample path of the associated CTMC model. (a) Stochastic solution indicates HIV infection can be cleared. (b) The displayed stochastic solution is consistent with the ODE solution showing the establishment of infection. In the upper right corner of (b), it zooms in the corresponding sample path of the CTMC model over

$[0, 8]$ days, where the

$y$ axis is displayed on a log scale. Here

$T_s(0) = T_r(0) = 1$ ,

$\epsilon_{RT}^s = 0.5$ ,

$\epsilon_{RT}^r = 0.02$ ,

$\epsilon_{PI}^s = 0.1$ ,

$\epsilon_{PI}^r = 0.01$ and the rest of the initial condition and the based values of the parameters are the same as in Figure 1. The ODE model is simulated by using the explicit Runge–Kutta method with the time step

$0.01$ .

DownLoad: Full-Size Img PowerPoint

Based on the CTMC model, for instance if the infection is initiated with a productively infected cell in both strains (i.e. $T_s(0) = T_r(0) = 1$ ), the probability of the clearance of HIV infection is $42.3 \%$ in the presence of ART, whereas the corresponding clearance probability is $18.4 \%$ in the absence of ART. This indicates that the treatment tends to increase the likelihood for clearance of the infection in patients.

Furthermore, we study the time to the clearance and the time to the establishment of infection by using the CTMC model. The simulated result is illustrated in . Particularly, the left panel displays the conditional probability distribution of the time to the clearance given the infection vanishes during a year. The right panel shows the distribution of the time to the establishment of infection provided the infection is successfully established during a year. Our result shows that the average time to the clearance of the infection is about $2.0$ days with a standard deviation of $1.9$ days. On the other hand, the mean time to the establishment of infection is about $6.3$ days with a standard deviation of $2.6$ days. Compared to the result in the absence of the ART (Section 4.1), we see that the time to the establishment of infection is expected to take longer in the presence of the ART.

Figure 6. Conditional probability distribution of time to (a) the clearance of HIV infection given the eradication of infection occurs during

$[0,365]$ days; (b) the establishment of infection given that the infection is successfully established during

$[0,365]$ days. The initial condition and the parameter values are the same as in Figure 5. This histogram is generated based on

$10,000$ sample paths of the CTMC model. The parameter values are the same as in Figure 5.

DownLoad: Full-Size Img PowerPoint

5. Discussion

In this paper, using a two-strain mathematical model we study the accumulation of mutations that may confer HIV-1 drug resistance before therapy (i.e. the drug efficacies are all zero) and how the mutant strain is selected and competes with the wild-type strain during treatment. We formulate a continuous-time Markov chain (CTMC) model based on the ODE models proposed and analyzed in [,]. We then apply a multitype branching process to approximate the CTMC at the IFE to obtain a theoretical estimate of the establishment of infection probability. Numerical simulations are carried out using the CTMC model. Our result shows that when the infection is initiated with one productively infected cell in the sensitive strain, the chance for the establishment of infection would be about $69\%$ . In contrast, when the infection is initiated with one productively infected cell in the resistance strain, the establishment of infection probability will be reduced to $44\%$ . On the other hand, in the case of infection initiated by the virus in either strain, the clearance of HIV infection happens almost immediately with probability $1$ . This is consistent in that the clearance rate is $23$ per day, which will always ensure that when the infection is initiated by virus with the density sufficiently low, clearance is done almost immediately. However, in the case of infected cells, cell death is $1$ per day, which will give a chance for infected cells to produce virions and infect other cells. The deterministic case shows that when $\mathcal{R}_0 > 1$ (deterministic threshold), there is always going to be disease persistence. However, our stochastic results show that this is not always the case. Under ART considered in this work, we see that the resistant strain dominates when $\mathcal{R}_s < \mathcal{R}_r$ . In this case, we see that even with $\mathcal{R}_0 > 1$ , stochastic dynamics indicate that there is a positive probability of clearance of the infection in patients. In addition, our stochastic investigations show that in the presence of antiretroviral therapy, the chance of the eradication of infection tends to increase although drug resistance is likely to emerge.

There are limitations in the current study. Our model predicts that either the infection is cleared or successfully established. However, current treatment cannot eradicate the virus and drug resistance may not be a major reason for HIV persistence during long-term effective therapy [56]. A major obstacle to viral elimination is the reservoir of latently infected cells. These cells are not affected by the treatment or immune response but they can be activated to produce new virions. The latent reservoir may also archive all the variants of HIV generated during treatment. Mathematical models have been developed to study the latent infection, low viral persistence and viral blips during antiretroviral therapy (e.g. see a review in [57]). The purpose of this paper is not studying HIV persistence under treatment. Here we used a simple stochastic model with one mutation to obtain some analytical results on the emergence of drug resistance and compare the results with the deterministic model. In addition, our model does not include the immune responses to HIV infection, which may play a critical role in explaining post-treatment control in some patients [27] and have important implications for the treatment with broadly neutralizing antibodies [58]. A more comprehensive model will be to consider cases such as cell-to-cell transmission, multiple mutations, incorporating the latently infected cells and immune responses, and studying these in the presence of antiretroviral therapy.

Data availability

The numerical codes used can be supplied by the authors upon request.

Acknowledgements

We thank the editor and reviewers for their helpful suggestions. L. Rong is supported by the NSF grant DMS-1950254.

Conflict of interests

The authors declare that they have no conflict of interests.

Appendix: Derivation of $\mathcal{R}_0$ of Model (2.1)

To calculate the reproduction number $\mathcal{R}_0$ for system (2.1), we apply the next generation matrix method [44].

It follows from direct calculation that the new infection matrix $F$ and transmission matrix $V$ are given by

$\begin{equation*} F = \begin{bmatrix} 0 &\frac{(1-u){(1-\epsilon_{RT}^s)}k_s \lambda}{d}&0&\frac{z{(1-\epsilon_{RT}^r)}k_r \lambda}{d} \\ 0&0&0&0\\ 0&\frac{u{(1-\epsilon_{RT}^s)}k_s \lambda}{d}&0&\frac{(1-z){(1-\epsilon_{RT}^r)}k_r \lambda}{d} \\ 0&0&0&0 \end{bmatrix} \quad V = \begin{bmatrix} \delta &0&0&0 \\ -{(1-\epsilon_{PI}^s)}N_s \delta&c&0&0\\ 0&0&\delta&0 \\ 0&0&-{(1-\epsilon_{PI}^r)}N_r \delta&c \end{bmatrix}. \end{equation*}$

By definition, $\mathcal{R}_0$ is the spectral radius of the matrix $FV^{-1}$ ; i.e.,

$\begin{equation*} \mathcal{R}_0 = \rho( FV^{-1}) = \begin{bmatrix} \frac{(1-u)N_s \lambda k_s {(1-\epsilon^s)}}{dc}& \frac{(1-u)\lambda k_s{(1-\epsilon_{RT}^s)}}{dc}&\frac{zNr \lambda k_r{(1-\epsilon^r)}}{dc}&\frac{z{(1-\epsilon_{RT}^r)}k_r \lambda}{dc}\\ 0&0&0&0\\ \frac{uN_s \lambda k_s{(1-\epsilon^s)}}{dc}&\frac{u\lambda k_s{(1-\epsilon_{RT}^s)}}{dc}& \frac{(1-z) N_r \lambda k_r{(1-\epsilon^r)}}{dc} &\frac{(1-z) \lambda k_r{(1-\epsilon_{RT}^r)}}{dc}\\ 0&0&0&0 \end{bmatrix}. \end{equation*}$

and hence we have

$\begin{array}{rcl} \mathcal{R}_0& = & \dfrac{\lambda}{2cd} \Big[ \big((1-u)k_s N_s{(1-\epsilon^s)} + (1-z)k_r N_r{(1-\epsilon^r)}\Big] \\ &&+\dfrac{\lambda}{2cd} \sqrt{ \bigg((u-1)k_s N_s{(1-\epsilon^s)} + (z-1)k_r N_r {(1-\epsilon^r)}\bigg)^2 - 4 k_r k_s N_r N_s {(1-\epsilon^s)(1-\epsilon^r)}} \end{array}$

References

[1]	WHO, Antibiotics resistance, 2020. Available from: https://www.who.int/news-room/fact-sheets/detail/antibiotic-resistance.
[2]	WHO, Antibiotics resistance, 2014. Available from: http://www.who.int/mediacentre/factsheets/fs194/en/.
[3]	Y. A. Luqmani, Mechanisms of drug resistance in cancer chemotherapy, Med. Princ. Pract., 14 (2005), 35–48. doi: 10.1159/000086183. doi: 10.1159/000086183
[4]	S. Falkow, Infectious multiple drug resistance, Pion Ltd, 1975.
[5]	F. M. Walsh, S. G. B. Amyes, Microbiology and drug resistance mechanisms of fully resistant pathogens, Curr. Opin. Microbiol., 7 (2004), 439–444. doi: 10.1016/j.mib.2004.08.007. doi: 10.1016/j.mib.2004.08.007
[6]	P. R. Hsueh, M. L. Chen, C. C. Sun, W. H. Chen, H. J. Pan, L. S. Yang, et al., Antimicrobial drug resistance in pathogens causing nosocomial infections at a university hospital in Taiwan, 1981–1999, Emerg. Infect. Dis., 8 (2002), 63–68. doi: 10.3201/eid0801.000454. doi: 10.3201/eid0801.000454
[7]	G. A. Curt, N. J. Clendeninn, B. A. Chabner, Drug resistance in cancer, Cancer. Treat. Rep., 68 (1984), 87–99.
[8]	G. Housman, S. Byler, S. Heerboth, K. Lapinska, M. Longacre, N. Snyder, et al., Drug resistance in cancer: an overview, Cancers, 6 (2014), 1769–1792. doi: 10.3390/cancers6031769. doi: 10.3390/cancers6031769
[9]	N. L. Komarova, D. Wodarz, Drug resistance in cancer: principles of emergence and prevention, P. Natl. Acad. Sci. USA., 102 (2005), 9714–9719. doi: 10.1073/pnas.0501870102. doi: 10.1073/pnas.0501870102
[10]	S. Blower, T. Porco, G. Darby, Predicting and preventing the emergence of antiviral drug resistance in HSV-2, Nat. Med., 4 (1998), 673–678. doi: 10.1038/nm0698-673. doi: 10.1038/nm0698-673
[11]	T. Cohen, C. Dye, C. Colijn, B. Williams, M. Murray, Mathematical models of the epidemiology and control of drug-resistant TB, Expert Rev. Respir. Med., 3 (2009), 67–79. doi: 10.1586/17476348.3.1.67. doi: 10.1586/17476348.3.1.67
[12]	AIDSinfo, Classes of drugs for HIV/AIDS, 2020. Available from: https://aidsinfo.nih.gov/drugs.
[13]	L. Rong, Z. Feng, A. S. Perelson, Emergence of HIV-1 drug resistance during antiretroviral treatment, Bull. Math. Biol., 69 (2007), 2027–2060. doi: 10.1007/s11538-007-9203-3. doi: 10.1007/s11538-007-9203-3
[14]	S. G. Deeks, Treatment of antiretroviral-drug-resistant HIV-1 infection, Lancet, 362 (2003), 2002–2011. doi: 10.1016/S0140-6736(03)15022-2. doi: 10.1016/S0140-6736(03)15022-2
[15]	J. M. Kitayimbwa, J. Y. Mugisha, R. A. Saenz, The role of backward mutations on the within-host dynamics of HIV-1, J. Math. Biol., 67 (2013), 1111–1139. doi: 10.1007/s00285-012-0581-2. doi: 10.1007/s00285-012-0581-2
[16]	A. A. Howard, J. H. Arnsten, Y. Lo, D. Vlahov, J. D. Rich, P. Schuman, et al., A prospective study of adherence and viral load in a large multi-center cohort of HIV-infected women, AIDS, 16 (2002), 2175–2182. doi: 10.1097/00002030-200211080-00010. doi: 10.1097/00002030-200211080-00010
[17]	B. M. Kane, HIV/AIDS treatment drugs, Infobase Publishing, 2008.
[18]	J. J. Eron, S. L. Benoit, J. Jemsek, R. D. MacArthur, J. Santana, J. B. Quinn, et al., Treatment with lamivudine, zidovudine, or both in HIV-positive patients with 200 to 500 CD $4^+$ cells per cubic millimeter, New. Engl. J. Med., 333 (1995), 1662–1669. doi: 10.1056/NEJM199512213332502. doi: 10.1056/NEJM199512213332502
[19]	A. R. McLean, M. A. Nowak, Competition between zidovudine-sensitive and resistant strain of HIV, AIDS, 6 (1992), 71–79. doi: 10.1097/00002030-199201000-00009. doi: 10.1097/00002030-199201000-00009
[20]	D. E. Kirschner, G. Webb, Understanding drug resistance for monotherapy treatment of HIV infection, Bull. Math. Biol., 59 (1997), 763–785. doi: 10.1007/BF02458429. doi: 10.1007/BF02458429
[21]	M. A. Nowak, S. Bonhoeffer, G. M. Shaw, R. M. May, Anti-viral drug treatment: dynamics of resistance in free virus and infected cell populations, J. Theor. Biol., 184 (1997), 203–217. doi: 10.1006/jtbi.1996.0307. doi: 10.1006/jtbi.1996.0307
[22]	R. M. Ribeiro, S. Bonhoeffer, M. A. Nowak, The frequency of resistant mutant virus before antiviral therapy, AIDS, 12 (1998), 461–465. doi: 10.1097/00002030-199805000-00006. doi: 10.1097/00002030-199805000-00006
[23]	T. B. Kepler, A. S. Perelson, Drug concentration heterogeneity facilitates the evolution of drug resistance, P. Natl. Acad. Sci. USA., 95 (1998), 11514–11519. doi: 10.1073/pnas.95.20.11514. doi: 10.1073/pnas.95.20.11514
[24]	J. M. Murray, A. S. Perelson, Human immunodeficiency virus: quasi-species and drug resistance, Multiscale Model. Sim., 3 (2005), 300–311. doi. 10.1137/040603024.
[25]	N. K. Vaidya, L. Rong, V. C. Marconi, D. R. Kuritzkes, S. G. Deeks, A. S. Perelson, Treatment-mediated alterations in HIV fitness preserve CD $4^+$ T cell counts but have minimal effects on viral load, PLoS Comput. Biol., 6 (2010). doi: 10.1371/journal.pcbi.1001012.
[26]	S. Moreno-Gamez, A. L. Hill, D. I. S. Rosenbloom, D. A. Petrov, M. A. Nowak, P. S. Pennings, Imperfect drug penetration leads to spatial monotherapy and rapid evolution of multidrug resistance, P. Natl. Acad. Sci. USA., 112 (2015), E2874–E2883. doi 10.1073/pnas.1424184112.
[27]	J. M. Conway, A. S. Perelson, Post-treatment control of HIV infection, P. Natl. Acad. Sci. USA., 112 (2015), 5467–5472. doi: 10.1073/pnas.1419162112. doi: 10.1073/pnas.1419162112
[28]	J. M. Conway, D. Coombs, A stochastic model of latently infected cell reactivation and viral blip generation in treated HIV patients, PLoS Comput. Biol., 7 (2011), e1002033. doi:10.1371/journal.pcbi.1002033. doi: 10.1371/journal.pcbi.1002033
[29]	J. M. Conway, B. P. Konrad, D. Coombs, Stochastic analysis of pre-and postexposure prophylaxis against HIV infection, SIAM J. Appl. Math., 73 (2013), 904–928. doi: 10.1137/120876800. doi: 10.1137/120876800
[30]	J. M. Conway, A. S. Perelson, Early HIV infection predictions: role of viral replication errors, SIAM J. Appl. Math., 78 (2018), 1863–1890. doi: 10.1137/17M1134019. doi: 10.1137/17M1134019
[31]	A. L. Hill, D. I. Rosenbloom, F. Fu, M. A. Nowak, R. F. Siliciano, Predicting the outcomes of treatment to eradicate the latent reservoir for HIV-1, P. Natl. Acad. Sci. USA., 111 (2014), 13475–13480. doi: 10.1073/pnas.1406663111. doi: 10.1073/pnas.1406663111
[32]	A. L. Hill, D. I. S. Rosenbloom, E. Goldstein, E. Hanhauser, D. R. Kuritzkes, R. F. Siliciano, et al., Real-time predictions of reservoir size and rebound time during antiretroviral therapy interruption trials for HIV, PLoS Pathog., 12 (2016), e1005535. doi: 10.1371/journal.ppat.1005535. doi: 10.1371/journal.ppat.1005535
[33]	H. C. Tuckwell, E. Le Corfec, A stochastic model for early HIV-1 population dynamics, J. Theore. Biol., 195 (1998), 451–463. doi: 10.1006/jtbi.1998.0806. doi: 10.1006/jtbi.1998.0806
[34]	D. Wick, S. G. Self, Early hiv infection in vivo: branching-process model for studying timing of immune responses and drug therapy, Math. Biosci., 165 (2000), 115–134. doi:10.1016/s0025-5564(00)00013-4. doi: 10.1016/s0025-5564(00)00013-4
[35]	A. Kamina, R. W. Makuch, H. Zhao, A stochastic modeling of early HIV-1 population dynamics, Math. Biosci., 170 (2001), 187–198. doi:10.1016/S0025-5564(00)00069-9. doi: 10.1016/S0025-5564(00)00069-9
[36]	T. Feng, Z. Qiu, X. Meng, L. Rong, Analysis of a stochastic HIV-1 infection model with degenerate diffusion, Appl. Math. Comput., 348 (2019), 437–455. doi:10.1016/j.amc.2018.12.007. doi: 10.1016/j.amc.2018.12.007
[37]	J. E. Pearson, P. Krapivsky, A. S. Perelson, Stochastic theory of early viral infection: continuous versus burst production of virions, PLoS Comput. Biol., 7 (2011), e1001058. doi: 10.1371/journal.pcbi.1001058. doi: 10.1371/journal.pcbi.1001058
[38]	D. B. Reeves, M. Rolland, B. L. Dearlove, Y. Li, M. L. Robb, J. T. Schiffer, et al., Timing HIV infection with a simple and accurate population viral dynamics model, J. R. Soc. Interface., 18 (2020), 20210314. doi: 10.1098/rsif.2021.0314. doi: 10.1098/rsif.2021.0314
[39]	M. Bofill, G. Janossy, C. A. Lee, D. Maconald-burns, A. N. Phillips, C.Sabin, et al., Laboratory control values for CD4 and CD8 T lymphocytes. implications for HIV-1 diagnosis, Clin. Exp. Immunol., 88 (1992), 243–252. doi: 10.1111/j.1365-2249.1992.tb03068.x. doi: 10.1111/j.1365-2249.1992.tb03068.x
[40]	H. Mohri, S. Bonhoeffer, S. Monard, A. S. Perelson, D. D. Ho, Rapid turnover of T lymphocytes in SIV-infected rhesus macaques, Science, 279 (1998), 1223–1227. doi: 10.1126/science.279.5354.1223. doi: 10.1126/science.279.5354.1223
[41]	A. S. Perelson, R. M. Ribeiro, Modeling the within-host dynamics of HIV infection, BMC. Biol., 11 (2013), 96. doi: 10.1186/1741-7007-11-96. doi: 10.1186/1741-7007-11-96
[42]	M. Markowitz, M. Louie, A. Hurley, E. Sun, M. D. Mascio, A. S. Perelson, et al., A novel antiviral intervention results in more accurate assessment of human immunodeficiency virus type 1 replication dynamics and T-cell decay in vivo, J. Virol., 77 (2003), 5037–5038. doi: 10.1128/jvi.77.8.5037-5038.2003. doi: 10.1128/jvi.77.8.5037-5038.2003
[43]	B. Ramratnam, S. Bonhoeffer, J. Binley, A. Hurley, L. Zhang, J. E. Mittler, et al., Rapid production and clearance of HIV-1 and hepatitis c virus assessed by large volume plasma apheresis, Lancet, 354 (1999), 1782–1785. doi: 10.1016/S0140-6736(99)02035-8. doi: 10.1016/S0140-6736(99)02035-8
[44]	Van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48. doi:10.1016/S0025-5564(02)00108-6. doi: 10.1016/S0025-5564(02)00108-6
[45]	N. Chitnis, J. M. Hyman, J. M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol., 70 (2008), 1272. doi: 10.1007/s11538-008-9299-0. doi: 10.1007/s11538-008-9299-0
[46]	Y. Yuan, L. J. Allen, Stochastic models for virus and immune system dynamics, Math. Biosci., 234 (2011), 84–94. doi: 10.1016/j.mbs.2011.08.007. doi: 10.1016/j.mbs.2011.08.007
[47]	L. J. Allen, G. E. Lahodny Jr, Extinction thresholds in deterministic and stochastic epidemic models, J. Biol. Dynam., 6 (2012), 590–611. doi:10.1080/17513758.2012.665502. doi: 10.1080/17513758.2012.665502
[48]	K. B. Athreya, P. Jagers, Classical and modern branching processes, Springer Science & Business Media, 2012.
[49]	K. S. Dorman, J. S. Sinsheimer, K. Lange, In the garden of branching processes, SIAM Rev., 46 (2004), 202–229. doi. 10.1137/S0036144502417843.
[50]	T. E. Harris, The theory of branching process, Berlin Springer, 1963.
[51]	P. Jagers, Branching processes with biological applications, Wiley, 1975.
[52]	H. W. Watson, F. Galton, On the probability of the extinction of families, J. Anthropological Inst. G. B. Irel., 4 (1875), 138–144.
[53]	A. Berman, R. J. Plemmons, Nonnegative matrices in the mathematical sciences, SIAM, 1994.
[54]	L. J. Allen, P. van den Driessche, Relations between deterministic and stochastic thresholds for disease extinction in continuous-and discrete-time infectious disease models, Math. Biosci., 243 (2013), 99–108. doi: 10.1016/j.mbs.2013.02.006. doi: 10.1016/j.mbs.2013.02.006
[55]	D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, J. Phys. Chem., 81 (1977), 2340–2361. doi: 10.1021/j100540a008. doi: 10.1021/j100540a008
[56]	D. B. Reeves, E. R. Duke, T. A. Wagner, S. E. Palmer, A. M. Spivak, J. T. Schiffer, A majority of HIV persistence during antiretroviral therapy is due to infected cell proliferation, Nat. Commun., 9 (2018), 1–16. doi: 10.1038/s41467-018-06843-5. doi: 10.1038/s41467-018-06843-5
[57]	L. Rong, A. S. Perelson, Modeling HIV persistence, the latent reservoir, and viral blips, J. Theor. Biol., 260 (2009), 308–331. doi: 10.1016/j.jtbi.2009.06.011. doi: 10.1016/j.jtbi.2009.06.011
[58]	D. B. Reeves, Y. Huang, E. R. Duke, B. T. Mayer, E. F. Cardozo-Ojeda, F. A. Boshier, et al., Mathematical modeling to reveal breakthrough mechanisms in the HIV antibody mediated prevention (AMP) trials, PLoS Comput. Biol., 16 (2020), e1007626. doi:10.1371/journal.pcbi.1007626. doi: 10.1371/journal.pcbi.1007626
[59]	L. J. Allen, Stochastic population and epidemic models, Mathematical Biosciences Lecture series, Stochastics in Biological Systems, 2015.

This article has been cited by:

Partha Sarathi Mandal, Sunil Maity, Impact of demographic variability on the disease dynamics for honeybee model, 2022, 32, 1054-1500, 083120, 10.1063/5.0096638

Reader Comments

Your name:*

Email:*
© 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Mathematical Biosciences and Engineering

3.9

Metrics

Article views(2989) PDF downloads(174) Cited by(1)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(6) / Tables(6)

Mathematical Biosciences and Engineering

Stochastic investigation of HIV infection and the emergence of drug resistance

Related Papers:

Abstract

1. Introduction

2. Deterministic model and analysis

2.1. Basic reproduction number

2.1.1. Case 1: Forward mutation ( $z = 0 \text{ and }0 < u < 1$ )

2.1.2. Case 2: Forward and backward mutation ( $0 < u, z < 1$ )

2.2. Sensitivity Analysis

3. Stochastic models

3.1. Continuous-time Markov Chain (CTMC)

3.2. Branching process approximation

4. Numerical simulations

4.1. Viral dynamics without the presence of drug

4.2. Dynamics of drug resistance under treatment

5. Discussion

Data availability

Acknowledgements

Conflict of interests

Appendix: Derivation of $\mathcal{R}_0$ of Model (2.1)

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

Mathematical Biosciences and Engineering

Stochastic investigation of HIV infection and the emergence of drug resistance

Related Papers:

Abstract

1. Introduction

2. Deterministic model and analysis

2.1. Basic reproduction number

2.1.1. Case 1: Forward mutation (z=0 and 0<u<1 z = 0 \text{ and }0 < u < 1 )

2.1.2. Case 2: Forward and backward mutation (0<u,z<1 0 < u, z < 1 )

2.2. Sensitivity Analysis

3. Stochastic models

3.1. Continuous-time Markov Chain (CTMC)

3.2. Branching process approximation

4. Numerical simulations

4.1. Viral dynamics without the presence of drug

4.2. Dynamics of drug resistance under treatment

5. Discussion

Data availability

Acknowledgements

Conflict of interests

Appendix: Derivation of R0 \mathcal{R}_0 of Model (2.1)

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog

2.1.1. Case 1: Forward mutation ( $z = 0 \text{ and }0 < u < 1$ )

2.1.2. Case 2: Forward and backward mutation ( $0 < u, z < 1$ )

Appendix: Derivation of $\mathcal{R}_0$ of Model (2.1)