Linking immunological and epidemiological dynamics of HIV: the case of super-infection

In this paper, a two-strain model that links immunological and epidemiological dynamics across scales is formulated. On the within-host scale, the two strains eliminate each other with the strain with the larger immunological reproduction persisting. However, on the population scale superinfection is possible, with the strain with larger immunological reproduction number super-infecting the strain with the smaller immunological reproduction number. The two models are linked through the age-since-infection structure of the epidemiological variables. In addition, the between-host transmission and the disease-induced death rate depend on the within-host viral load. The immunological reproduction numbers, the epidemiological reproduction numbers and invasion reproduction numbers are computed. Besides the disease-free equilibrium, there are two population-level strain one and strain two isolated equilibria, as well as a population-level coexistence equilibrium when both invasion reproduction numbers are greater than one. The single-strain population-level equilibria are locally asymptotically stable suggesting that in the absence of superinfection oscillations do not occur, a result contrasting previous studies of HIV age-since-infection structured models. Simulations suggest that the epidemiological reproduction number and HIV population prevalence are monotone functions of the within-host parameters with reciprocal trends. In particular, HIV medications that decrease within-host viral load also increase overall population prevalence. The effect of the immunological parameters on the population reproduction number and prevalence is more pronounced when the initial viral load is lower. AMS Subject Classification: 92D30, 92D40


Introduction
Mathematical models have a long history of enriching our knowledge of immunological and epidemiological processes within their separate scales. Viral and bacterial dynamics mathematical models have contributed immensely to our understanding of the within-host interplay of the pathogen with the host immune system in diseases such as HIV, Hepatitis C and malaria (see e.g. [14,21,30] and the references therein). Population-level dynamics of many diseases has been studied by many authors and numerous books and review articles have been devoted to the between-host interactions (see e.g. [7,18] and the references therein). Epidemiological models of infectious diseases have been used in the past hundred years to answer critical public health questions.
The importance of linking mathematical immunology and mathematical epidemiology cannot be underestimated particularly in the light of pathogen evolution and competition [3,17]. Since the seminal paper of Gilchrist and Sasaki [15] that nests a simple within-host model within a susceptible-infected (SI)-resistant epidemic model, a number of authors have considered the connection between these two scales, accounting for the epidemiological scale mostly though the epidemiological reproduction number [2,19,24].Another approach that commonly has been applied is to link the within-host dynamics with the between-host dynamics when both are described by ordinary differential equation models. Such an approach has been used in a number of articles, but it typically does not separate the different time scales on which immunological and epidemiological processes progress [9]. Of course, many articles have followed the nested approach of Gilchrist and Sasaki [4,22].
In this paper, we use the dynamical approach in [15] to formulate a multi-scale immunoepidemiological model of HIV [13]. A different multi-scale approach to immuno-epidemiological modelling of HIV has been proposed by Milner and Sega [25]. HIV's high virulence and mutation rates make that virus the preferred choice for the study of the evolution of pathogens across the scales [20]. Coombs et al. [8] were the first to use a unified approach to modelling competition of parasites at both the within-and between-host levels. In their model, however, infected individuals are continuously structured by the initial inoculum and an explicit trade-off mechanism for competition between the strains on population level is not modelled. Consequently, the authors find that in the absence of mutation on the within-host scale, the strain that maximizes the betweenhost reproduction number out-competes the other strain.
In the present article we investigate a scenario where the two strains are subjected to a complete competitive exclusion on the within-host scale. In the context of competitive exclusion, the strain with the maximal immunological reproduction number dominates. The two strains might not have necessarily emerged through mutation from one another. They differ in many of their characteristics and may belong to different subtypes of HIV [39,40]. On the between-host scale, the strain with the higher immunological reproduction number can super-infect individuals infected with the strain with the lower reproduction number. Because in our framework, the strain with the higher within-host reproduction number takes over almost instantaneously, individuals in the population are infected with only one of the strains. Super-infection is a coexistence mechanism considered previously predominantly on an epidemiological level in the context of evolution of virulence [5]. It is closely related to co-infection and some authors argue that super-infection follows from co-infection when the co-infected class equilibrates much faster [1,26]. In that sense, our results derived here for super-infection may also apply to co-infection.
We use the immuno-epidemiological model that we formulate to address two questions: First, an early epidemiological model of HIV structured by age-since-infection suggested that the variability of infectivity with age-since-infection which presumably follows the within-host viral load, may destabilize the epidemiological dynamics and lead to sustain oscillations [37]. Our model takes the dependence of the infectivity on the viral load explicitly. Consequently, we address the question: Do single-strain periodic solutions exist in the immuno-epidemiological model of HIV? As we show, for the simplest dependence of the transmission and disease-induced death rate on the dynamic within-host viral load, the single-strain equilibrium is locally asymptotically stable and Hopf bifurcation does not occur. We hypothesize that the single-strain model does not exhibit oscillations. We point out that the result that rules out oscillations concerns the boundary (semi-trivial) equilibria, and not the coexistence equilibria.
The second question that we address is: How do the within-host parameters affect the epidemiological reproduction number, prevalence of HIV, and invasion capabilities of the two strains?
We find that the epidemiological reproduction number is increasing with within-host parameters that increase the within-host reproduction number, and decreasing with those that decrease the within-host reproduction number. The within-host reproduction number does not depend on the infected cell death rate, and it is not immediately clear how the viral load depends on this parameter. We find that the epidemiological reproduction number is an increasing function of the infected cell death rate. The prevalence seems to depend in an exactly opposite manner on the within-host parameters. That suggests that HIV medications which decrease healthy cells infection rate β and/or the number of virions produced by an infectious cell ν in fact increase the population level of the prevalence of HIV -a trade-off also observed in practice. Furthermore, our results imply that the sensitivities of the epidemiological reproduction number and prevalence with respect to the immunological infection parameters are decreasing functions of the initial viral load, suggesting that by decreasing the amount of virus transmitted HIV medications also increase the sensitivity of the epidemiological parameters with respect to immunological parameters. Finally, we find that the epidemiological invasion reproduction number of strain one decreases with the increase in the immunological reproduction number of strain two, while the invasion reproduction number of strain two increases with increase in the immunological reproduction number of strain two. Although, which strain persists and which strain dies out on the population level is governed by the epidemiological invasion numbers, in the majority of the simulations this relationship depends strongly on the relationship between the within-host reproduction numbers. That is, we find out that the most common outcome is that the strain with the largest immunological reproduction number dominated on the population level.
We formulate our multi-scale model of HIV in the next section. We use a standard immunological and a standard SI epidemiological models. The two models are linked through age-sinceinfection and through the epidemiological parameters which depend on the within-host viral load. In Section 3 we discuss the trivial and semi-trivial equilibria and their stability. The local stability of the isolated-strain equilibria suggests that sustained oscillations in the single-strain dynamics do not occur. In Section 4 we establish the existence of an interior equilibrium under the condition that both epidemiological invasion numbers are larger than one. Section 5 is devoted to elucidating the connection between within-host viral dynamics and the epidemiology of HIV.

A two-strain immuno-epidemiological model with superinfection
In this section we formulate a two-strain immuno-epidemiological model of HIV with superinfection. Individuals infected with strain i experience within-host dynamics, given by a standard immunological single-strain model of HIV. Within a host, the time variable is denoted by τ and signifies time-since-infection. Upon infection, HIV attacks the CD4+ T lymphocytes, which we call the target cells and we denote by x when healthy and by y i when infected by the ith strain of the virus. The number of free virions of strain i is denoted by V i .
In this model, uninfected target cells are created at a constant rate r and are infected by interactions with free virions of strain i according to a mass action law with rate constant β i . Productively infected cells, y i , produce ν i > 1 new virions of strain i. The healthy target cell death rate is μ, and cells infected with strain i have death rate d i . The parameter δ i denotes the strain i viral clearance rate, and s i is the shedding rate of strain i. Shedding rate is given by the amount of virus that leaves an infected individual and may serve to infect another individual. The within-host model of strain i is given by the following system of ordinary differential equations: All parameters and dependent variables of this within-host model and their definitions can be found in Table 1. Variations of model (1) have been studied extensively and have been commonly used in the investigation of HIV [31,34] and Hepatitis C dynamics [27]. Mathematically, model (1) has been completely analysed [12].
The within-host reproduction number of strain i in model (1) is computed as follows: The reproduction number of strain i gives the numbers of secondary infections that one virus i-infected cell will produce in entirely healthy cell population during its lifespan as infected. Model (1) always has an infection-free equilibrium E 0 = (r/μ, 0, 0), and if R i > 1, then the infected equilibrium , We now introduce a two-strain epidemiological model with superinfection in the host. We denote by S(t) the number of susceptible human individuals, where t is the chronological time. We structure the infected individuals by age-since-infection τ . Let i j (τ , t) be the density of individuals infected by strain j. Individuals in the class i j (τ , t) experience the same within-host dynamics given by model (1) for strain j. We assume that the epidemiological dynamics of the HIV infection is driven by the following two-strain model with super-infection: In model (3) is the birth/recruitment rate and m 0 = m j (0) is the natural death rate of susceptible hosts at zero viral load. Here, we assume that all susceptible individuals have approximately the same equilibrium level of healthy CD4+ T cells. Furthermore, the infected hosts die at a variable rate dependent on their viral load m j (V j (τ )). We will assume the simplest form of m j (V j (τ )), that is, where m j V j (τ ) gives the additional host mortality due to the virus. It should be noted that there is an intrinsic issue with this assumption. While the increased viral load during the AIDS phase definitely impacts mortality, this is not the case for the early peak in the viral load. Immediately after infection the viral load grows exponentially, peaks at high levels and then drops low for the asymptomatic phase. This early peak is captured by our within-host model and would imply increased mortality. However, that is not what is observed in practice, where the acute phase causes only mild flu-like symptoms. One way to remedy this is to connect the disease-induced mortality to the T-cell count but this approach also has issues. Another is to assume that m j is, say, a step function, which is zero (or very small) during the acute phase, and a constant afterwards. We obtain our results for a constant value of m j , and note that these results may not extend if m j is in fact a step function. The transmission coefficient of HIV β j (V j (τ )) is also dependent on the within-host viral load. We may assume that β j (V j (τ )) is proportional to the viral load at a given age-since-infection τ . Hence, where V j (τ ) is the viral load of the strain j, q j represents the probability of successful transmission and s j is the shedding rate of strain j. We notice that in the most general case we assume that not only the two strains have different within-host dynamics but they also transmit differently, shed differently and cause different virulence. HIV is a highly mutable pathogen. It is very likely that there are significant differences between two strains which potentially may be from different groups and subtypes. The total population size is given by the sum of all classes: We assume that individuals are subjected to super-infection [29] in the sense that individuals infected with strain one can get into contact with individuals infected with strain two and get themselves infected by strain two. Vice versa, individuals infected with strain two can get into a contact with individuals infected with strain one and get themselves infected by strain one. The super-infection occurs at a reduced rate where the coefficient of reduction η i satisfies η i < 1. We note that individuals infected with strain i who get super-infected with strain j have a coefficient of reduction η i .
We assume that one of the strains 'takes over' almost instantaneously, but during the short span when the super-infected individuals are actually infected with both strains, the within-host dynamics of the two variants is governed by the following within-host model: The long-term outcome of the competition of the two strains in the above model is competitive exclusion. The reproduction numbers of the two strains are again given by Equation (2). The strain with the larger reproduction number wins and eliminates the strain with the smaller reproduction number. Model of the form (4) has first been used in [6,16]. Competitive exclusion has first been proven rigorously by De Leenheer and Pilyugin [11]. The parameters P j take on values zero and one and control the direction of the super-infection. The values P j depend on the size of the reproduction numbers R 1 and R 2 corresponding to virus strain one and strain two. They can be expressed as follows: We note that if R 2 > R 1 strain two super-infects individuals infected with strain one, that is P 1 = 1 and vice versa. The parameters P 1 and P 2 are not probabilities, but a binary tool to capture two symmetric models in one.
The epidemiological reproduction number of strain j in system (3) is given by the following expression: The reproduction number of strain j gives the number of secondary infections that one strain j-infected individual will produce in an entirely susceptible population during its lifespan as infectious. This is the reproduction number of the single-strain model and does not depend on the interaction between the strains, that is it includes no super-infection effects.
In the next section we compute explicit expressions for the equilibria and establish their local stability.

Disease-free and strain dominance equilibria and their local stability
System (3) has one disease-free equilibrium and two single-strain equilibria. The disease-free equilibrium, in which infection is not present, always exists and is given by In addition, if the basic reproduction number R j > 1, for each j there is a corresponding singlestrain equilibrium E j given by The total population size N * j at the strain j dominance equilibrium is We note that the system also may have an interior equilibrium point. We discuss this option in Section 4. Next, we turn to the linearized equations of the disease-free equilibrium. To introduce the linearization at the disease-free equilibrium E 0 , let We look for solutions of the form For each j we obtain the following eigenvalue problem: Note that S * 0 /N * 0 = 1. Solving the differential equation, we obtain and replacingz j in the equation forz j (0), we obtain the following characteristic equation for strain j: Now we are ready to establish the following result.
then the disease-free equilibrium is locally asymptotically stable. If max{R 1 , Consider λ with λ ≥ 0. For such λ and each j we have Hence, the equations G j (λ) = 1 do not have a solution with nonnegative real part and the diseasefree equilibrium is locally asymptotically stable. Now assume For that fixed j and λ real we have G j (0) = R j > 1. Furthermore, lim λ→∞ G j (λ) = 0. Hence, the equation G j (λ) = 1 has a real positive root. Therefore, the disease-free equilibrium is unstable.
To investigate the local stability of the endemic equilibrium E j for a fixed j, we recall that strain j dominance equilibrium E j exists if and only if the basic reproduction number corresponding to this strain R j is greater than 1. So we assume that R j > 1 for a fixed j. The local stability of the strain j dominance equilibrium depends on the invasion reproduction numbers R 2 1 and R 1 2 for each of the strains. By definition, the invasion reproduction number of the first strain R 2 1 gives the number of secondary infections that one infected individual with the first strain will produce in a population in which the second strain is at equilibrium. The invasion reproduction numbers of the first and the second strains in our case are given by where The results on local stability of single-strain equilibrium E j are summarized below Proposition 3.2 Assume for a fixed j, R j > 1. Then the single-strain equilibrium E j is locally Proof To simplify the presentation, we will assume without loss of generality that j = 1, that is R 1 > 1. We study the linearized equation around the single-strain equilibrium E 1 . We introduce the following notation for the perturbations An approach similar to [23] (see Appendix B in [23]) can show that the linear stability of the system is in fact determined by the eigenvalues of the linearized system (10). To investigate the point spectrum, we look for exponential solutions (see the case of the disease-free equilibrium) and obtain a linear eigenvalue problem. The equation for strain two separates from the whole system. We note that The system for the perturbations of strain two takes the form: Let From Equation (11) it follows that Substitutingz 2 (τ ) in the equation for Z 2 , dividing Z 2 from both side of the resulting equation, we get the following characteristic equation: We denote the right-hand side of that equation with L(λ), that is It can be seen that L(λ) is a decreasing function of λ for real λ. Furthermore, lim λ→∞ L(λ) = 0 and L(0) = R 1 2 . If R 1 2 > 1, then Equation (12) has at least one positive root, and therefore E 1 is unstable. If R 1 2 < 1, for λ with λ ≥ 0, we have Hence, the equation L(λ) = 1 has no solution with positive real part and all the eigenvalues of this equation have negative real part. Therefore, the stability of E 1 depends on the eigenvalues of the following system: First, we solve the differential equation Before we continue, we adopt the following notations From the first and third equation in Equation (13) we obtain Therefore, Linearizing the equation for the total population size It follows from the third equation of Equations (13) and (15) that Substituting Equations (15) and (16) in the first equation in Equation (13), and then cancelling z 1 (0) from both sides of the resulting equation, we obtain the following characteristic equation Here, we use the fact that That is, Evaluated for λ = 0 the above equation also gives Substituting above relations in Equation (17) we can rewrite the characteristic equation in the following form: For λ with λ ≥ 0 we have 1 On the other hand, R 1 > 1 implies m 1 /c 1 < 1. The quotients m i /c i are important quantities and will be used often. Hence, we adopt the following notation for them: We have that for λ with λ ≥ 0, the following inequality holds: Hence, for λ with λ ≥ 0, the characteristic equation H 1 (λ) = 1 has only solutions with negative real parts. Thus, the single-strain equilibrium E 1 is locally asymptotically stable. This concludes the proof.
In conclusion we notice that the results in this section imply that in the single-strain case the endemic equilibrium will be locally stable. Furthermore, the results can be extended to an n-strain scenario, although an appropriate assumptions for the super-infection among n strains have to be made.

Existence of the interior equilibrium
In the previous sections we discussed the existence of disease-free equilibrium and strain dominance equilibria and established the local stability of the equilibria under biologically interpretable conditions. In this section, we consider the existence of an interior equilibrium for which both strains are present. We show an interior equilibrium exists under the conditions that the two basic reproduction numbers R i and the two invasion reproduction numbers R j i corresponding to each strain are all greater than one. For simplicity, and without loss of generality, we assume that the two basic reproduction numbers in the immunological system satisfy R 1 > R 2 , in this case, P 1 = 0 and P 2 = 1. Assume E * = (S * , i * 1 (τ ), i * 2 (τ )) is the co-existence equilibrium of the system, then it satisfies the following equations: Setting Then system (19) becomes Solving the first differential equation in (20) and replacing i * 2 (τ ) in the resulting equation from (21), we obtain Substituting Equation (22) in the expression of Q 1 and canceling Q 1 , we obtain Substituting Equation (21) in the expression of Q 2 and then cancelling Q 2 , we have From the first equation in (20) and the fact that it follows that To derive the equation for N * in terms of Q 1 and Q 2 we integrate the equations for i * 1 (τ ) and i * 2 (τ ) in system (20) and add them. We obtain Hence, This gives the following expression of N * in terms of Q 1 and Q 2 : With expressions (25) and (27) for S * and N * in terms of Q 1 and Q 2 we obtain the following system in Q 1 and Q 2 : To show the existence of an interior equilibrium, we have to show that this system has a nonzero solution (Q 1 , Q 2 ). System (28) is nonlinear and cannot be explicitly solved. To show the existence of nonzero solution, we introduce the following notation: From Equations (24) and (25) (or the second equation in (28)) it follows that Q 2 can be expressed as a continuous function of Q 1 : We use this expression for Q 2 to eliminate Q 2 from the first equation in system (28). Consequently, we obtain an equation in Q 1 only: The existence of an interior equilibrium for Equation (19) is now reformulated as the existence of at least one positive Q 1 satisfying G(Q 1 ) = 1 for which Q 2 as computed from Equation (29) is nonzero. In what follows, we will prove this fact. Let Q 1 = 0, from Equation (29) we get the correspondingQ * 2 and it can be expressed as We note that in the expression above we have 1 − α 2 R 2 = m 0 ρ 2 > 0. Hence,Q * 2 is well defined and positive. It is easy to show that Let Q 2 = 0. From Equation (23) we obtain that the correspondingQ * 1 can be expressed as We note that in the expression above, we have 1 − α 1 R 1 = m 0 ρ 1 > 0. Hence,Q * 1 is well defined and positive. If Q 2 = 0, then from Equation (29) we have

Now we define a function
It is easy to see that F(0) = g(0) = R 2 > 1, F (Q 1 ) < 0, and also Therefore, there exists a uniqueQ * 1 such that Note that Hence, is a decreasing function, the above inequality implies thatQ * 1 <Q * 1 . Then, we have Combining Equations (31) and (32) we conclude that the equation G(Q 1 ) = 1 has a positive solution. It is easy to establish that for that value of Q 1 , the corresponding Q 2 is also positive. This establishes the existence of positive solution of system (28), and hence the existence of an interior equilibrium. We summarize the result in the following proposition.

How does within-host infection affect population level of disease?
Within-host variability of disease progression in HIV is significant. How does the within-host disease transmission affect the epidemiological prevalence of HIV? This is one of the most critical question that can be addressed with immuno-epidemiological models. The main approach to this question, however, is primarily through simulations as explicit solution of the immune HIV system is not feasible. In the next subsection, we estimate the parameters of the linked immuno-epidemiological model.

Parameter estimation
Two groups of parameters need to be estimated: the epidemiological parameters and the immunological parameters. We set the time unit of t and τ to be in days. A number of articles estimate parameters for the standard within-host HIV model [31,36], particularly viral clearance rate and infected T cells life-span. Considerable variability exits in estimates of the viral clearance rate. Perelson et al. [31] give the viral clearance rate δ ≈ 3 days −1 . We assume the baseline value for the viral clearance of strain one δ 1 . Faster viral clearance rates from the plasma have also been suggested, for instance δ 2 ≈ 23 days −1 [10]. We assume this faster clearance rate for strain two. These two very different viral clearance rates correspond to different viral production numbers.
Snedecor [35] uses 250 particles per infected cell, which we adopt for the first strain ν 1 = 250, but De Boer et al. [10] suggest much larger numbers. Given the much larger clearance rate, we take ν 2 = 1000.
Perelson et al. [31] also estimate the infected cell life-span given by 1/d. It gives a value of d 1 ≈ 0.5. This is the value we assume for strain one. Now it is widely assumed that the infected cell life-span is about 1 day. Consequently, we assume for strain two d 2 = 1. The uninfected CD4+ T-cell population in the human blood is given to be 5 * 10 9 in [35]. Stafford et al. [36] give the death rate of uninfected T cells at about μ = 0.01. Since the equilibrium of uninfected T cells in the model is given by r/μ, then r = 5 * 10 7 T cells per day.
The shedding rate is different for the acute and chronic stage of infection. Several references give information about the viral load in copies per millilitre in the vaginal tract of HIV-infected women and in the semen of infected men. Estimating the total amount of virus shed is possible for the males. Pilcher et al. [32] give the mean plasma viral load in the chronic stage as approximately 100,000 copies/mL. Average 190 lb male has approximately 6000 mL of blood so the total viral load is approximately 6 * 10 8 . Pilcher et al. [32] also give the viral load in the semen as approximately 8000 copies/mL. The mean amount of semen for an average male is 5 mL so that the total amount of virus in the semen is 48,000 copies RNA. That gives a shedding rate of about s 1 = 0.00008, which we adopt for strain one. For strain two, we adopt twice as big shedding rate s 2 = 0.00015. The amount of virus in the semen can help estimate the initial viral load V 0 . If all the amount of virus in the semen is absorbed and passed into the blood, then the initial V 0 ≈ 20,000 viral particles, considering that one viral particle has two RNA copies. We vary V 0 within the range (1000, 100,000). These parameters give an immune reproduction number of strain one of R 1 = 6.22 and of strain two of R 2 = 6.95. These reproduction numbers are similar to the median within-host reproduction number of 5.7 reported by Stafford et al. [36]. With these parameters, the equilibrium viral load in the blood is computed to be 3.6 * 10 9 which is consistent with what is reported in [35].
To estimate the epidemiological parameters, we take the probability of transmission q 1 = q 2 = 1, so that c 1 = 0.00008 and c 2 = 0.00015. We define m 1 and m 2 so that the resulting epidemiological reproduction number falls within the reasonable for HIV range. Nishiura [28] reports epidemiological reproduction numbers in the range 3.5 − 4 for west European countries. That gives an estimates for m 1 = 0.00002 and for m 2 = 0.00004.
To determine the population demographic parameters, we obtain the average human lifespan from [33] to be 70 years. Hence, m 0 = 1/(365 * 70). Population Reference Bureau also gives the total human population in 2011 as nearly 7 billion people. That gives ≈ 275,000. Parameter values are listed in Table 2.

How do the within-host transmission parameters impact the epidemiological reproduction number and HIV prevalence?
To answer this question we plot the epidemiological reproduction number as a function of the immunological and epidemiological parameters. In this subsection, we focus on strain one only. The dependence of the viral load V j (τ ) on the infection rate β j is not evident by inspection of the within-host system. Plotting the epidemiological reproduction number as a function of β j reveals that the epidemiological reproduction number is an increasing function of β j (Figure 1). Furthermore, Figure 1 reveals that even for very small β j , the epidemiological reproduction number is well over one. At first this observation seems counterintuitive. To understand this result, we notice that after integration by parts the epidemiological reproduction number can be written as  Assume now that β j = 0 and y j ≈ 0. Then, In Equation (33) under the integral we have an exponent of the following integral: even when V j (0) = 500. Raising the exponent to that number results in a small number. Consequently, the value of the integral in formula (33) is small. This means that the epidemiological reproduction number does not depend in a significant way on the within-host dynamics, and we have which is what we observe in Figure 1. Plotting the epidemiological reproduction number against each within-host parameter reveals that the reproduction number is a monotone function of the parameter when the remaining parameters are fixed as strain one in Table 1. The reproduction number is also an increasing function of r, ν and d. However, the reproduction number is a decreasing function of δ and μ. To produce Figure 1 (left), we have taken β as a parameter, solved the within-host single-strain system, obtaining V (τ ). Now V (τ ) also depends on the parameter β and from this, we have computed R 0 as a function of β. The simulation is executed with the computer algebra system Mathematica.
For the parameters of the within-host system, Figure 1 also shows that the reproduction number is bounded and tends to a finite limit as an immunological parameter increases. The value of the limit is not sensitive to the virus initial amount V 0 . We notice that the epidemiological reproduction number changes over a very small range as the within-host parameter varies. Maximum change of the epidemiological reproduction number occurs at a critical value of the parameter which makes the immunological reproduction number equal to one.
The epidemiological reproduction number is much more sensitive to the epidemiological parameters. It increases when the shedding rate increases and does so almost linearly, without an upper bound ( Figure 2). R 0 increases with s and decreases with m 0 and m 1 . Figure 2 shows that as a function of s, the epidemiological reproduction number does not depend on the virus initial amount V 0 .
The dependence of the HIV prevalence on each immunological parameter is also monotone when the remaining parameters are fixed at values listed in Table 1 as strain one. We define the total prevalence at an equilibrium as The dependence of the total prevalence on within-host parameters, however, seems reversed compared with the reproduction number. For instance, HIV total prevalence is decreasing with increasing β (Figure 3). This dependence implies the somewhat surprising result that increasing β, which increases the within-host viral load, decreases the equilibrial prevalence of the disease. In fact this dependence mimics correctly the current situation with HIV, where various drugs reduce the within-host viral load, including through reduction of β, and extend the lifespan of infected individuals. As a result, even though the incidence is declining, since infected individuals live longer, the overall prevalence is rising [38]. The immuno-epidemiological model in this article captures this trade-off between treatment and prevalence in HIV. We believe this results is predicated on our assumption that the disease-induced mortality is governed by the within-host viral load so that when the within-host viral load is reduced, so is the disease-induced mortality. We believe this effect will be retained even if the disease-induced mortality is dependent on the T cells rather than the viral load. Reducing β will reduce the viral load, keep the T-cell load larger and consequently the disease-induced mortality small.
Elasticity is a concept that measures the responsiveness of a quantity Q with respect to a parameter p. Elasticity is defined by the following formula: The elasticity ε Q p gives the percent change in the quantity Q given 1% change in p. The elasticity is positive if Q is increasing with p, and negative if Q is decreasing with p. We investigate the elasticity of the epidemiological reproduction number and the prevalence with respect to the within-host infection-related parameters: ν, β, d and δ. The elasticities of R 0 with respect to ν, β, d are positive, and the elasticity of R 0 with respect to δ is negative. The elasticities of the prevalence I 1 with respect to ν, β, d are negative, and the elasticity of I 1 with respect to δ is positive. To compare the magnitude of the elasticities with respect to the different within-host parameters, we plot the absolute values of the elasticities in Figure 4. Furthermore, we plot the magnitude of the elasticities with respect to the initial viral load V 0 . Because the initial viral load  potentially varies within a wide range, and may change as medication is applied, we would like to know how the magnitude of the sensitivities changes with respect to the average initial viral load. Figure 4 shows that the epidemiological reproduction number R 0 is less sensitive to the withinhost infection parameters than the prevalence. Furthermore, the magnitudes of the elasticities are generally decreasing functions of the initial viral load V 0 . In other words, as the average initial viral load in the population decreases, the epidemiological reproduction number and the prevalence become more sensitive to the within-host parameters. Put in the public health perspective, this suggests that HIV medications that decrease ν or β in fact act in two ways: by decreasing the individual viral load, and by increasing the sensitivity of key epidemiological quantities such as the prevalence and the incidence to the change of these parameters. In other words, decreasing the average initial viral load in the population leads to higher impact of these medications on these key epidemiological parameters.

How does within-host infection affect the invasion numbers of each strain?
The population-level prevalence of each strain is governed by the population invasion reproduction numbers given by expressions (9). In general, the numbers suggest the following outcomes of the population-level competition: (1) If R 2 1 < 1 and R 1 2 > 1, strain two excludes strain one; (2) if R 2 1 > 1 and R 1 2 < 1, strain one excludes strain two; (3) if R 2 1 > 1 and R 1 2 > 1, the two strains coexist; (4) if R 2 1 < 1 and R 1 2 < 1, one of the strains excludes the other based on the initial conditions. To determine the impact of within-host competition of the strains on the populationlevel competition, we investigate the impact of the within-host competition on the population-level invasion reproduction numbers. Figure 5 shows that the invasion reproduction number R 2 1 is a decreasing function of the withinhost reproduction number of strain two R 2 while R 1 2 is an increasing function of the within-host reproduction number of strain two R 2 . Furthermore the invasion reproduction numbers exhibit the same trends of change when the initial viral load for strain two increases. Both reproduction numbers experience a jump at the critical value of R 2 when R 2 = R 1 . The increase in R 2 affects much more significantly R 1 2 when R 2 < R 1 but when R 2 > R 1 R 2 1 is most sensitive to the change of R 2 . More importantly, if R 2 < R 1 then R 2 1 > 1 and R 1 2 < 1 signifying that strain one excludes strain two on within-host levels as well as on the population level. If, however, R 2 > R 1 then R 2 1 < 1 and R 1 2 > 1 leading to elimination of strain one on within-host levels as well as on the population level. This scenario, depicted in Figure 5, where the within-host reproduction numbers control the population level of the competition between the strains, was the most common scenario produced in simulations. As η 1 and η 2 vary other scenarios are also possible. For instance, another scenario   Figure 6. Plots of the invasion reproduction number R 2 1 and the invasion reproduction number R 1 2 with respect to within-host reproduction number of strain two R 2 when R 2 < R 1 . Here, η 1 = 0.02 and η 2 = 0.027. It is clear that for R 2 > 3 both invasion numbers are above one.
we produced was one where strain one remains dominant as the immunological reproduction number of strain two increases so that R 2 > R 1 . One uncommon and difficult scenario to produce was the case when both invasion reproduction numbers are larger than one with the parameters in Table 1 or their variations. However, we were able to find a case in which R 1 2 > 1 and R 2 1 > 1 with the parameters in Table 1 and η 1 = 0.02, η 2 = 0.027. This scenario is illustrated in Figure 6. Thus, Figure 6 demonstrates that coexistence in our model indeed occurs. The location of that coexistence equilibrium in the (Q 1 , Q 2 ) plane is given by the point (644, 6702).