Bacillus cereus strains from donor human milk and hospital environment: uncovering a putative common origin using comparative analysis of toxin and infra-red spectroscopy profiles

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 [19] 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 [20] 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:

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.

2.1.   Basic reproduction number

The dynamics of the system were studied by Kitayimbwa et al. in [15]. 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 [44], it follows from the analysis in [15] that $\mathcal{R}_0$ associated with system (2.1) is

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 [13]. 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

The second one is a boundary steady state

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

where

Here

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

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

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 [15] 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

with

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:

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 Tables 2 and 3, 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 Table 2 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, Table 3 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.

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)$,

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.

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

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

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 Table 5. Using this, we derive the offspring pgf in terms of $T_S, V_s, T_r, V_r$ and they are given by

By the theory of branching process [59], 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,

Then

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.

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.,

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

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

Consequently, the probability of the establishment of infection is

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. Figure 1(a), (b) 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 Figure 2, 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$. Figure 2(a) (resp. Figure 2(b)) 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 Figure 2(a) 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, Figure 2 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.

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. Figure 3 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. Figure 4 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.

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.

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 Figure 6. 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.

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 [13,15]. 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

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

and hence we have