A hybrid stochastic-deterministic approach to explore multiple infection and evolution in HIV

To study viral evolutionary processes within patients, mathematical models have been instrumental. Yet, the need for stochastic simulations of minority mutant dynamics can pose computational challenges, especially in heterogeneous systems where very large and very small sub-populations coexist. Here, we describe a hybrid stochastic-deterministic algorithm to simulate mutant evolution in large viral populations, such as acute HIV-1 infection, and further include the multiple infection of cells. We demonstrate that the hybrid method can approximate the fully stochastic dynamics with sufficient accuracy at a fraction of the computational time, and quantify evolutionary end points that cannot be expressed by deterministic models, such as the mutant distribution or the probability of mutant existence at a given infected cell population size. We apply this method to study the role of multiple infection and intracellular interactions among different virus strains (such as complementation and interference) for mutant evolution. Multiple infection is predicted to increase the number of mutants at a given infected cell population size, due to a larger number of infection events. We further find that viral complementation can significantly enhance the spread of disadvantageous mutants, but only in select circumstances: it requires the occurrence of direct cell-to-cell transmission through virological synapses, as well as a substantial fitness disadvantage of the mutant, most likely corresponding to defective virus particles. This, however, likely has strong biological consequences because defective viruses can carry genetic diversity that can be incorporated into functional virus genomes via recombination. Through this mechanism, synaptic transmission in HIV might promote virus evolvability.


Introduction
Human Immunodeficiency Virus-1 (HIV-1) evolution is considered to be a driving factor of infection progression, which creates great potential for the virus to escape from immune system responses, form drug resistance, and overcome vaccine candidates [9,16,22]. HIV-1 has a relatively high mutation rate (3 × 10 −5 per base pair) [21], which along with a quick turnover of the virus population [32,10,26], is thought to contribute to virus evolution in vivo.
As mutation is an important component of HIV-1 infection and evolution, mathematical models are a useful tool in order to study how virus mutation and evolution occurs in vivo [23,34,27,7,1], especially where experimental tracking and manipulation of mutant virus strains can be difficult. Mutant strains are key for the ability of an infection to escape from immune responses as well as to the evolution of drug resistance during anti-viral therapy, as viral variants can escape from drug effects as well as immune cell clones.
Since mutation is rare and random, the dynamics of infection contain important stochastic effects. However, since it is often not computationally practical to stochastically simulate infections at realistic population sizes, many mathematical models of virus infection turn to ordinary differential equation (ODEs) models. These models describe the overall effects of how populations interact with each other, and thus can only describe the "average" generation and subsequent growth of mutant strains of the virus.
Multiple infection has been documented in HIV-1 infected patients [14], and multiple infection has been shown to be promoted by direct cell-to-cell transmission through virological synapses [19] [11,3]. Multiple infection occurs when a single cell is infected with multiple (possibly genetically distinct) strains of virus. When a cell is infected with multiple copies of virus, complex social interactions can take place, such as viral complementation, viral interference, recombination, and others [20,6,30,8,2]. Synaptic transmission occurs when copies are transferred through virological synapses between cells. This typically results in multiple copies of virus being tranferred during one infection event. This has a large effect on the multiple infection of cells, and additionally is likely to result in the direct cotransmission of different virus strains from one cell to another [5]. This is in contrast to free virus transmission, which generally results in the infection of cells with only one copy of virus.
In this paper, in order to study the stochastic effects of multiple infection and synaptic transmission, we implement a hybrid stochastic-deterministic algorithm, which was introduced in [29]. This algorithm was developed for population dynamic and evolutionary biology systems, and allows us to stochastically simulate the trajectories of many mutant strains of the virus. We use an ODE mathematical model in order to include an arbitrary number of mutant strains, and then describe the hybrid algorithm used to simulate these systems, which is both computationally reasonable and includes important stochastic effects of mutant populations at low numbers. This hybrid algorithm is then used to analyze the effect of multiple infection, by modulating the maximum infection multiplicity from a minimum of one (single infection) to higher integers (multiple infection). We also study the effect of synaptic transmission on the existence and spread of different mutant populations throughout infection, by varying the percent of different virus transmission strategies (free virus and cell-to-cell). We also use the hybrid algorithm to compare differences between the completely deterministic mathematical model and the stochastic simulations of it, to demonstrate the impact of stochastic effects on basic viral infection in vivo.

Mathematical model description
In this section we derive a mathematical model in order to simulate infection with a large number of mutant strains.

Single virus strain model
We begin with a model for HIV-1 infection in vivo, which includes both free virus transmission and synaptic cell-to-cell transmission [17]. We assume that cells are sufficiently well-mixed, such that an ordinary differential equation (ODE) model is enough to capture the relevant dynamics. To include the possibility of multiple infection, we let x i (t) represent the number of cells infected with i copies of the virus at time t, where i takes on whole numbers between 0 (uninfected) and N (maximum infection multiplicity). Descriptions of the model parameters can be found in Table 1.
With only one viral strain, the model iṡ We assume that N is large enough to not result in a large amount of cells near the end of the cascade. In the above equations, cell free virus transmission happens at rate β, and synaptic transmission at rate γ. Mathematically the two processes are similar, except the number of viruses passed during an infection event is one for free virus transmission and S for synaptic transmission. Therefore, the total number of infected cells is independent of transmission mode when the total rate of transmission, β + γ = c, is fixed. However, varying rates of free virus and synaptic transmission such that β + γ = c will have a large effect on the multiplicity of infection (average number of virus copies in infected cells) and the relative fraction of different strains in infected cells.
Let us define the number of infected cells which is the sum of all of the infected subpopulations. This model, and all subsequent models, have a virus free steady state and an infection steady state. The virus free steady state is The infection steady state is The stability of these steady states depends on the basic reproduction ratio, R 0 = λ(β+γ) ad . If R 0 < 1 the virus free steady state is stable and if R 0 > 1 then the infection steady state is stable. This is detailed and further explored in [17].

Including mutant strains
Mutation is an important component of HIV-1 infection. Therefore, we generalize the preceding ODE model to include an arbitray number of mutant strains. We assume that forward and back point mutations can occur independently at k locations, all with rate µ. This creates 2 k viral strains. In the case k = 3, we have the wild type (abc), three strains with a single point mutation (Abc, aBc, and abC), three strains with a double point mutation (ABc, AbB, and aBC), and a triple mutant strain (ABC). In the general case of k mutations, we denote the wild type (abc. . . ) as strain S 1 and the k-mutant (ABC. . . ) as the S 2 k strain. Single mutant through k − 1 mutant strains are denoted by some S i for i = 2, 3, . . . , 2 k − 1. Similarly as in previous models, upon each infection event, mutation can happen at each location. Therefore, any strain can mutate into any other strain, and the probability of such an event is based on the number of point mutations in sum of all infected populations at time t, . . , k} as the number of point mutations where strain S i and S j differ. The probability that strain S i mutates into strain S j (or vice versa) is then ..,s 2 k (t) be the number of cells that are infected with the respective number of viruses at time t. Let Z 1 (t), Z 2 (t), . . . , Z 2 k (t) be the sum of the fraction of all subpopulations infected with the relevant strain at time t, that is where we define |i| = s 1 +s 2 +. . .+s 2 k . To include mutation, when calculating each Z i we include the relevant mutation factors µ D(S i ,S j ) , which describe the likelihood of each strain mutating into each other. Examples of this are included in the Supplementary Information Section 1.1 (single mutants) and Section 1.2 (double mutants). Again we let Z(t) be the sum of all infected subpopulations. We have that We have that S ≥ 1 is the constant number of viral copies transferred per synaptic infection event. Therefore, for each synaptic infection event, we must choose S viral copies from the infected cell to pass on to the target cell. Assume that the infecting cell is infected with s 1 , s 2 , . . . , s 2 k copies of the respective number of viruses, and that we transmitŝ 1 ,ŝ 2 , . . . ,ŝ 2 k copies of the respective strains. Synaptic transmission can result in infection by any combination of these strains, as long aŝ s 1 +ŝ 2 + . . . +ŝ 2 k = S is satisfied. The likelihood for each of these synaptic infection events is based on the probability to pick each strain for infection (Z i for i = 1, 2, . . . 2 k ). For each transmission event, we multiply the probability to transmit the respective number of each strain Zŝ 1 1 Zŝ 2 2 . . . Zŝ 2 k 2 k by the number of ways to transmit that number of each strain, which is the multinomial coefficient Ŝ s 1 ,ŝ 2 ,...,ŝ 2 k = S! s 1 !ŝ 2 !...ŝ 2 k ! . Therefore, forŝ 1 +ŝ 2 + . . . +ŝ 2 k = S, we define the synaptic infection intensity as The ODEs are theṅ x 0,0,...,0 = λ − βZx 0,0,...,0 − γZx 0,0,...,0 − dx 0,0,...,0 , where any term with a negative index is 0 and cells can be infected with a maximum of N total copies of virus.

Including advantageous/disadvantageous mutants
There are many ways to model fitness differences between distinct viral strains. If all the strains are neutral, and a coinfected cell is chosen for an infection event, we assume that the probability to transmit a given virus strain is proportional to the fraction of the strain among all viruses in the cell. If different viral strains have different fitness, the difference is implemented during the infection event, which can correspond to different rates of reverse transcription. That is, an infecting strain is again chosen based on its fraction in the cell. A disadvantageous (advantageous) mutant would then have a lower (higher) rate of transmission upon infection. In this way, a strain's fitness is independent of whether or not it is contained within a coinfected cell. This method of modeling fitness is one choice among many, others are explored further in [17].
This can be implemented in a straightforward way, by modifying the previous models. Here, the mean fractions of each population are modified by a fitness parameter F i ∈ [0, ∞) for i = 1, 2, . . . , 2 k . In the non-neutral case (F i ≡ 1) we have that for instance Equation (11) becomes Therefore, higher values of F i are correlated with higher fitness for the i th strain. Also note that when non-neutral fitness landscapes are considered, Z(t) becomes the sum of all infected subpopulations weighted by their respective fitness and is no longer equal to just the sum of all infected subpopulations.

Model parameters
Based on experimental evidence, we estimate the kinetic parameters in the model. Recall from Section 2.1 that in the neutral case the basic reproduction ratio R 0 for these systems is R 0 = λ(β+γ) ad , the infection steady state for infected cells is λ a − d (β+γ) , and the virus free steady state for uninfected cells is λ d .
We have an infection steady state value of 3.1 × 10 7 infected cells [4]. Assuming that the death rate of infected cells is approximately 2.2 days [26], we assume that a = 0.45 days −1 . Further, we set R 0 = 8 [28,24]. Finally, we introduce parameter ν such that we have a total body population of 10 ν lymphocytes at virus free steady state. Then we have the system of equations These equations uniquely imply that λ ≈ 1.59 × 10 7 days −1 , d = λ × 10 −ν days −1 , and β + γ = 3.6 × 10 −ν days −1 . We assume S = 3 and to get an average multiplicity of infection of around 3-4 [15,12,13,31], we vary the rate of free virus transmission (β) and cell-to-cell transmission (γ). Widely varying estimates for average infection multiplicities have been published, and there is some uncertainty about that. Further, we choose ν = 9, and therefore have d ≈ 0.0159 days −1 . For initial conditions, we assume uninfected steady state and a single infected cell (infected with the wild type) [25], that is x 0,0,...,0 = λ d and x 1,0,...,0 = 1.
At these parameters, triple and further mutant strains are negligible at peak infection, that is on average less than half a cell is infected with these strains when the number of infected cells is maximized (around seven days into infection). After peak infection, each strain will tend toward selection-mutation balance. These patterns can be seen by running the model with triple mutation first to peak infection, and then many days past peak infection.

Hybrid algorithm
In this section we introduce the hybrid algorithm, which allows for stochastic simulations at large population sizes with many cell subpopulations.

Theoretical basis
Since our system includes mutations, which are rare, the dynamics will contain important stochastic effects. Because of the large population sizes and numerous cell subtypes, it is not computationally practical to simulate the entire system with a fully stochastic method such as the Gillespie algorithm. We therefore implement a hybrid stochastic-deterministic algorithm, introduced in [29]. This algorithm was developed specifically to describe biological population dynamics in evolutionary systems, and allows to simulate the trajectories of both large and small populations simultaneously.
The hybrid algorithm is based on the idea that if a cell population is sufficiently large, an ODE representation can provide a good approximation of most stochastic trajectories of the population. We can write the ODEs systems as a single vector equation dV/dt = F(V), where V is a vector that contains all the cell subpopulations. Let M > 0 be a given size threshold. We classify the cell population x i as small at time t if x i (t) < M, or large otherwise. We simulate the small populations stochastically using the Gillespie algorithm and use the ODEs for the large populations. Further details of the hybrid method are given in the Supplementary Information Section 2.1 and 2.2.

Implementation
The size threshold M is a very important parameter in the hybrid algorithm. If M = 0, then at each time point each non-zero population is classified as large and the hybrid algorithm is identical to the deterministic solution of the ODEs. If M is very large, that is larger than all populations at all times, then the hybrid algorithm is the same as the completely stochastic Gillespie simulation of the model and can be extremely computationally inefficient. For intermediate M > 0, the hybrid algorithm is computationally efficient and the averages over many hybrid simulations go from approximating the deterministic predictions to converging to the stochastic averages as M increases. Therefore, in order to efficiently approximate the completely stochastic implementation of the model, we need to choose a small intermediate M such that the results are close to the fully stochastic implementation.
We can achieve this by comparing the hybrid averages over many simulations to completely stochastic averages over many simulations for simplified models, such as assuming a constant large number of uninfected cells or using similar parameters that result in smaller and more computationally manageable population sizes. For these models, completely stochastic simulations can be carried out and allow us to determine what size threshold M is reasonable for the related models. Specifically, since the averages over many hybrid simulations start from the deterministic prediction (M = 0) and converge to the completely stochastic average, similarly to [29] we i) set some difference threshold ε > 0, ii) test multiple size thresholds M, and iii) choose the smallest M such that the hybrid average is within ε of the completely stochastic average. Table 2 contains approximate average computer simulation run times for the completely deterministic ODE system, the hybrid method, and the completely stochastic Gillespie algorithm. Each system is run for the single mutation, double mutation, and triple mutation models. All simulations include only free virus transmission with multiple infection (N = 3), represent established infections only (we ignore stochastic simulations in which the infection dies out), and are stopped once the infected cell population reaches 10 8 cells. The times for the ODE and hybrid simulations also depend on the ODE deterministic step size, h (here h = 10 −5 ). In general, the number of strains per model is 2 k and the number of equations (subpopulations) per model is Because the parameters chosen for the simulations in Table 2 Table 2: Approximate average run times for a single simulation for the completely deterministic ODE system, the hybrid method with different threshold values (M), and the completely stochastic Gillespie algorithm (rows) Each system is run for the single mutation, double mutation, and triple mutation models (columns). In each system we assume all strains are neutral (F i = 1 for all i). The other parameters are N = 3, h = 10 −5 , µ = 3 × 10 −5 , λ = 1.59 × 10 7 , β = 3.60 × 10 −9 , γ = 0, d = 0.016, and a = 0.45.

Choosing a size threshold M
In the case of virus dynamics, we have developed an analytical method for finding a lower bound on size threshold M, which is based on the notion of R 0 . This method does not depend on the number of mutations, infection multiplicity, fitness landscape, etc, and provides an analytic lower bound on M. Depending on the details of the model, it is possible that a higher M is needed to get accurate descriptions of mutant dynamics. In general, we can always confirm that a chosen M is large enough using the ε test described in the preceding section and in [29].
The basic reproduction ratio, R 0 , is the average number of newly infected cells per single infected cell at the beginning of the infection. Therefore, infections with larger R 0 will lead to quicker and more successful growth of the overall virus population. For R 0 > 1, the dynamics of the infected subpopulations at the beginning of the infection are driven by two effects: i) largely deterministic exponential growth, and ii) stochastic fluctuation. Stochastic fluctuation is most important at initial infection, when the number of infected cells is low. As the number of infected cells grows, the dynamics transition to more deterministic exponential growth [25]. Furthermore, as R 0 increases, the amount of fluctuating in the infected cell subpopulations decreases and deterministic growth dominates, which accounts for decreased stochasticity in the system and better prediction from the ODEs. Therefore, as R 0 increases, lower hybrid size threshold M is needed in order to approximate the fully stochastic dynamics.
In the deterministic system, if R 0 > 1 then the infection will never go extinct. This is not the case in the stochastic setting, where even if R 0 is quite large, a single infected cell can die out before successfully infecting other cells. The rate at which infections stochastically go extinct, 1 R 0 , can be approximated from the deterministic ODEs (see Supplementary Information Section 2.3 for details).
The rate of stochastic extinction of the hybrid method should match the rate of stochastic extinction of the simulations in the fully stochastic case. We define a successful infection as an infection that reaches at least K > 0 cells. For R 0 > 1, the hybrid size threshold M and size of successful infection K play similar roles in hybrid simulations. This is because once an infected cell subpopulation reaches M cells, the ODEs will take over and it will grow deterministically (i.e. the infection will never die out). This provides a nice lower bound for size threshold M for the hybrid method. Let δ > 0 be some small difference threshold. We define the lower bound size threshold, M, asM where . denotes the ceiling function. Details of this calculation are described in Supplementary Information Section 2.3.
Since increasing R 0 monotonically leads to a higher rate of successful infections, we see that increasing R 0 leads to decreasingM. This makes sense because as R 0 increases, the amount of fluctuating in the infected cell population decreases, resulting in better prediction from the deterministic ODEs. Therefore, for higher R 0 , lower M is needed in the hybrid method in order to reasonably approximate the fully stochastic trajectories.
While this analytic method only provides a lower bound for the size threshold parameter for the hybrid method to approximate the completely stochastic dynamics, it also transforms the more abstract choice of a size threshold M into the more concrete choice of a difference threshold δ which explicitly measures the difference between the rate of stochastic extinction and the potentially much lower rate of extinction for hybrid simulations. This is especially important for studying statistics of the stochastic runs (such as the variation in the numbers of mutants). (a) Deterministic prediction versus stochastic averages with standard error bars (too small to be seen). The insert represents the deterministic prediction versus stochastic averages with standard deviation bars. (b) Percent error between the deterministic predictions and stochastic averages shown in (a). We see that increasing R 0 generally leads to a better prediction, both in terms of the difference between the deterministic prediction and the stochastic average and also in the percent difference between the two. We also see that the variance of the stochastic simulations decreases with increasing R 0 . We have R 0 = λ(β+γ) ad , and the parameters are N = 3, µ = 3 × 10 −5 , λ = 1.59 × 10 7 , β = 3.60 × 10 −9 , γ = 0, and d = 0.016. The infected cell death rate a is adjusted to achieve the required R 0 . All strains are neutral.
4 Comparison of the ODE and hybrid models in the context of infection dynamics

The expected number of neutral mutants
We start by predicting the number of neutral mutant viruses. Since infections are always established in the deterministic setting for R 0 > 1, the ODEs do a better job of predicting stochastic averages over many runs if we restrict to only simulations in which a successful infection is established. If we calculate the averages of mutants over all stochastic runs, then the average for uninfected (infected) cells will be much larger (smaller) than the deterministic prediction. This is because when a stochastic extinction event occurs, the number of uninfected cells returns to the large λ d steady state and the number of infected cells is 0.
Even if we only take into account the stochastic runs, in which infections are successfully established and maintained for the relevant timespan, it is possible for the deterministic predictions and averages over many stochastic simulations to disagree significantly. As described in Section 3.3, for R 0 closer to 1, there is a higher degree of stochastic fluctuation, which is difficult for the ODE model to capture and predict. Accordingly, the variance of the stochastic simulations also decreases with increasing R 0 .
While HIV-1 infection is characterized by a relatively high value of R 0 ≈ 8, bacterial and other infections can have values of R 0 very close to 1. For these systems, it is even more important to consider how stochasticity affects the system dynamics. Figure 1 compares the deterministic prediction versus the stochastic average for the number of cells infected with a mutant strain when the infected cell population has reached 10 4 cells. The deterministic predictions are in blue and the stochastic averages are in yellow. We see that increasing R 0 leads to a better prediction, both in terms of the difference between the deterministic prediction and the stochastic average and also in the percent difference between the two. We also see that the variance of the stochastic simulations decreases with increasing R 0 . For larger infected cell populations closer to peak infection, the absolute difference between the ODE predictions and stochastic averages are much larger, but the percent error is similar. For instance an infected cell population of 10 7 gives similar results, with a relative error at R 0 = 1.1 of about 11%. Past peak infection, the wild type and mutant population will fluctuate around steady state values, which the ODEs correctly predict (and which can be determined analytically).
We now turn to studying the dynamics of single mutants that are advantageous or disadvantageous. These mutant strains are characterized by larger or smaller fitness values compared to the wild type viruses.

The expected number of advantageous/disadvantageous mutants
As in the neutral case, the ODEs do not give accurate predictions for the mean numbers of mutants, especially if R 0 is close to 1. Figure 2(a) compares the deterministic prediction versus the stochastic average (similarly to Figure 1) for the number of cells infected with an advantageous mutant strain when the infected cell population has reached 10 4 cells. Here, the fitness of the wild type is 1 and the fitness of the advantageous mutants is 1.1. Figure 2(b) compares the deterministic prediction versus the stochastic average for the number of cells infected with a disadvantageous mutant strain, where the fitness of the wild type is 1 and the the fitness of the disadvantageous mutants is 0.9.
In both cases (advantageous or disadvantageous) we again see that increasing R 0 leads to a better prediction, both in terms of the difference between the deterministic prediction and the stochastic average and also in the percent difference between the two. We also see that the variance and percent error between the ODEs and the stochastic simulations decreases with increasing R 0 . For advantageous mutants, the percent error for R 0 = 1.2 is 1116% and for R 0 = 8 is 51%. For disadvantageous mutants, the percent error for R 0 = 1.2 is 18% and for R 0 = 8 is 7%.
In the advantageous case, we see that as in the neutral case, the deterministic prediction is an overestimate of the number of cells infected with the mutant. This is because low numbers of cells infected with the mutant can experience stochastic fluctuation before taking off and spreading rapidly. However, for the disadvantageous mutant we see that the deterministic prediction is an underestimate of the number of cells infected with the disadvantageous mutant. This is because the disadvantageous mutant can achieve fluctuating spikes of larger numbers before being out-competed by the wild type.

Probability distributions of mutant numbers
The deterministic setting tells us nothing about the distribution of mutant numbers. Many biological contexts, such as the evolutionary fate of drug resistant or immune escape mutants, require information not only on the mean, but on the whole distribution. The hybrid method is especially useful here, because it allows for efficient computations at large population sizes. For relatively large R 0 values and large numbers of uninfected cells, the hybrid method produces very accurate results even for relatively small threshold values M (and thus requires little computational time). This is because with a large R 0 and large initial number of uninfected cells, the mutant dynamics depend largely on when the mutant is generated and grows to exceed a small threshold. Figure  3 shows an example of the distribution of the number of cells infected with the single and double mutant strain of the virus when the infected cell population reaches 10 6 cells for hybrid simulations with size threshold M = 50 and M = 500. We find that for both single and double mutants, the probability distributions with M = 50 and M = 500 are statistically the same, and thus that M = 50 is large enough to capture the stochastic dynamics.
Furthermore, the hybrid method also provides a good approximation of the stochastic distribution and variance for smaller R 0 > 1. Here again, the mutant dynamics depend largely on when the mutant is generated and grows to exceed a threshold. However, since there is more stochasticity, larger size threshold sizes M are needed. Figure 4 shows the distribution of the number of cells infected with the single mutant strain of the virus when the infected cell population reaches 10 6 cells for hybrid simulations with varying size thresholds. Here we see larger size thresholds M are needed for a small basic reproduction ratio R 0 . In summary, Figure 3 shows that for R 0 = 8 (less stochasticity) M = 50 is large enough, but Figure 4 shows that for R 0 = 1.5 (more stocahsticity) M = 50 is not large enough to capture the stocahstic dynamics, and larger size thresholds are needed.

The timing of mutant emergence
An important benchmark in the evolution of HIV infection is the emergence of mutant strains of the virus. Using the continuous deterministic ODE system, the best we can do to predict when    Figure 1 with R 0 = 8. the first single/double mutant virus will emerge is by finding where the number of single/double mutant viruses first reaches 1. However, with the hybrid model, which is stochastic and discrete for small populations, we can get a full distribution of when we expect the first single/double mutant virus to appear, both in terms of time and infected cell population size, which could be important for treatment. An example of such a distribution can be seen in Figure 5. Here, the deterministic model underestimates the mean time of single mutant generation by 13%, and underestimates the mean number of infected cells at first single mutant generation by 75%. Similarly, the deterministic model underestimates the mean time of double mutant generation by 14%, and underestimates the mean number of infected cells at first double mutant generation by 94%. In general, the ODE error will become larger as more mutations are introduced, because of the increasing importance of stochastic generation for more mutational steps.

Impact of multiple infection and synaptic transmission on mutant evolution
Mutant evolution is impacted by multiple infection [33]. Here, we use the hybrid algorithm to simulate stochastic infections, to determine to what extent multiple infection makes a difference under different evolutionary scenarios. To this end, we compare simulations where only single infection is taken into account (maximum infection multiplicity N = 1) and simulations containing multiple infection (maximum infection multiplicity N > 1). Multiple infection, and coinfection (when a cell is infected with distinct virus strains), is especially important because interactions between distinct mutant strains can occur within a cell (this includes the phenomenon of viral recombination, which goes beyond the scope of this paper [18]). In the context of multiple infection, synaptic transmission (where multiple copies of virus are transferred per transmission event) is also a factor. As stated before, the spread of mutant virus depends on two factors i) stochastic generation of the mutant strain, and ii) the subsequent spread of the mutant strain once it has been generated (which is more stochastic for small population sizes, and more deterministic for larger populations sizes).

Neutral mutants
We first study the dynamics of neutral mutants, and then consider more complex fitness landscapes.

Effect of multiple infection in the absence of synaptic transmission
First, we assume that only free virus transmission takes place. In this model, for R 0 = 8, most cells do not become infected with more than 11 copies of virus by the time of peak infection (about 6.2 × 10 8 infected cells). Therefore, for the multiple infection case we set maximum infection multiplicity N = 11. Figure 6 shows the distributions of cells infected with single and double mutants for simulations with single and multiple infection, close to peak infection. We record the number of mutants at 6×10 8 infected cells, as it is close to the peak and most stochastic simulations will reach this threshold. Here, we see that when multiple infection is allowed (N = 11), there is a larger number of both single and double mutant viruses at peak infection. The number of single mutants is higher because of two factors: (i) When multiple infection is included, there are more target cells the single mutant strains can infect. When limited to single infection, the wild type initially spreads, limiting the number of cells in which the mutants, which are generated later, can infect (as only uninfected cells can be infected). These effects of multiple infection contribute to increased mutant virus spread during early infection. In particular, near the peak infection, we observe that in the presence of multiple infection, the number of cells infected with a mutant virus is about 2.4 higher than that under single infection only; the standard deviation is also larger by a factor of about 4.2, figure 6(a). Both of these factors also apply to the double mutant strain, which also benefits from larger numbers of single mutants at peak infection, as higher numbers of single mutants will lead to higher numbers of double mutants through mutation. Under multiple infection, the mean number of cells infected with double-hit mutants is abut 2.0 times higher than that under single infection only, and the standard deviation is about 3.2 times higher, figure 6(b).
Finally, we tested the prevalence of triple mutants, by comparing simulations with N = 1 and N = 3 (not shown in the figure). While with or without multiple infection, we expect a significant number of single and double mutants at peak infection, the number of triple mutants is very low (less than 0.2 on average). The probability to have at least one cell infected by a triple-mutant is 3.3% under multiple infection, which is about 1.7 times higher than that under single infection (1.9%) (this result is significant with p = 2.5 × 10 −3 by the z-test, with 2.7 × 10 4 runs under single infection and 1.9 × 10 3 runs under multiple infection).
Multiple infection can also have complex effects that can hinder mutant virus spread. In the absence of mutation, if for instance an infecting cell is singly (or only) infected with the mutant, then it is guaranteed to pass on the mutant strain. However, if an infecting cells is infected with a single copy of the mutant but many copies of the wild type, then the mutant strain is unlikely to spread. In the case of only free virus transmission, this is not much of a factor as the multiplicity of infection at peak infection is relatively small (1 when N = 1 and about than 1.9 when N = 11).

Effect of multiple infection in the presence of synaptic transmission
Next, we study the role of synaptic transmission in the evolutionary dynamics of infection. Again, we compare simulations with single infection only (N = 1), and simulations that allow multiple infection, where we now include synaptic cell-to-cell transmission.
During cell-to-cell transmission, we assume that S ≥ 1 viruses are transferred to the target cell per successful infection event. We allow simultaneous transfer of all possible non-negative integer combinations of the different strains such that the total number of virus copies transferred is S. During an infection event, upon selection of each virus, the viral copy has a probability to mutate. This greatly increases the ability of synaptic transmission to generate mutant strains of virus. To see this, consider infection from a cell that is infected only by the wild type. During a free virus transmission event, there is only one chance for mutation into a mutant strain. However, during a synaptic transmission event, there are S chances for mutation into a mutant strain.
Including synaptic transmission results in increased multiplicity of infection. Therefore, we need to adjust the maximum infection multiplicity used in the model, to avoid spurious accumulation of cell numbers at the end of the infection cascade. For our model, for R 0 = 8, S = 3, and 100% synaptic transmission, it is rare for a cell to become infected with more than 25 copies of virus by the time of peak infection. Therefore, we set N = 26 for all simulations that include synaptic transmission.  Figure 1 with R 0 = 8.
We use varying rates of free virus transmission (β) and cell-to-cell transmission (γ) such that β + γ = c. In this way, we can compare the effect of synaptic transmission at different levels of contribution, while preserving R 0 , the steady state values, and the overall rate of transmission. We will refer to different values of γ (under a constant total transmission rate, c) as the different transmission strategies. Figure 7 shows the distributions for multiple transmission strategies for the single mutation model, close to peak infection. Here, we see that when synaptic transmission is included, the average number of mutants is maximized with only synaptic transmission. This is because synaptic transmission is better for generating mutants earlier and more often (here S = 3). As we decrease the amount of synaptic transmission, we greatly reduce the generation of mutant strains of virus, which results in a smaller number of cells infected with the mutant strain at peak infection. Not surprisingly, in the presence of synaptic transmission, the result of Section 5.1.1 remains valid: simulations with multiple infection correspond to larger numbers of mutants compared to simulations with single infection only.

Disadvantageous/advantageous mutants, complementation and interference
In this section, we examine the effect of multiple infection and synaptic transmission on mutant strains with varying fitness. We first consider the dynamics of non-neutral (disadvantageous and advantageous) single mutants. Then we consider scenarios where the mutant strain is disadvantageous when it is alone in a cell, but when a cell is coinfected with both the wild type strain and the mutant strain then both strains are neutral (complementation). Finally, we also consider the scenario of interference, when a single mutant (alone in a cell) is advantageous, but this advantage is removed once the cell is coinfected with a wild type virus.

Effect of multiple infection in the absence of synaptic transmission
We first assume only free virus transmission. Figure 8 shows the histograms of numbers of cells infected with disadvantageous mutants in the presence of complementation (a) and advantageous mutants in the presence of interference (b). Simulations with and without multiple infection were compared. As before, the number of cells was recorded near peak infection, and the maximum multiplicity parameter was for N = 1 for simulations with single infection only, and N = 11 for the model including multiple infection.
Again, we see that when multiple infection is included, there is a larger number of cells infected with the mutant virus at peak infection. We observe a 2.1-fold increase in the mean number of cells infected with the mutant, and a 3.3-fold increase in the standard deviation when multiple infection is included. In the case of disadvantageous mutations, when multiple infection is included, the effect of complementation comes into play by decreasing the disadvantage of the mutants (it does not in the single infection case as cells cannot be coinfected). Complementation helps disadvantageous mutant spread, as the disadvantageous mutant becomes neutral in coinfected cells.
In the case of advantageous mutants, we observe a 1.9-fold increase in the mean number of cells infected with the mutant, and a 1.5-fold increase in the standard deviation when multiple infection is included. The effect of interference here is the opposite, as it hinders advantageous mutant spread, because the advantageous mutant loses its fitness advantage in coinfected cells. However, the benefits of multiple infection to mutant spread outweighs this effect because the multiplicity of infection is still relatively small just before peak infection (about 1.9 copies of virus per infected cell on average), and therefore multiple infection still results in larger number of cells infected with the mutant virus.

Effect of multiple infection in the presence of synaptic transmission
As before, with multiple infection we can include synaptic transmission. Figures S2 and S4 show the number of cells infected with a disadvantageous/advantageous mutant near peak infection, and Figures S3 and S5 show the number of cells infected with a disadvantageous mutant with complementation/advantageous mutant with interference near peak infection. Furthermore, Figure 9 shows the mean number of mutants for the different transmission strategies and all fitness landscapes.
We first see that in all fitness landscapes considered, only synaptic transmission is the best for mutant evolution. This is again because synaptic transmission is superior at generating mutants, regardless of fitness landscape. We also see that at peak infection there are more advantageous mutants (without and with interference), followed by neutral mutants, followed by disadvantageous mutants (with and without complementation), because of fitness differences.
Disadvantageous mutants. For the disadvantageous mutant landscapes (Figures S2, S3, and 9(b-c)), we see that the mean number of cells infected with the mutant does not decrease monotonically with increasing contribution of free virus transmission. The reason is that only free virus is better at spreading the mutant, once it is generated, as free virus transmission can result in cells singly infected with the mutant, rather than coinfected with multiple copies of the wild type. This effect is most important for disadvantageous mutants, which have the most difficulty spreading after generation.
Since many cells infected with the mutant are coinfected with the wild type as well in the presence of synaptic transmission, in the scenario with complementation ( Figures S3 and 9(c)), The T test between the two distributions gives a p-value less than 10 −6 . Histograms represent 20 thousand hybrid simulations. Simulations in which infections are not established (or in the rare case a simulation does not reach the infected size threshold) are discarded. Simulations are stopped when the infected cell population is close to peak infection (6 × 10 8 cells). The other parameters are as in Figure 1 with R 0 = 8. we see that synaptic transmission becomes more powerful. This is because with complementation, the disadvantageous mutant gains fitness in coinfected cells, and larger contributions of synaptic transmission are correlated with higher number of coinfected cells. Therefore synaptic transmission is beneficial for complementation because it increases the multiplicity of infection and rate of coinfection (so mutant viruses are more likely to be with wild type virus in a cell). For only free virus transmission, complementation has a very small effect because not a lot of cells are coinfected with both the wild type and mutant.
Advantageous mutants. For the advantageous mutant landscapes (Figures S4, S5, and 9(d-e)), initial mutant generation is very important, since advantageous mutants will spread efficiently after generation. Synaptic transmission is better than free virus at generating mutants, so again we see that only (or mostly) synaptic transmission is optimal, and that the average number of cells infected with the mutant strain decrease monotonically with increasing contribution of free virus transmission.
However, since many cells infected with the mutant are coinfected with the wild type as well with high contributions of synaptic transmission, in the scenario with interference ( Figures S5 and  9(e)), we see that synaptic transmission hinders mutants spread with interference compared to the non-interference case. This is because with interference, the mutant virus loses fitness in cells coinfected with the wild type. For only free virus transmission, interference has a very small effect as there are few cells that are coinfected with both the wild type and mutant.

Discussion
In this study, we developed a hybrid stochastic-deterministic framework to simulate in vivo HIV dynamics in order to assess the effect of multiple infection and synaptic cell-to-cell transmission on virus evolution. We analyzed the effect of multiple infection and synpatic transmission for different fitness landscapes, including disadvantageous mutants with complementation and advantageous mutants with interference. Furthermore, we also describe and characterize the differences between the stochastic simulations and deterministic predictions, which is an important but understudied phenomena that has large implications in mathematical modeling of biological processes and virus dynamics, as many classical models of virus dynamics only make use of ODEs.
Stochastic effects of virus dynamics are especially important in the context of rare mutation, because small cell populations can experience large amounts of fluctuation. However, large cell numbers and multiple subpopulations make fully stochastic methods such as the Gillespie method difficult and computationally intractable. Therefore, we adapted a hybrid stochastic-deterministic algorithm [29] which treats small populations stochastically and large populations deterministically in order to feasibly simulate evolutionary infection dynamics of this complex system. This algorithm has large potential for applications to other viral dynamics models as well as other models in the fields of population dynamics and evolutionary biology, where large population sizes also make fully stochastic simulations or agent-based models difficult.
An important parameter in the algorithm is the size threshold M, which determines which populations are treated stochastically and and which ones can be considered large enough to be treated deterministically by the algorithm. Larger size thresholds capture more stochastic effects, but come at a large potential computational cost. A priori estimations of the size threshold are difficult, but here we introduce an analytical lower bound on the size threshold, using the fact that the rate of stochastic extinction using the hybrid algorithm should be the same as the rate of stochastic extinction in a completely stochastic setting. This is an important step toward being able to determine a reasonable size threshold a priori.
We started by using the hybrid algorithm to compare the population averages over many stochastic runs with deterministic predictions. Focusing on the prediction of the number of cells infected with mutant virus, we found that for neutral and advantageous mutants, ODEs overpredict the number of infected cells, and for disadvantageous mutants (with the disadvantage above a threshold) ODEs underpredict the number of infected cells. The relative error depends on system parameters but generally decreases with R 0 . Further, ODEs could be used to estimate the timing of mutant generation (by scaling the population variable to correspond to cell numbers and finding the time when the population of cells infected by the mutant virus reaches unity). We find that such ODE estimates underestimate the time of single-and double-mutant generation compared to the stochastic model.
Our study demonstrates that multiple infection is an important part of HIV infection. In the presence of multiple infection, higher number of mutants are expected. The reason is as follows. When multiple infection is included, there are more target cells in which to infect, which also results in faster generation of mutants relative to the total number of infected cells. We showed that at the peak infection, the number of cells infected with single-or double-mutant viruses is about twice as high in the simulations with multiple infections, compared to simulation with single infections only.
We also show the effect of synaptic transmission on the dynamics of an infection near peak infection. In general, the respective contributions of the two virus transmission pathways depend on the following trade-off: synaptic transmission is more effective at mutant generation, and free virus transmission and lower multiplicities are better at facilitating mutant spread, as the mutant viruses do not have to compete with the wild type strain (which is initially well-established) within a cell to spread. We investigated a spectrum of different transmission strategies ranging from purely synaptic to purely free-virus transmission, to estimate the number of cells infected with mutant viruses at infection peak. We found that in the case of neutral and advantageous viruses, purely synaptic transmission results in the largest number of cells infected with the mutant virus. This effect remains even in the case when the advantage of mutants is countered by interference (when coinfection with a wild type virus reverses the advantage conferred by the mutation, making the mutant virus neutral.) For disadvantageous mutants, we observe that both purely synaptic and purely free-virus transmission result in an increased number of cells infected by the mutant virus (the relative numbers depend on specific parameters). In the presence of complementation, when co-infection of a disadvantageous mutant with a wild type virus removes the disadvantage, we observe an increase in the mutant number for purely synaptic transmission.
While the hybrid method is capable of handling a large number of subpopulations, the number of "reactions" included in the stochastic part of the algorithm increases with (i) the number of different virus strains, (ii) the maximum multiplicity N , and (iii) the number of virus transferred per synapse S. If these parameters are too large, the number of reactions for the Gillespie algorithm can becomes too high to be computationally feasible (even if only small populations are handled stochastically). In general, the number of strains per model is 2 k and the number of differential equations (subpopulations) per model is N +2 k 2 k . The events in the model are uninfected cell generation (only one reaction), cell death ( N +2 k of reactions is on the order of 10 4 , each simulation becomes very computationally expensive, which happens, for example, if we consider triple mutants in the presence of synaptic transmission.
Another limitation of our study is that synaptic transmission is non-spatial, which makes it relatively more powerful for promoting mutant generation. In the spatial case, synaptic transmission is restricted to spread to neighboring cells only, and so it is less dominant [18]. In vivo data from humanized mice indicate that synaptic transmission can result in spatial clusters of infected cells [19]. One way to incorporate this in future work is to extend our current model to include a system of spatial subdivided regions of cells where synaptic cell-to-cell transmission can only occur between cells within the same region, but free virus transmission can occur between all regions. Each region would be simulated stochastically with our hybrid stochastic-deterministic framework.