The role of cell-to-cell transmission in HIV infection: insights from a mathematical modeling approach

1.   Introduction

As the causative agent of acquired immunodeficiency syndrome (AIDS), human immunodeficiency virus (HIV) infection remains a serious health problem in the world. At the end of 2021, there were approximately 38.4 million people living with HIV, about two-thirds of whom lived in Sub-Saharan Africa ^[1]. Since its discovery in 1981, the death toll of AIDS and the rapid spread of HIV have motivated a significant body of scientific and medical research. Over twenty antiretroviral drugs have been approved by the US Food and Drug Administration ^[2]. Currently, highly active antiretroviral therapy (HAART) containing three to five antiretroviral drugs has been successful in prolonging patients' survival and controlling viral replication. After 3–6 months of treatment, the viral load can be reduced below the limit of detection (i.e., 50 RNA copies per ml by standard assays), but supersensitive assays can reveal continual low-level viremia in patients under suppressive HAART ^[3,4].

The mechanisms underlying low-level viral persistence are not yet fully understood. Earlier work ^[5,6] suggested that viral persistence could be explained by the reactivation of cells from long-lived and stable latent reservoirs. Residual ongoing viral replication is also considered as a possible cause of viral persistence due to physiological sites or cellular compartments that have poor drug penetration ^[7,8]. Zhang and Perelson ^[9] showed that the kinetics of leakage of virus from follicular dendritic cells (FDC) is consistent with some clinical data, suggesting that FDC might also contribute to low-level viral persistence.

HIV can spread through the population of target cells by two basic routes ^[10,11]: (i) in cell-free virus infection, infected cells can produce free viruses, which infect susceptible target cells; (ii) in cell-to-cell spread, the virus can spread directly from an infected donor cell to a susceptible target cell, possibly through the formation of virological synapses ^[10,11,12]. There is experimental evidence that cell-to-cell spread is more effective than infection with free-floating viruses ^[13,14], and it could play an important role in HIV dynamics. During cell-to-cell viral transmission, the probability that at least one of the virions escapes from antiretroviral inhibition increases as the number of transmitted virions increases ^[15]. Therefore, we will explore if cell-to-cell viral spread can explain the low HIV persistence during long-term therapy in this paper.

Cell-to-cell transmission can lead to the transfer of tens to hundreds of virus particles in a single synapse, a certain fraction of which successfully integrates into the target cell genome ^[12,16], resulting in the multiple infection of cells. Multiply infected cells are frequently observed in experiments ^[17,18,19]. For example, Jung et al. ^[17] found that the HIV-infected splenocytes of two patients harbored between one and eight proviruses (a provirus is referred to as the HIV DNA integrated into the host cell's DNA), with an average of three to four proviruses per cell. During acute infection of macaques by simian immunodeficiency virus (SIV), infected CD4+ memory T cells carried an average of 1.5 viruses per cell ^[18]. Rebecca et al. ^[19] showed, using fluorescence in-situ hybridization, that an average of 3.4 proviral copies occurred per cell as a result of cell-to-cell viral transmission.

Mathematical models have been used to gain an understanding of HIV multiple infection ^{[20,21,22,23]}. Dixit and Perelson ^[20] first proposed probabilistic models of multiple infection and calculated the rate of infection when including only cell-free virus infection or both cell-free virus infection and cell-to-cell transmission. In a subsequent paper ^[21], they showed that the number of cells infected with $i$ viruses correlates with the $i$th power of the number of singly infected cells under certain conditions. Ito et al. ^[22] developed an ordinary differential equation (ODE) model with the heterogeneity of the target cell population, aimed at analyzing the HIV co-infection dynamics in the context of cell-free virus infection. They found that two different susceptible sub-populations were enough to reproduce the experimental data and around 98% of co-infected cells emerged from the more susceptible target cells. In a recent mathematical model of multiple infection ^[23], various numbers of virions were assumed to be spread with distinct probabilities.

Here we develop and examine a mathematical modeling framework incorporating multiple infected cell groups and two routes of viral spread. In order to investigate how multiple infection is formed, we compare the models with various scenarios of multiple infection to the HIV proviral copy data ^[17]. Simulating the model with the best-fit parameter estimates, we study the dynamics of infected cell populations, and evaluate the respective contributions of the two transmission routes to the singly and multiply infected cells. To determine if cell-to-cell viral spread can account for low HIV persistence in infected individuals receiving effective treatment, we also assess the sensitivity of the basic reproduction number and the level of HIV to the effectiveness of the drug.

2.   Models and mathematical results

2.1.   A model with two transmission routes

Based on basic models of within-host virus dynamics ^[24,25], we formulate the following model incorporating cell-to-cell viral spread alongside cell-free viral infection.

In the model, $T$ and $V$ denote the concentrations of uninfected CD4+ T cells and free virus, respectively. $T^*$ denotes the concentration of CD4+ T cells infected by free virions, and $I_i$ ($i = 1, 2, \dots, n$) represents cells infected through cell-to-cell spread with $i$ viruses transmitted. Uninfected cells are generated at rate $s$, and die at rate $d$. $\beta$ and $k$ represent the infection rates of uninfected cells via direct cell-free viral infection and cell-to-cell spread, respectively. The parameter $f_i$ ($i = 1, 2, \dots, n$) denotes the chance that cell-to-cell spread transmits $i$ viral genomes. Because the sum of all the possibilities is equal to 1, we have that $\sum _{i = 1}^{n} f_{i} = 1$. Infected cells die at the rate $\delta$ or $\delta_i$, corresponding to their respective class. It is assumed that virions are generated at a rate $N\delta$ from the cells infected by the cell-free virus, and at rate $N\delta_i$ from cells infected by cell-to-cell spread, where $N$ is the viral burst size, that is the total amount of virions created by one infected cell during its lifetime (i.e., $1/\delta$ for $T^*$ and $1/\delta_i$ for $I_i$). Free virions are assumed to be removed from circulation at the rate $c$. The parameter $\epsilon$ represents the drug's effectiveness in inhibiting cell-free virus infection, and $\epsilon_i \; (i = 1, 2, \dots, n)$ is the effectiveness of antiretroviral treatment in inhibiting the infection by $i$ virions via cell-to-cell transmission. Sigal et al. ^[15] suggested that when more virions are spread in each cell-to-cell transmission, the therapy is less likely to inhibit the infection induced by all the transmitted virions. Therefore, it is assumed that $\epsilon_1 > \epsilon_2 > \dots > \epsilon_n$. The schematic diagram of the model (2.1) is shown in Figure 1, and all parameters and their values are summarized in Table 1.

2.2.   Model analysis

We will use model (2.1) to analyze the HIV multiple infection distribution data in two patients who did not take any antiretroviral therapy ^[17]. Thus, we first let all the drug efficacies be 0, i.e., $\epsilon = \epsilon_i = 0$, for $i = 1, 2, \dots, n.$ The model with drug treatment will be discussed later. The model without antiretroviral therapy is given by

Any steady state of model (2.2) must satisfy the following equalities

It is clear that model (2.2) always admits an infection-free equilibrium $E_0 = (T_0, 0, 0, \cdots, 0)$ with $T_0 = s/d$, where there is no infection.

Using the next-generation method ^[30], we derive the basic reproduction number, given by (see Appendix A).

The basic reproduction number is a key measure in estimating the ability of a disease to spread in epidemiological research. It is defined as the average number of secondary transmissions from one infected person in a wholly susceptible environment. When it is greater than 1, the epidemic is growing; when it is less than 1, the disease is predicted to die out. In the context of the within-host viral dynamic model (2.2), $R_0$ represents the average number of newly infected cells (or virions) generated by one infected cell (or one virion) introduced into an environment where all cells are susceptible. The complexity of the expression, in particular the square root, makes it hard to explain the biological meaning. However, with straightforward algebraic computation, we find that the following $R'_0$ is an equivalent measure but it is easy to explain the biological meaning.

It can be easily shown that $R_0 > 1$, $R_0 = 1$, and $R_0 < 1$ are equivalent to $R'_0 > 1$, $R'_0 = 1$, and $R'_0 < 1$, respectively. $R'_0$ consists of three parts, each representing the contribution of infection by free virus, by cell-to-cell transmission from the infected cell populations $I_i$, and by cell-to-cell spread from the infected cell population $T^*$.

In addition to the infection-free equilibrium $E_{0}$, model (2.2) has an infected equilibrium $E_1 = (\hat{T}, \hat{T^*}, \hat{I_1}, \hat{I_2}, \dots, \hat{I_n}, \hat{V})$, where

and $\hat{T^*}$ is the unique positive solution of the following equation

with

To study the local stability of the infection-free equilibrium $E_{0}$ of model (2.2), we linearize model (2.2) and obtain the following Jacobian matrix evaluated at $E_{0}$

Similarly, we can obtain the corresponding Jacobian matrix of model (2.2) at the infected equilibrium $E_1$, given by

With distinct death rates of infected cells, the analysis of the model for local stability is challenging. However, with the assumption of equal death rates (i.e., $\delta = \delta_1 = \dots = \delta_n$), the infection-free equilibrium $E_0$ and new infection matrix $\mathcal F_{n+2}$ remain unchanged, and the basic reproduction number becomes $R_0 = \frac{N\beta}{c} T_0+\frac{k}{\delta}T_0$ (Appendix A). In this case, when $R_0 > 1$, there exists an infected equilibrium $E_1 = (\hat{T}, \hat{T^*}, \hat{I_1}, \hat{I_2}, \dots, \hat{I_n}, \hat{V})$, where

It is clear that the infected equilibrium exists if and only if $R_{0} > 1$. For the local stability of the infection-free equilibrium $E_0$ and the infected equilibrium $E_1$, we have the following results, with the proof given in Appendixes B and C, respectively.

Theorem 2.1. The infection-free equilibrium $E_0$ is locally asymptotically stable if $R_{0} < 1$ and unstable if $R_{0} > 1$.

Theorem 2.2. The infected equilibrium $E_1$ is locally asymptotically stable whenever it exists, i.e., when $R_{0} > 1$.

3.   Data analysis and computer simulations

3.1.   Multiple infection data analysis

Model (2.2) will be used to analyze the multiple HIV infection data in Jung et al. ^[17]. The data showed a bimodal frequency distribution for the number of HIV-1 proviral copies per infected cell in two patients who were not on antiretroviral therapy. The proviral copy number was between one and eight per cell, with a mean of about 3.2. Thus, we have $n = 8$ in our model. We calculate the fraction of cells containing $i$ proviruses in the overall infection when the system is at equilibrium (because the two patients were in the chronic stage of infection), and the proportion is denoted by $\psi_{i} \; (i = 1, 2, \cdots, n)$. We also reorganize the data by dividing the number of cells containing a given number of proviruses by the total number of cells collected from each patient ^[17]. In this way, we obtain the frequency of multiple infection and compare it with our model prediction.

Because of the limited frequency data in these two patients, we keep the assumption that all the infected cells have the same death rate. We also set some parameter values based on existing research ^{[23,26,27,28,29,31]}. In the absence of infection, CD4+ T cell levels are approximately $10^{6} \; ml^{-1}$ ^[31]. The death rate of uninfected cells $d$ is 0.01 $day^{-1}$ ^[23,26]. Thus, the generation rate of uninfected cells $s$ is $10^6\times0.01 = 10^4 \; ml^{-1}day^{-1}$. The death rate of infected cells $\delta$ is fixed to 1 $day^{-1}$, which is the same as in ^[23,27]. The clearance rate of virus $c$ is chosen to be 23 $day^{-1}$ ^[28,29]. We fit the calculated frequency $\psi_{i}$ (which will be derived later for each sub-model) to the proviral distribution data of each patient using the nonlinear least square method and provide estimates for the other parameters.

Multi-infection may be formed by cell-to-cell transmission originating from the directly virus-infected cells $T^{*}$, cell-to-cell transmission from the cells $I_i$, or a combination of both. In order to discriminate the potential mechanism of multiple infection formation and determine the respective contribution from each transmission route, we consider the following two sub-models and the full model.

Sub-model A: Only virus-directly-infected cells transmit viruses

If only the virus-directly-infected cells $T^{*}$ transit viruses, then the sub-model is given by

This simplified sub-model A has an infected equilibrium $E_A = (T_{A}, T^*_{A}, I_{1A}, I_{2A}, \dots, I_{nA}, V_{A})$, where

According to this equilibrium, we compute the fraction of infected cells containing $i$ proviruses, given by

The chance of spreading $i$ virions during cell-to-cell transmission ($f_{i}$, $i = 1, 2, \cdots, 8$) remains unknown. Inspired by the work of Dixit and Perelson ^[20], we first assume that $f_{i}$ takes the Poisson distribution $f_{i} = \frac{e^{-\theta}\theta^{i}}{i!}$, where $\theta$ is the mean value of viral genomes spread in each cell-to-cell transmission. Thus, $\sum _{i = 0}^{\infty}\frac{e^{-\theta}\theta^{i}}{i!} = 1$, which leads to $\sum _{i = 1}^{\infty}\frac{e^{-\theta}\theta^{i}}{i!} = 1-e^{-\theta}$, i.e., $\sum _{i = 1}^{\infty}\frac{e^{-\theta}}{1-e^{-\theta}}\frac{\theta^{i}}{i!} = 1$. We modify the Poisson distribution to $f_{i} = \frac{e^{-\theta}}{1-e^{-\theta}}\frac{\theta^{i}}{i!}$ to rule out the event that cell-to-cell transmission does not transmit any virus. The data (blue circles) and the best fit are displayed in Figure 2 upper panels. For both patients, the model fitting is similar in the distribution of proviral copy numbers. Only singly infected cells exist (Figure 2 upper panels). This is not consistent with the data that one cell has a mean of 3.2 proviruses ^[17]. Thus, model (3.1) with a modified Poisson distribution of $f_{i}$ fails to reproduce the multiple infection data ^[17].

We also use other distributions of the probability $f_{i}$ in the data fitting. In Figure 2 lower panels, we use a binomial distribution, which was also used in some other models ^[16,32,33]. We assume the chance of spreading $i$ genomes is

where $\gamma$ is the probability of having a successful infection from one virus. Using the binomial distribution, we have the same fitting (Figure 2 lower panels). Therefore, model (3.1) cannot explain the patient multiple infection data using either a Poisson distribution or binomial distribution. We also try some other distributions such as Gaussian and log-normal distributions and get similar fitting results.

Sub-model B: Only cell-to-cell spread transmits viruses

In this sub-model, we assume there is no direct cell-free virus infection and there is only cell-to-cell transmission. The sub-model is given by

The steady state of infected cells with $i$ proviruses is $I_{iB} = f_{i}(\frac{s}{\delta}-\frac{d}{k})$ and its frequency in the total infected cell population is given by

We compare the frequency with the same data set ^[17], as shown in Figure 3 upper panels. It was found that model (3.3) with the Poisson distribution of $f_{i}$ provides a reasonable data fitting, except for missing the first peak (i.e., infected cells with one proviral copy). Using a binomial distribution, the model also provides a reasonable fit to the data (see Figure 3 lower panels). Therefore, we conclude that model (3.3) assuming only cell-to-cell spread without direct cell-free virus infection provides a good fit to the patient data. Next, we investigate which distribution is better for reproducing the multiple infection data. For the comparison of the fitting using the two distributions, we calculate SSR (see below), the sum of squared residuals, which measures the error between model prediction and experimental data. We also compute AIC, the Akaike information criterion, which also evaluates the relative quality of different mathematical models in fitting a given set of data.

where $n$ is the number of data points from each patient, $\sigma$ is the number of parameters used in data fitting, $\psi_{i}$ is the predicted frequency of infected cells containing $i$ proviral copies calculated from our model and $\hat{\psi_{i}}$ is the corresponding patient data. Based on the best-fitting results, we present the SSR and AIC values in Table 2. It turns out that model (3.3) with the Poisson distribution of $f_{i}$ ($model_{im}$ of Table 2) has smaller SSR and AIC than the binomial distribution (i.e., $model_{ib}$ of Table 2). This indicates that the model with Poisson distribution is better at fitting the data.

Although the model prediction with only cell-to-cell spread (Eq (3.3)) agrees with the patient data, the fitting did not capture the bimodal proviral peaks at one and three proviral copies. We tried some other distributions, e.g., Gaussian distribution and log-normal distribution, but found that they cannot explain the two peaks either. The best fitting always shows a peak for infected cells containing three proviruses (similar to the fitting shown in Figure 3). Because the model assuming only virus-directly-infected cells can transmit viruses (i.e., Eq (3.1)) always exhibits a single peak for singly infected cells in the best fitting (Figure 2), we speculate that both transmission routes are required to capture the two-peak distribution of multi-infection. In the following, we use the full model to fit the same data set and further evaluate the relative contributions from each transmission mode.

Full model C: Both transmission routes can transmit viruses

To study if both the two transmission routes are required to capture the two-peak distribution of proviruses, we use the following full model

For the full model, the infected steady state $E_1$ is given in Eq\ (2.9). The frequency of infected cells containing $i$ proviruses is calculated from this steady state as follow:

Figure 4 shows the model fitting of the calculated frequency to the multiple infection distribution data. The upper panels use the modified Poisson distribution and the lower panels use the binomial distribution for the probability $f_i$. We find that both can capture the bimodal distribution at one and three proviral copies. Therefore, model (3.5) improves the fitting to the proviral distribution data. For a further comparison of the data fitting using the full model, SSR and AIC values were calculated for each distribution (see Table 2). The full model with the binomial distribution, denoted by $model_{bb}$, has smaller SSR and AIC values than $model_{bm}$, the model (3.5) with the Poisson distribution.

In conclusion, we find that the two routes of viral spread are both required to capture the two-peak data of multiple infection shown in experiment ^[17]. Model comparison from Table 2 suggests that the full model (3.5) with a binomial distribution for cell-to-cell transmission yields the smallest SSR and AIC values. The best-fitting parameter values using the full model and binomial distribution are provided in Table 3.

3.2.   Respective contributions to multiple infection

Data fitting to the multiple infection data in the last section suggests that the two transmission routes are both required to explain the two-peak proviral distribution. In this section, we evaluate the relative contribution of these two transmission routes to the multiple infection and total infection. Using the best-fitting parameter values of the full model, we simulate the dynamics of all of the infected cell populations. Figure 5(a) shows that all the infected cells containing a different number of proviruses have similar dynamics despite different magnitudes. For each class of $I_i$, the cell population increases rapidly, then decreases, oscillates, and eventually converges to a steady state. Using a different model, Dixit and Perelson ^[21] found that the abundance of the infected cell subpopulation at equilibrium scales with the number of proviruses present in cells, i.e., singly infected cells are most abundant, followed by cells containing two, three, or four viruses, etc. However, in our prediction, the cells containing three proviruses are most abundant, followed by cells infected with four, one, two, five, six proviruses, etc (see the zoom-in figure in Figure 5(a)). Cell-to-cell transmission that is included in our model may explain the difference.

To elaborate on the contribution to the total infection from single infection and multiple infection, we separate infected cells into two categories, single infection $T^*(t)+I_{1}(t)$ and multiple infection $\sum\limits_{i = 2}^{8}I_{i}(t)$. The infection can be established either from cell-free virus infection or cell-to-cell spread. The relative contributions from these two transmission routes are given by

Therefore, the contribution to the total infection from single infection and from cell-free virus infection, from single infection and from cell-to-cell spread, from multiple infection and from cell-free virus infection, and from multiple infection and from cell-to-cell spread is respectively given by

Using the same best-fitting parameter values as in Figure 5(a), we plot the above relative contributions in Figure 5(b). We observe that at the beginning of infection, cells infected with a single virus from cell-free virus infection are the most common (zoomed-in figure in Figure 5(b)). As the infection continues, the contributions to the total infection from multi-infection and from cell-to-cell spread increase rapidly. All curves have some minor fluctuation before converging to the respective steady states. Multiple infection from cell-to-cell transmission accounts for the majority of the total infected cells (about 70%). The contribution from multiple infection initially via cell-free virus infection accounts for approximately 10% of the total infection. The remaining infection is singly infected cells, in which cell-to-cell spread and direct viral infection account for 18 and 2%, respectively. This shows that cell-to-cell viral spread plays a crucial role in developing multi-infection and also accounts for most cell infection. The percentages may be slightly different due to the potential variation in parameter estimates. However, this should not affect the major conclusion in view of their relative magnitude.

3.3.   Influence of cell-to-cell spread on HIV persistence during therapy

To determine if cell-to-cell spread can explain the low HIV persistence despite long-term effective antiretroviral treatment, we study the viral change when drug efficacy increases. We start with the sensitivity of the basic reproduction number during treatment, denoted by $R_0^T$, with respect to the drug effectiveness against direct viral infection $\epsilon$ and the drug effectiveness against cell-to-cell spread $\epsilon_i \; (i = 1, 2, \dots, 8)$.

In the presence of therapy, the basic reproduction number given in (2.4) becomes

For each drug efficacy, we calcualte the normalized forward sensitivity index of $R_0^{T}$ as

where

Using the best-fitting parameter estimates in Table 3, we plot the sensitivity index in Figure 6(a). All the indices are less than zero, indicating a negative correlation between the basic reproduction number and drug efficacy. The basic reproduction number is relatively sensitive to the drug effectiveness against cell-free virus infection (i.e., $\epsilon$). However, the basic reproduction number is not sensitive to some drug effectiveness against cell-to-cell spread, such as $\epsilon_6, \epsilon_7,$ and $\epsilon_8$. This shows that cell-to-cell viral spread may explain the low viral load persistence during therapy, particularly when the drug effectiveness against cell-free virus infection is fixed.

We further study how the steady-state viral load changes as the drug efficacy increases. Like the sensitivity of the basic reproduction number, the viral load is relatively responsive to the alteration in the drug effectiveness against cell-free virus infection. However, it is not sensitive to the drug's effectiveness against cell-to-cell spread. (see Figure 6(b)). Taken together, these analyses show that the viral load is not sensitive to some drug efficacies against cell-to-cell transmission. This suggests that cell-to-cell viral spread can explain low HIV persistence despite long-term combination treatment.

4.   Conclusions and discussion

Increasing evidence from in vivo and in vitro studies suggests the high prevalence of HIV multi-infection ^{[17,18,19,34,35]}. However, the biological factors responsible for HIV multiple infection remain unclear, especially the two-peak distribution data in ^[17]. In this study, we construct a mathematical modeling framework including both infection from cell-free viral infection and cell-to-cell viral spread to explain HIV multiple infection. By fitting mathematical models to the distribution data of multiple infection, we discover that there are only single infected cells if the model only includes cell-to-cell spread from $T^*$ (Figure 2). Including cell-to-cell transmission alone without cell-free viral infection can generate the peak for cells containing three proviruses (see Figure 3). The results in Figure 4 show that both the infected cells $T^*$ and $I_{i}$ are required in our model to explain the two-peak proviral distribution. This suggests that the two modes of transmission both contribute to the formation of HIV multi-infection.

We did not include sequential cell-free virus infection in the model because the cell can down-regulate its surface CD4 molecules following the initial infection, rendering further infections difficult ^[36,37]. However, the CD4 down-modulation takes about one day, during which infected cells still remain susceptible to further infections ^[38]. In this case, a certain amount of co-infection takes place and may lead to multiple infection. The exact shape of the provirus distribution and which cell population is most abundant may depend on the choice of the rate of receptor down-modulation ^[39]. When the rate is relatively slow, two peaks occur in the distribution: one for cells infected with a single virus, and one for cells infected with four viruses (see Figure 7c in ^[39]). The distribution can be altered for different rates of receptor down-modulation, as demonstrated in ^[39]. The sequential virus infection with a decreasing infection rate due to CD4 down-regulation was included in a model and it was shown that the sequential virus infection alone cannot explain the multiple infection ^[40].

From some experiments and modeling work ^[15,33], cell-to-cell viral transmission, e.g., via the formation of virological synapses between cells, renders the infection less susceptible to antiretroviral drugs. In our work, we assume that the effect of treatment in inhibiting cell-to-cell viral spread decreases as the number of transmitted virions during cell-to-cell spread increases. Sensitivity analyses of the basic reproduction number and the viral load during therapy show that the model prediction is less sensitive to the drug effectiveness of inhibiting cell-to-cell spread, although it is still sensitive to the drug effectiveness of blocking cell-free viral infection (Figure 6). Therefore, cell-to-cell spread may play a critical role in explaining the low HIV persistence during suppressive therapy.

Some studies showed that the number of single infected cells is higher than the total number of multiply infected cells ^[21,41]. However, in our work and also in the experiment ^[17], the population of multiply infected cells is more frequent than the population of single infected cells. This difference can be attributed to the spread of the virus from one cell to another. In in vitro experiments and their modeling studies ^[21,41], the number of virions is much greater than the number of cells by several orders of magnitude. Therefore, cell-free viral infection dominates and cell-to-cell viral spread may be neglected. However, in in vivo experiments such as Jung et al. ^[17] and this study, viral transmission from cell to cell is highly possible. Model predictions based on the data on multiple cell infections in ^[17] show that more than 80% of the total infection is due to cell-to-cell virus transmission. This viral mode of transmission also leads to about 70% of multiple infections (Figure 5(b)), further suggesting that cell-to-cell viral spread plays a critical role in the development of HIV multi-infection in vivo.

Another reason we introduce cell-to-cell viral transmission is that infected CD4 T cells are isolated from the spleens of two individuals ^[17], where the cells are relatively densely packed and have a high probability of synapsing. In contrast, the cells in the blood mix better and may inhibit the formation of synaptic connections. This can lead to the frequent single infection of target cells ^[42]. The concentration of target cells is higher in the lymphatic tissues than in the blood, so cells may be infected several times either through sequential viral infection ^[39] or cell-to-cell spread discussed here. In summary, the routes of spread of the virus and thus the types of infected cells (single or multiply infected cells) may also vary depending on the tissue site.

HIV recombination is a crucial consequence of multiple infection. Recombination between genetically different viruses in the same infected cell has been observed in the spleens of HIV-infected patients ^[17]. Although the effect of recombination on HIV evolution is complex, the integration of different alleles into a single genome can confer drug resistance to antiretroviral drugs and facilitate evasion of the host immune responses ^{[43,44,45,46,47]}. The differential equation system employed here to examine the dynamics of multiple HIV infections may not be suitable for investigating HIV recombination, as this typically involves the use of bioinformatics tools and phylogenetic methods to analyze genetic sequence data and reconstruct the evolutionary relationships between different HIV lineages.

Elucidating the mechanisms underlying multiple infections and HIV persistence has important implications for understanding HIV pathogenesis and developing more effective treatment strategies to control infection. Our modeling study, combined with experimental data ^[17], shows that cell-to-cell virus transmission can play an important role in the development of multiple infection and HIV persistence despite long-term effective combination antiretroviral therapy. Therefore, all treatment strategies that can effectively block cell-to-cell transmission, e.g., by blocking the establishment of viral synapses that might be required for cell-to-cell transmission, would be crucial for the control or even elimination of the disease. This work also provides a novel mathematical modeling framework that can examine some other issues in HIV infection, such as HIV latent infection, which is considered another important obstacle to viral elimination in current antiretroviral treatment. HIV can establish latent infection very rapidly following initial infection, and the resulting reservoir of virus remains relatively stable during therapy but can be activated to release new viral particles. If this latent reservoir accounts for the majority of the persistent virus, it could explain the lack of viral evolution observed in patients during suppressive therapy ^[48,49].

Acknowledgments

T. Guo was supported by the Changzhou Scientific and Technological Program grant (CJ20220134), the Natural Science Foundation of Jiangsu Higher Education (22KJB110007), and the National Natural Science Foundation of China (12201077, 12071217). X. Wang was supported by the National Natural Science Foundation of China (No. 12171413), the Natural Science Foundation of Henan Province (222300420016), and the Program for Innovative Research Team (in Science and Technology) in Universities of Henan Province (21IRTSTHN014).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.   Derivation of the basic reproduction number

Using the next-generation method ^[30], we define the new infection matrix $\mathcal F_{n+2}$ and the transfer matrix $\mathcal V_{n+2}$, given by the following expressions. The infection matrix in an epidemiological model is a matrix that describes the rates at which individuals in different compartments of the model (e.g., susceptible, infected, recovered) become infected when they come into contact with individuals in other compartments. The transfer matrix, on the other hand, is a matrix that describes the rates at which individuals move between compartments in the model.

It is clear that $\mathcal F_{n+2}$ is non-negative and $\mathcal V_{n+2}$ is non-singular. Thus, the basic reproduction number can be calculated by the spectral radius of the next generation matrix $\mathcal M = \mathcal F_{n+2} \mathcal V_{n+2}^{-1}$, given by

Therefore, the basic reproduction number is

B.   Proof of Theorem 2.1

From the Jacobinan matrix (2.7), we obtain the characteristic equation for the infection-free steady state $E_{0}$ as follows:

where $\lambda$ is the eigenvalue and $I$ is the identity matrix with the same dimension as $J(E_0)$. Expanding (B.1), we obtain

Equation (B.2) always has a negative root $\lambda = -d$ and the negative root $\lambda = -\delta$ with multiplicity $n$. Thus, the local stability of $E_{0}$ is determined by the roots of following equation

Let

When $R_{0} > 1$, $F(\lambda) = 0$ has at least one positive real root. When $R_{0} < 1$, all roots of $F(\lambda) = 0$ have negative real parts. Thus, the infection-free equilibrium $E_{0}$ is locally asymptotically stable when $R_{0} < 1$ and unstable when $R_{0} > 1$. This completes the proof.

C.   Proof of Theorem 2.2

From the Jacobian matrix (2.8), we obtain the following characteristic equation at the infected steady state $E_{1}$

where $\lambda$ is the eigenvalue, and $\hat{T}, \hat{T^*}, \hat{I_i}$ and $\hat{V}$ are given in Eq (2.9). Expanding the above determinant and regrouping terms, we obtain the following characteristic equation

where

It's clear that Eq (C.2) always has a negative root $\lambda = -\delta$ with multiplicity $n$. When $R_{0} > 1$, the inequalities $a_{2} > 0$ and $a_{3} > 0$ hold. Moreover, we have

It follows the Routh-Hurwitz criterion that all the roots of (C.2) have negative real parts. Thus, the infected equilibrium $E_1$ is locally asymptotically stable when it exists, i.e., when $R_{0} > 1$. This completes the proof.