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

Abstract

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.

Keywords:

• immuno-epidemiology
• within-host dynamics
• age-since-infection
• multi-strain
• coexistence
• reproduction numbers
• invasion reproduction numbers

AMS Subject Classification:

• 92D30
• 92D40

1. 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 immuno-epidemiological 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 between- host 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-since-infection 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.

2. 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].

Table 1. List of parameters.

Notation	Meaning
x	Number of healthy cells
yi	Number of infected cells with i
Vi	Number of virions i
r	Production rate for healthy cells
βi	Infection rate of healthy cells infected by virus i
μ	Natural death rate for healthy cells
di	Death rate of infected cells with virus i
δi	Clearance rate for virus i
si	Shedding rate for virus i
νi	Number of virions produced by a single infected cell with strain i

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₀=(r/μ, 0, 0), and if ℝ_i>1, then the infected equilibrium E_i=(x_i*, y_i*, V_i*) exists with

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₀=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_jV_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 ℝ₁ and ℝ₂ corresponding to virus strain one and strain two. They can be expressed as follows: We note that if ℝ₂>ℝ₁ strain two super-infects individuals infected with strain one, that is P₁=1 and vice versa. The parameters P₁ and P₂ are not probabilities, but a binary tool to capture two symmetric models in one.

Model (3) is equipped with the following initial conditions: S(0)=S₀ and i_j(a, 0)=ϕ_j(a). All parameters are nonnegative and Λ>0, c_j>0, m₀>0.

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.

3. 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 where In addition, if the basic reproduction number , for each j there is a corresponding single-strain equilibrium given by where and 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 , let S(t)=S₀+x(t), i_j(τ, t)=z_j(τ, t). The linearized system becomes We look for solutions of the form x(t)=x̄ e^{λ t}, z_j(τ, t)=z̄_j(τ) e^{λ t}. For each j we obtain the following eigenvalue problem: Note that S₀*/N₀*=1. Solving the differential equation, we obtain and replacing z̄_j in the equation for z̄_j(0), we obtain the following characteristic equation for strain j: where Now we are ready to establish the following result.

Proposition 3.1

If then the disease-free equilibrium is locally asymptotically stable. If it is unstable.

Proof

Assume Consider λ with ℜ λ≥0. For such λ and each j we have Hence, the equations do not have a solution with nonnegative real part and the disease-free equilibrium is locally asymptotically stable.

Now assume For that fixed j and λ real we have . Furthermore, . Hence, the equation has a real positive root. Therefore, the disease-free equilibrium is unstable.   ▪

To investigate the local stability of the endemic equilibrium for a fixed j, we recall that strain j dominance equilibrium exists if and only if the basic reproduction number corresponding to this strain is greater than 1. So we assume that for a fixed j. The local stability of the strain j dominance equilibrium depends on the invasion reproduction numbers and for each of the strains. By definition, the invasion reproduction number of the first strain 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 are summarized below

Proposition 3.2

Assume for a fixed j, . Then the single-strain equilibrium is locally asymptotically stable if for i, j=1, 2, and i≠j. The single-strain equilibrium is unstable if .

Proof

To simplify the presentation, we will assume without loss of generality that j=1, that is . We study the linearized equation around the single-strain equilibrium . We introduce the following notation for the perturbations , , i₂(τ, t)=z₂(τ, t), N=N₁*+n. The system for the perturbations becomes 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 Substituting z̄₂(τ) in the equation for Z₂, dividing Z₂ from both side of the resulting equation, we get the following characteristic equation: We denote the right-hand side of that equation with , that is It can be seen that is a decreasing function of λ for real λ. Furthermore, and . If , then Equation (12) has at least one positive root, and therefore is unstable. If , for λ with ℜλ≥0, we have Hence, the equation has no solution with positive real part and all the eigenvalues of this equation have negative real part. Therefore, the stability of 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 we obtain 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₁(0) from both sides of the resulting equation, we obtain the following characteristic equation for where Here, we use the fact that From β₁(V₁(τ))=c₁V₁(τ) and m₁(V₁(τ))=m₀+m₁V₁(τ) it follows that That is, Evaluated for λ=0 the above equation also gives Note that Substituting above relations in Equation (17) we can rewrite the characteristic equation in the following form: For λ with ℜλ≥0 we have On the other hand, implies m₁/c₁<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 has only solutions with negative real parts. Thus, the single-strain equilibrium 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.

4. 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 and the two invasion reproduction numbers 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 ℝ₁>ℝ₂, in this case, P₁=0 and P₂=1. Assume is the co-existence equilibrium of the system, then it satisfies the following equations: Setting Then system (19) becomes First, we solve the second differential equation in (20) Solving the first differential equation in (20) and replacing i₂*(τ) in the resulting equation from (21), we obtain Substituting Equation (22) in the expression of Q₁ and canceling Q₁, we obtain Substituting Equation (21) in the expression of Q₂ and then cancelling Q₂, 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₁ and Q₂ we integrate the equations for i₁*(τ) and i₂*(τ) in system (20) and add them. We obtain Hence, This gives the following expression of N* in terms of Q₁ and Q₂: With expressions (25) and (27) for S* and N* in terms of Q₁ and Q₂ we obtain the following system in Q₁ and Q₂: To show the existence of an interior equilibrium, we have to show that this system has a nonzero solution (Q₁, Q₂). 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₂ can be expressed as a continuous function of Q₁: We use this expression for Q₂ to eliminate Q₂ from the first equation in system (28). Consequently, we obtain an equation in Q₁ only: , where The existence of an interior equilibrium for Equation (19) is now reformulated as the existence of at least one positive Q₁ satisfying for which Q₂ as computed from Equation (29) is nonzero. In what follows, we will prove this fact. Let Q₁=0, from Equation (29) we get the corresponding Q̄₂* and it can be expressed as We note that in the expression above we have . Hence, Q̄₂* is well defined and positive. It is easy to show that Let Q₂=0. From Equation (23) we obtain that the corresponding Q̄₁* can be expressed as We note that in the expression above, we have . Hence, Q̄₁* is well defined and positive. If Q₂=0, then from Equation (29) we have Now we define a function It is easy to see that , F′(Q₁)<0, and also Therefore, there exists a unique Qˆ₁* such that Note that Hence, Because F(Q₁) is a decreasing function, the above inequality implies that Q̄₁*<Qˆ₁*. Then, we have

Combining Equations (31) and (32) we conclude that the equation has a positive solution. It is easy to establish that for that value of Q₁, the corresponding Q₂ 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.

Proposition 4.1

Assume and . Then (3) has at least one co-existence equilibrium.

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

5.1. 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⁻¹. We assume the baseline value for the viral clearance of strain one δ₁. Faster viral clearance rates from the plasma have also been suggested, for instance δ₂≈ 23 days⁻¹ [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 ν₁=250, but De Boer et al. [10] suggest much larger numbers. Given the much larger clearance rate, we take ν₂=1000.

Perelson et al. [31] also estimate the infected cell life-span given by 1/d. It gives a value of d₁≈ 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₂=1. The uninfected CD4+ T-cell population in the human blood is given to be 5*10⁹ 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⁷ 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⁸. 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₁=0.00008, which we adopt for strain one. For strain two, we adopt twice as big shedding rate s₂=0.00015. The amount of virus in the semen can help estimate the initial viral load V₀. If all the amount of virus in the semen is absorbed and passed into the blood, then the initial V₀≈ 20,000 viral particles, considering that one viral particle has two RNA copies. We vary V₀ within the range (1000, 100,000). These parameters give an immune reproduction number of strain one of ℝ₁=6.22 and of strain two of ℝ₂=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⁹ which is consistent with what is reported in [35].

To estimate the epidemiological parameters, we take the probability of transmission q₁=q₂=1, so that c₁=0.00008 and c₂=0.00015. We define m₁ and m₂ 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₁=0.00002 and for m₂=0.00004.

To determine the population demographic parameters, we obtain the average human lifespan from [33] to be 70 years. Hence, m₀=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.

Table 2. Values of the within-host and epidemiological parameters.

Parameter	Value	Units	Refs
μ	0.01	day^−1	[36
r	5 * 10^7	Number of cells/day	[35
d1	0.5	day^−1	[31
δ1	3	day^−1	[31
β1	0.15*10^−10	1/{virion*day}	[35
ν1	250	virions	[35
s1	0.00008	day^−1	[32
q1	1	unitless	[19
m1	0.00002	1/virion*day	See text
η1	varies	unitless	—
d2	1	day^−1	[10
δ2	23	day^−1	[10
β2	0.15*10^−10	1/{virion*day}	[35
ν2	1000	virions	[35
s2	0.00015	day^−1	[32
q2	1	unitless	[19
m2	0.00004	1/virion*day	see text
η2	varies	unitless	—
Λ	275000	individuals/day	see text
m0	1/(365*70)	day^−1	[33

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

Figure 1. Epidemiological reproduction number vs. the infection rate of blood T cells β and the viral clearance rate δ.

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 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₀. 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). increases with s and decreases with m₀ and m₁. Figure 2 shows that as a function of s, the epidemiological reproduction number does not depend on the virus initial amount V₀.

Figure 2. Epidemiological reproduction number vs. the shedding rate.

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.

Figure 3. HIV prevalence vs. the infection rate of blood T cells β. HIV prevalence decreases with increased infection rate and viral load.

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 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 with respect to ν, β, d are positive, and the elasticity of with respect to δ is negative. The elasticities of the prevalence I₁ with respect to ν, β, d are negative, and the elasticity of I₁ 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₀. 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. Magnitude of the elasticities of ℛ₀ (left figure) and prevalence I₁ (right figure) with respect to within-host infection parameters. The elasticities are plotted as a function of the initial viral load V₀ and are colour-coded. Red, pink colours correspond to positive values, black and blue to negative values.

Figure 4 shows that the epidemiological reproduction number is less sensitive to the within-host infection parameters than the prevalence. Furthermore, the magnitudes of the elasticities are generally decreasing functions of the initial viral load V₀. 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.

5.3. 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 and , strain two excludes strain one; (2) if and , strain one excludes strain two; (3) if and , the two strains coexist; (4) if and , 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 population-level 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 is a decreasing function of the within-host reproduction number of strain two ℝ₂ while is an increasing function of the within-host reproduction number of strain two ℝ₂. 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 ℝ₂ when ℝ₂=ℝ₁. The increase in ℝ₂ affects much more significantly when ℝ₂<ℝ₁ but when ℝ₂>ℝ₁ is most sensitive to the change of ℝ₂. More importantly, if ℝ₂<ℝ₁ then and signifying that strain one excludes strain two on within-host levels as well as on the population level. If, however, ℝ₂>ℝ₁ then and 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 η₁ and η₂ vary other scenarios are also possible. For instance, another scenario we produced was one where strain one remains dominant as the immunological reproduction number of strain two increases so that ℝ₂>ℝ₁. 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 and with the parameters in Table 1 and η₁=0.02, η₂=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₁, Q₂) plane is given by the point (644, 6702).

Figure 5. Plots of the invasion reproduction number (left figure) and the invasion reproduction number (right figure) with respect to within-host reproduction number of strain two ℝ₂. Several variants of the invasion reproduction numbers are plotted for different values of the initial viral load of strain two V₀₂. Here, η₁=0.1 and η₂=0.005.

Figure 6. Plots of the invasion reproduction number and the invasion reproduction number with respect to within-host reproduction number of strain two ℝ₂ when ℝ₂<ℝ₁. Here, η₁=0.02 and η₂=0.027. It is clear that for ℝ₂>3 both invasion numbers are above one.

Acknowledgements

Maia Martcheva acknowledges partial support from NSF grant DMS-0817789 and grant DMS-1220342, Xue-Zhi Li acknowledges partial support from NSF of China grants 11271314 and 10911120387. This work was performed when Xue-Zhi Li was visiting the Department of Mathematics, University of Florida from January to March, 2012. Xue-Zhi Li greatly appreciates the hospitality of his hosts.