key: cord-0861475-83lffwrm
authors: Li, Xue-Zhi; Gao, Shasha; Fu, Yi-Ke; Martcheva, Maia
title: Modeling and Research on an Immuno-Epidemiological Coupled System with Coinfection
date: 2021-10-13
journal: Bull Math Biol
DOI: 10.1007/s11538-021-00946-9
sha: 7eeb5644bef82249b2cb97efaaeecf6a9d76d2b2
doc_id: 861475
cord_uid: 83lffwrm

In this paper, a two-strain model with coinfection that links immunological and epidemiological dynamics across scales is formulated. On the with-in host scale, the two strains eliminate each other with the strain having the larger immunological reproduction number persisting. However, on the population scale coinfection is a common occurrence. Individuals infected with strain one can become coinfected with strain two and similarly for individuals originally infected with strain two. The immunological reproduction numbers [Formula: see text] , the epidemiological reproduction numbers [Formula: see text] and invasion reproduction numbers [Formula: see text] are computed. Besides the disease-free equilibrium, there are strain one and strain two dominance equilibria. The disease-free equilibrium is locally asymptotically stable when the epidemiological reproduction numbers [Formula: see text] are smaller than one. In addition, each strain dominance equilibrium is locally asymptotically stable if the corresponding epidemiological reproduction number is larger than one and the invasion reproduction number of the other strain is smaller than one. The coexistence equilibrium exists when all the reproduction numbers are greater than one. Simulations suggest that when both invasion reproduction numbers are smaller than one, bistability occurs with one of the strains persisting or the other, depending on initial conditions.

the prevalence of measles and other infectious diseases reduced the number of Indians in Mexico from 30 million to 3 million. The black death (bubonic plague) had four large-scale epidemics in Europe (AD 600, 1346 (AD 600, -1350 (AD 600, , 1665 (AD 600, -1666 (AD 600, , 1720 (AD 600, -1722 , and each time caused a significant decline in the European population. And the outbreak of novel coronavirus at the end of 2019 in early 2020 has brought great impact to the global economy and people's lives. For a long time, mankind has waged an indomitable struggle against various infectious diseases and made brilliant achievements. For example, smallpox, which has been rampant for nearly a thousand years, has finally been eliminated. Further, leprosy, poliomyelitis, diphtheria, measles, pertussis and other diseases have been curbed.

However, completely conquering epidemics is still a task for the future. The world health report issued by the World Health Organization shows that infectious diseases are still the primary killer of humans. A variety of infectious diseases, such as cholera, dysentery, malaria, sexually transmitted diseases, tuberculosis, influenza and so on, are still prevalent all over the world. In addition, with the rapid development of livestock and poultry breeding industries and the prosperity of international and domestic trade, a variety of infectious diseases (such as avian influenza, mouth kick disease, etc.) which transmit between animals and affect humans are widely prevalent in the world, which has aroused the extensive attention of governments, medical experts and mathematicians.

In recent years, the research progress in infectious disease dynamics modeling, virus dynamics modeling and immune system dynamics modeling has been very rapid throughout the world. A large number of mathematical models are used to analyze the spread of various infectious diseases, virus infection, immune response and other issues. Most of these models are suitable for studying the general laws of various infectious diseases, viruses and immunity. Infectious disease dynamics mainly takes the total population of a country or a region as the research object. The population of the country or region is divided into several disjoint subclasses, such as susceptible class, infected class, moving out class and so on. Using a flow chart, the corresponding disease transmission model is established and studied. Immune system dynamics mainly takes an individual as the research object. Considering the infection of a virus of certain cells in this individual, this kind of cells in this individual is divided into healthy cells and infected cells, and the corresponding immune system dynamics is established and studied.

In reference (Levin and Pimentel 1981) , Levin and Pimental introduced the possibility of superinfection: when an individual infected by virus 1 enters into a contact with an individual who has been infected by virus 2, it is likely that, say, the first individual becomes infected also with virus 2; furthermore, it is assumed that in this case, the individual super-infected with virus 2 will lose the original virus; in other words, no individuals are infected with more than one virus at the same time. Intuitively, when virus 1 has the smaller reproduction number, superinfection can occur. A multi-scale immuno-epidemiological model of super-infection has been studied in Martcheva and Li (2013) . However, a more realistic situation is that individuals can be infected with more than one virus at the same time, that is, coinfection may occur. When an individual with virus 1 contacts an individual infected by virus 2, this individual will become coinfected by virus 1 and virus 2. Such infected individuals may continue to carry two viruses or become only one virus due to the competition between the two viruses. Coinfection and super-infection models both within-host and between-host have been studied before (Alizon 2013; Candela et al. 2009 ).

Immuno-epidemiology has been drawing a significant interest lately (Hellriegel 2001) , and the first immuno-epidemiological model of nested type has been developed by Gilchrist and Sasaki (2002) . Since then, a number of different multi-scale immuno-epidemiological models have been developed and investigated (Shen et al. 2015; Ratchford and Wang 2019; Hu and Lou 2016; Mideo et al. 2008; Feng et al. 2012; Xue and Xiao 2020; Martcheva and Pilyugin 2006; Martcheva 2011; Cai et al. 2017) . Typical questions addressed with immuno-epidemiological models are the impact of the within-host dynamics on the between-host dynamics (Martcheva 2011; Cai et al. 2017 ) and the questions related to the evolution of the pathogens (Alizon 2013; Mideo et al. 2008; ). Evolutionary questions often lead to multi-strain frameworks where the competitive exclusion (Dang et al. 2016 ) and coexistence (Martcheva and Li 2013) of strains is investigated. Super-infection (Martcheva and Li 2013) , coinfection (Candela et al. 2009 ) and mutation ) are all mechanisms that can lead to coexistence of pathogen variants.

In this paper, we strive to extend our results in Martcheva and Li (2013) . We formulate a two-strain immuno-epidemiological model coupling the immune process in human body, which was introduced by Webb et al. (2005) , and consider the transmission process of infectious diseases in the population with coinfection between strains. Immuno-epidemiological modeling of infectious diseases is a technique that links a with-in host mathematical model with a between-host mathematical model. In the next section, we formulate the with-in host model with coinfection and discuss the basic reproduction numbers, the invasion reproduction numbers and the existence of the coexistence equilibrium within the host. In Sect. 3, we propose a two-strain epidemiological model with coinfection. We discuss the disease-free equilibrium and single-strain equilibrium and their stability. In Sect. 4, we establish the existence of the interior equilibrium under the condition that the reproduction numbers are larger than one. Simulations are conducted in Sect. 5. Finally, we discuss the dynamic relationship between immunology and epidemiology in Sect. 6.

In this section, we formulate a two-strain with-in host model with coinfection. Within a host, the time variable is denoted by τ and signifies time-since-infection. We denote healthy cells by x and cells infected by the ith strain of the virus by y i . The coinfected cells are denoted by y. The number of free virions of strain i is denoted by V i .

In this model, uninfected cells are created at a constant rate r and are infected by interaction with free virions of strain i according to a mass action law with constant rate β i . Productively infected cells y i produce γ i > 1 new virions of strain i. The parameter α i denotes new virions of strain i produced by cells coinfected with two strains. The death rates for healthy cells, cells infected with strain i and cells coinfected by two strains are μ, d i and d, respectively. Here we assume d i and d are greater than μ. The (blue) is the between-host model. The population is divided into four classes: susceptible, infected by strain 1, infected by strain 2 and co-infected, denoted by S, I 1 , I 2 and J , respectively. The lower panel (green and orange) is the within-host model. The cell population (green) is divided into four classes: healthy, infected by strain 1, infected by strain 2 and co-infected, denoted by x, y 1 , y 2 and y, respectively. V 1 and V 2 represent virions with strains 1 and 2, respectively (orange). They are produced by infected cells. The viral load of infected people can affect the between-host transmission of the disease (purple line) (Color figure online) parameters δ i and s i denote the clearance rate and shedding rate for virions with strain i, respectively. Shedding rate is given by the amount of virus that leaves an infected individual and may serve to infect another individual. The parameter η 1 denotes the rate at which cells infected by strain 2 are simultaneously infected by the free virus of strain 1, and η 2 has a similar meaning. This with-in host model with coinfection is given by the following system of ordinary differential equations. A schematic diagram of the model is shown in Fig. 1 (lower panel) .

All parameters and dependent variables of this within-host model with coinfection and their definitions are found in Table 1 .

In this within-host model with coinfection, we can determine the infection-free equilibrium E 0 = (r /μ, 0, 0, 0, 0, 0). Computing the Jacobian of the model (2.1) at the infection-free equilibrium, we have

, where x * = r /μ. Considering J − λI = 0, we can expand the matrix in terms of the first column. This will give us an eigenvalue λ 1 = −μ. Then, we can expand the remaining matrix around the last row and obtain another eigenvalue λ 2 = −d. The remaining eigenvalues are the eigenvalues of the following matrices:

We can see that the above two matrices are the same except that the subscripts of the parameters are different. Therefore, we only consider J 1 . Now we can apply the usual conditions that guarantee that the eigenvalues of J 1 have negative real parts. We want Tr J 1 < 0 and Det J 1 > 0. Tr J 1 < 0 is obvious, and Det J 1 > 0 leads to the following condition:

(2.2) For J 1 and J 2 , the above inequality (2.2) can be written as

Therefore, the within-host reproduction number of strain i in model (2.1) is given by

We notice that if R i < 1, then Det J i > 0. If R i > 1, Det J i < 0, which implies there exists a positive eigenvalue. Consequently, if both R 1 < 1 and R 2 < 1 hold, the infection-free equilibrium is locally asymptotically stable. If R i > 1 for at least one i, then the infection-free equilibrium is unstable. Next we discuss the strain dominance equilibria.

Taking strain 1 as an example. Let the strain 1 dominance equilibrium E 1 = (x * 1 , y * 1 , V * 1 , 0, 0, 0). It satisfies the following equations:

The solution is

It can be seen that when R 1 > 1, E 1 exists. We compute the invasion reproduction number R 1 2 , which gives the number of secondary infections that one infected individual with the second strain will produce in a population in which the first strain is at equilibrium. From finding the eigenvalues of the Jacobian matrix at E 1 , we can obtain the following equation

That is

Dividing both sides of the equation by (λ + d) λ + d 2 + η 1 β 1 V * 1 , we have

Then, we observe that if λ increases, f (λ) is decreasing and g(λ) is increasing. Therefore, if

then the equation f (λ) = g(λ) does not have roots with non-negative real parts. We can prove this by contradiction. Suppose λ = ξ + iη with ξ ≥ 0, then

.

| for all λ with non-negative real part. In this case, the equilibrium E 1 is locally asymptotically stable. On the other hand, if f (0) > g (0), then the equation f (λ) = g(λ) has a real positive root. Therefore, the equilibrium E 1 is unstable.

We rewrite the inequality f (0) < g(0) as follows:

Multiplying both sides of the inequality by d 2 + η 1 β 1 V * 1 , we have

Transposing and simplifying, we get

Therefore, the invasion reproduction number is given by

We can see that R 1 2 < 1 is equivalent to f (0) < g(0) and R 1 2 > 1 is equivalent to f (0) > g(0). Similar results can be obtained for the strain 2 dominance equilibrium E 2 = (x * 2 , 0, 0, y * 2 , V * 2 , 0). Its stability depends on the relation between the following invasion reproduction number and one :

(2.5)

We summarize the above results in the following proposition:

If we do not remove healthy cells infections with V i and infected cells infections by another strain infected simultaneously with V i from the class of free virions, that is to say, we ignore the last two terms in the V 1 and V 2 equation, then the model (2.1) becomes

The within-host reproduction number of this system is

And the strain i dominance equilibrium has components

We can also calculate the invasion reproduction number of strain i in model (2.6). They are

Next we discuss the existence of coexistence equilibrium in model (2.6). When the within-host reproduction numbers and the invasion reproduction numbers are greater than one, there may be a coexistence equilibrium. The coexistence equilibrium E = (x, y 1 , V 1 , y 2 , V 2 , y) of model (2.6) satisfies the following equations:

The first three equations of the system (2.10) can be solved as follows:

From the fourth equation of the system (2.10), we have

Substituting y 1 , y 2 , y into the last two equation of the system (2.10), we get

Simplifying and substituting into x, we have

(2.11)

We need to judge whether the system (2.11) has a solution. If there is a solution, it means that coexistence equilibrium of the model (2.6) exists. Otherwise, there is no coexistence equilibrium.

Let the left-hand sides of Eq. (2.11) divided by the right-hand side are f 1 (V 1 , V 2 ) and f 2 (V 1 , V 2 ), respectively. Then, we can rewrite the system (2.11) as follows:

Then,

and

Therefore, we only need to prove

If we want to get (2.12), we need to consider the following two cases. Case 1: If

Here we use the truth that the death rate for infected cells d 1 is bigger than the natural death rate for healthy cells μ and the probability η 2 ≤ 1. Therefore,

This indicates that f 2 (V 1 , 0) > 1 for ∀V 1 . So there is a unique solution for all V 1 . Case 2: If

2 is the single strain 2 equilibrium. In fact,

We have

has two positive real solutions in the interval (0, V * 1 ). We can rewrite the above equation about V 1 as a quadratic equation of one unknown:

Therefore, C/A < 0, this means that there is only one positive solution for the quadratic equation of one unknown. This conclusion contradicts the hypothesis. Thus,

Then, (2.6) has a co-existence equilibrium.

In this section, we introduce a two-strain epidemiological model with coinfection. S(t) denotes the number of susceptible human individuals, where t is the chronological time. We stratify the infectious individuals by age-since-infection τ . Let i k (τ, t) be the density of individuals infected by strain k. Individuals in the class i k (τ, t) experience the same within-host dynamics given by model (2.1) with the initial condition of the other strain and the coinfected class set to zero. We use j(τ, t) to denote the density of individuals coinfected by the two stains. These individual experience withinhost dynamics given by (2.1). The two-strain immuno-epidemiological model with coinfection is given by the following system. The flowchart is shown in Fig. 1 .

Descriptions of the variables and parameters of the above model are found in Table 2 . In model (3.1), is the recruitment rate, and m 0 = m k (0) is the natural death rate of susceptible hosts at zero viral load. The infected hosts die at a variable rate dependent on their viral load m k (V k (τ )), where V k (τ ) is the viral load of the strain k. Hence, we will assume that the mortality of the infected hosts is

where k = 1, 2 denotes strain k and m k V k (τ ) is the additional host mortality due to the virus. The mortality rate of hosts coinfected by the two strains is

is also dependent on the within-host viral load. We suppose that β k (τ ) is proportional to the viral load at a given age-since-infection τ , that is,

The parameters q k denote the probability of successful transmission. The total population size is given by the sum of all classes:

The parameters P k take on values zero and one and control the direction of the co-infection. They are the modeling tools for scaling up the outcome of within-host competition to the between-host interactions. The values of P k rely on the size of the invasion reproduction numbers R 2 1 and R 1 2 . They can be expressed as follows:

If R 2 1 > 1 and R 1 2 > 1, P 1 = P 2 = 0. The parameters P 1 and P 2 are not probabilities, but a binary tool to capture two symmetric models in one. Flowchart of the model is given in Fig. 1 .

Description of variables and parameters in model The epidemiological reproduction number of strain k in system (3.1) is given by the following expression:

The system (3.1) has one disease-free equilibrium and two strain dominance equilibria. The disease-free equilibrium always exists and is given by

We assume that the stain 1 dominance equilibrium is E 1 = (S * 1 , i * 1 , 0, 0). Then E 1 satisfies the following equations:

From the second and third equations in Eq. (3.3), we obtain

Eliminating i * 1 (0), we then have

From the first and third equations in (3.3), we obtain

Dividing both sides of the equation by N * 1 and substituting

The solution is

Therefore, the components of the strain 1 dominance equilibrium are

Similarly, the strain 2 dominance equilibrium is

We consider the equilibrium (S * , 0, 0, j * ). But through our analysis of the model, we can see that this is impossible.

Next, we study 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 x(t) =xe λt , z 1 (τ, t) =z 1 (τ )e λt , z 2 (τ, t) = z 2 (τ )e λt , w(τ, t) =w(τ )e λt . Then, from the last two equations of the system (3.4), we havew(τ ) = 0. In addition, we obtain

where S * 0 /N * 0 = 1. The solution of the differential equation is

Replacingz k (τ ) in the equation forz k (0), we obtain the following characteristic equation for strain k:

Therefore, we have the following result. Consider λ with λ ≥ 0 and for such λ and each k, we have

Hence, λ < 0, in other words, the solutions of the equations H k (λ) = 1 have negative real part. Thus, the disease-free equilibrium is locally asymptotically stable. Now we assume

If λ = 0, and k is fixed, then we have H k (0) = R k > 1. In addition,

Therefore, the equation H k (λ) = 1 has a real positive root. That is, the disease-free equilibrium is unstable.

To study the local stability of the strain k dominance equilibrium E k for a fixed k, we recall that the equilibrium E k exists if and only if the basic reproduction number corresponding to this strain R k is greater than 1. So we assume that R k > 1 for a fixed k. What is more, the local stability of the strain k dominance equilibrium depends on the invasion reproduction numbers R 2 1 and R 1 2 for each of the strains. We can give the expression of the invasion reproduction numbers of the first and the second strains:

The invasion number of the strain two at the equilibrium of strain one R 1 2 gives the number of secondary infections one individual infected with strain two will produce in a population in which strain one is at equilibrium during her/his lifetime as infectious.

The results on the local stability of the strain k dominance equilibrium E k are summarized as follows Proposition 3.2 Assume for a fixed k, R k > 1. If R k i < 1, for i, k = 1, 2, and i = k, then the strain k dominance equilibrium E k is locally asymptotically stable. If R k i > 1, the strain k dominance equilibrium E k is unstable.

Proof In order to simplify the proof process, we will assume without loss of generality that k = 1, that is, R 1 > 1. We study the linearized equations around the strain one dominance equilibrium E 1 . We introduce the following notation for the perturbations

The linearized system around the strain one dominance equilibrium E 1 becomes

We look for exponential solutions of the form x(t) = xe λt , z k (τ, t) = z k (τ )e λt , w(τ, t) = w(τ )e λt . We are considering the case when the immune system is converging to strain 1 equilibrium, that is, P 1 = 0 and P 2 = 0. Substituting z 2 (τ, t) = z 2 (τ )e λt , w(τ, t) = w(τ )e λt into the differential equations of z 2 and w, we can obtain the following linear ordinary differential system

By solving Eq. (3.8), we can get

Substituting z 2 (τ ), w(τ ) into Z and simplifying, we get

We can get rid of Z on both sides of this equation and obtain the following characteristic equation

We denote the right-hand side of the equation with L 1 (λ), that is,

We can find that L 1 (λ) is a decreasing function of λ for λ is real. Moreover, lim λ→∞ L 1 (λ) = 0, and L 1 (0) = R 1 2 . If R 1 2 > 1, then Eq. (3.9) has at least one positive real root. Hence, the strain one dominance equilibrium E 1 is unstable. If R 1 2 < 1, for λ with λ ≥ 0, then we have

Therefore, all the eigenvalues of the equation L 1 (λ) = 1 have negative real part. In addition, the stability of E 1 also depends on the eigenvalues of the following system

Solving the differential equation, we obtain

To facilitate the calculation, we apply the notation

From the first and third equations in Eq. (3.10), we obtain

The solution is

Moreover, linearizing the equation for the total population size

Substituting (3.11) and (3.12) in the equation for z 1 (0) and canceling z 1 (0) from both sides of the resulting equation, we can obtain the characteristic equation for λ:

From β 1 (τ ) = c 1 V 1 (τ ) and m 1 (V 1 (τ )) = m 0 + m 1 V 1 (τ ) and integration by parts, we have

Then we can obtain

For λ = 0, we have

Besides, from the previous calculation, we know that

Because R 1 > 1 implies m 1 /c 1 < 1, we have

On the other hand, for λ with λ ≥ 0,

Therefore, the characteristic equation F 1 (λ) = 1 only has solutions with negative real parts. That is, the strain one dominance equilibrium E 1 is locally asymptotically stable. This completes the proof.

From the above discussion, we saw that if for exactly one strain the corresponding basic reproduction number and invasion reproduction number are greater than one, the disease becomes endemic and the single-strain equilibrium will be locally asymptotically stable. If all reproduction numbers are less than one, all strains are eliminated and the disease dies out. Moreover, this conclusion can also extended to the case of multiple strains.

In this section, we discuss the existence of an coexistence equilibrium for which both strains are present with coinfection. We consider the coexistence equilibrium when the two basic reproduction numbers R k and the two invasion reproduction numbers R k j corresponding to each strain are all greater than one. In this case, P 1 = P 2 = 0. We assume E * = (S * , i * 1 (τ ), i * 2 (τ ), j * (τ )) is the coexistence equilibrium of the system (3.1). The coexistence equilibrium satisfies the following equations

then the system (4.1) becomes

We solve Eq. (4.2) and obtain

Next, we substitute i * 1 (τ ), j * (τ ) in the expression of Q 1 and obtain

Canceling Q 1 , we can get

Similarly, we can also obtain

To determine the expression of N * , we note m 1 = q 1 m 1 , m 2 = q 2 m 2 , then

Simplifying the above equation, we get

It follows that

The solution is

Therefore, the coexistence equilibrium satisfies the following equations

where

Setting

Then, we need to show that the equation F 1 (Q 1 , Q 2 ) = 1 and F 2 (Q 1 , Q 2 ) = 1 have real solutions. From the above equations, we can observe that F 1 (Q 1 , Q 2 ) decreases with respect to Q 1 and F 2 (Q 1 , Q 2 ) decreases with respect to Q 2 . From the expression of S * and N * , we have

Since Q * 2 = (R 2 − 1)/ρ 2 = B 2 , then

This function is continuous.

To show that F 1 (Q 1 , Q 2 ) < 1 for any Q 2 and for Q 1 large enough, letting Q 1 → ∞, we obtain

We can find that 1 − m 0 ρ 1 + m 0 ρ 1 · q 1 C 1 α 1 < 1 if q 1 C 1 α 1 < 1. In fact,

That is, we have F 1 (∞, Q 2 ) < 1. We can obtain that F 1 (Q 1 , Q 2 ) = 1 has a positive solution.

Define G 2 (Q 2 ) = F 2 (h(Q 2 ), Q 2 ). We can have F 1 (0, Q 2 ) > 1. In this case, h(Q 2 ) is always defined and not equal to 0. Thus,

When the limit of Q 2 tends to infinity, the limit of Q 1 also tends to infinity. In this case, we can easily conclude that

That is, we can conclude that the equation G 2 (Q 2 ) = 1 has a positive solution. It is easy to get Q 2 corresponding to the value of Q 1 . To summarize, this establishes the existence of positive solution of the system (4.3). This conclusion can be stated as follows Proposition 4.1 Assume R 1 > 1, R 2 > 1, and R 2 1 > 1, R 1 2 > 1. Then, the system (3.1) has at least one coexistence equilibrium.

In this section, we provide numerical simulations to confirm and extend our analytical results. Most of the parameter values are chosen from our previous study (Martcheva and Li 2013 ) (see Tables 1 and 2) . Between-host disease transmission is our major concern. It is influenced by the viral load of the infected individuals, which is determined by the parameters of the within-host model. Therefore, we are interested in Table 1 . The initial condition is (x 0 , y 0 1 , y 0 2 , y 0 , V 0 1 , V 0 2 ) = (100, 000, 6, 6, 6, 200, 100 ). (a) shows that the dynamics for V 1 are almost the same for model (2.1) and (2.6). (b) shows the same result for V 2 the impact of within-host parameters on the immuno-epidemiological reproduction numbers as well as the dynamics of the between host model. We studied analytically two within-host models, namely the original model (2.1) and its simplified form (2.6). Figure 2 shows that there is almost no difference of the virion populations for two models. Hence, we will only consider the simplified model (2.6) in the following simulations. In all simulations, we use the same within-host initial condition (x 0 , y 0 1 , y 0 2 , y 0 , V 0 1 , V 0 2 )=(100000, 6, 6, 6, 200, 100). First, we choose the clearance rate for virions with strain one δ 1 as an example. We start from looking at the influence of δ 1 on within-host model. Figure 3 shows that strain 1 dominant equilibrium is locally asymptotically stable for δ 1 = 2, 3, 4, 5 and strain 2 dominant equilibrium is locally asymptotically stable for δ 1 = 6. This agrees with the analysis results since δ 1 = 2, 3, 4, 5 correspond to the invasion reproduction numbers R 1 2 < 1, R 2 1 > 1, but δ 1 = 6 leads to R 1 2 > 1, R 2 1 < 1. From Fig. 3a , d, we can see that as δ 1 increases from 2 to 5, the peak of infected cells with strain 1 y 1 and the peak of virus with strain 1 V 1 are delayed. Figure 3d also shows that the peak of V 1 decreases as δ 1 increases from 2 to 5.

Then, we look at the impact of δ 1 on between-host level. From Fig. 4a , we can see that in general as δ 1 increases, the between host basic reproduction number R 1 and the invasion reproduction number R 2 1 decrease, while R 2 and R 1 2 increase. Since R 1 , R 2 , R 1 2 and R 2 1 are important indicators of disease transmission, this implies increasing δ 1 will depress the transmission of strain 1, but is helpful for the transmission of strain 2. Figure 4b -e shows that as δ 1 increases from 2 to 5, that increase causes only tiny difference in all prevalences. These four cases result in strain 1 persisting. Another case δ 1 = 6 leads to strain 2 persisting. This is in accord with our analytical results. To be more specific, δ 1 = 2, 3 correspond to R 1 > 1 and R 2 < 1, δ 1 = 4, 5 correspond to R 1 > 1, R 2 > 1, R 2 1 > 1, R 1 2 < 1 (see Fig. 4a ). In both scenarios, strain 1 dominant equilibrium is locally asymptotically stable. However, δ 1 = 6 corresponds to R 1 > 1, R 2 > 1, R 2 1 < 1, R 1 2 > 1 (see Fig. 4a ), which means strain 2 dominant equilibrium is Table 1 . The initial condition is the same as Fig. 2 . For δ 1 = 2, 3, 4, 5, the within-host basic reproduction numbers R 1 > 1, R 2 > 1, the invasion reproduction numbers R 1 2 < 1, R 2 1 > 1, strain 1 dominant equilibrium is locally asymptotically stable. For δ 1 = 6, R 1 > 1, R 2 > 1, R 1 2 > 1, R 2 1 < 1, strain 2 dominant equilibrium is locally asymptotically stable (Color figure online) Tables 1 and 2 . The within-host initial condition is the same as Fig. 2 . The between-host initial condition is (S 0 , i 0 1 (τ ), i 0 2 (τ ), j 0 (τ )) = (1,000,000, 5, 5, 5). (a) shows that in general as δ 1 increases, R 1 , R 2 1 decrease, and R 2 , R 1 2 increase. (b)-(e) show that as δ 1 increases from 2 to 5, only tiny difference in all prevalences. These five cases result in strain 1 dominant equilibria. Another case δ 1 = 6 leads to strain 2 dominant equilibrium. As δ 1 increases, the prevalence of individuals infected with strain 1 decreases, the prevalence of individuals infected with strain 2 increases and the prevalence of all infected individuals decreases (Color figure online) locally asymptotically stable. From Fig. 4e , we can see that strain replacement happens as δ 1 increases. It shows that the strain 2 dominant state has a lower prevalence of total infected individuals than strain 1 dominant states. However, the influence of δ 1 is small. Similar results are observed for δ 2 , β 1 and β 2 (see Figs. 6e, 7b, d) .

Zooming in Fig. 4b , c, e, we can see that as δ 1 increases, the prevalence of individuals infected with strain 1 decreases (see Fig. 4b ), the prevalence of individuals infected with strain 2 increases (see Fig. 4c ) and the prevalence of all infected individuals decreases (see Fig. 4e ). From Fig. 4d , we can see that δ 1 = 5 leads to a high peak for the prevalence of co-infected individuals. One possible explanation is that δ 1 = 5 results in a high peak for co-infected cells in within-host level (see Fig. 3c ) and δ 1 = 5 is close to the intersection of R 2 1 and R 1 2 (see Fig. 4a ). Hence, the transient dynamic is relatively slow.

Zooming in Fig. 4a , we can see there exist some values for δ 1 satisfying both invasion reproduction numbers R 2 1 and R 1 2 are less than one (see Fig. 5a ). Interestingly, even though both δ 1 = 5.3 and 5.4 lead to strain 1 dominant equilibria in the within-host level (see Fig. 5b ), in between-host level they result in strain 1 persistence and strain 2 persistence, respectively (see Fig. 5c , d). From Fig. 5b , we can see that even though V 2 goes to zero eventually in both cases, there are several high peaks in the beginning. Besides, R 2 1 and R 1 2 are less than one, and the order of R 2 1 and R 2 1 changes at some value between δ 1 = 5.3 and 5.4. Hence, the dominant strain in within-host level and between-host level might be different.

In general, when both R 2 1 and R 1 2 are less than 1, bi-stability may happen. In other words, which strain will dominate depends on the initial conditions. For example, in Fig. 5e , f, we choose the same value for δ 1 = 5.364, but different initial conditions for between-host model, the system could result in strain 1 dominant or strain 2 dominant. This bi-stability is further verified in Fig. 6 .

From Fig. 6a , we can see that, in general as δ 2 increases, R 1 , R 2 1 increase, and R 2 , R 1 2 decrease. This indicates that increasing δ 2 will depress the transmission of strain 2, but is helpful for the transmission of strain 1. Similarly, increasing β 1 will depress the transmission of strain 2, but is helpful for the transmission of strain 1 (see Fig. 7a ). However, increasing β 2 will depress the transmission of strain 1, but is helpful for the transmission of strain 2 (see Fig. 7c ).

In this paper, we have formulated a two-strain model with coinfection. The model links the dynamics of an immune system of two strains and an epidemic model of two strains. Through proper analysis, we know that in the with-in host model, there is only one infection-free equilibrium in the case of R i < 1, i = 1, 2. When the basic reproduction number R i > 1, in addition to the infection-free equilibrium, there are single-strain equilibria-one single strain equilibrium corresponding to each strain. When both the basic reproduction numbers and the invasion reproduction numbers are greater than 1, there is a coexistence equilibrium in the model. We define the invasion numbers and we establish that strain i dominance equilibrium is locally asymptotically stable if the invasion reproduction number of strain j is smaller than . e, f Between-host prevalences for δ 1 = 5.364 for different between-host initial conditions. All the other parameters are obtained from Tables 1 and 2 . For all panels, the within-host initial condition is the same as Fig. 2 . The between-host initial condition for (c)-(e) is the same as Fig. 4 . The between-host initial condition 2 for (f) is (S 0 , i 0 1 (τ ), i 0 2 (τ ), j 0 (τ )) = (1,000,000, 0.01, 50, 0.01). (a) shows that there exist some values of δ 1 satisfying both invasion reproduction numbers are less then 1. b shows that for both δ 1 = 5.3 and 5.4, strain 1 dominant in within-host level. c, d show that in between-host level, strain 1 dominant equilibrium is stable when δ 1 = 5.3 and strain 2 dominant equilibrium is stable when δ 1 = 5.4. e, f show that when δ 1 = 5.364, which stain will dominant depends on the initial condition (Color figure online) Fig. 6 a, b Epidemiological basic reproduction numbers and invasion reproduction numbers for the model (3.1) with respect to δ 2 . c, d Prevalences for δ 2 = 13.685 with different between-host initial conditions. e Prevalences of all infected individuals for different values of δ 2 . All the other parameters are obtained from Table 1 and 2. For all panels, the within-host initial condition is the same as Fig. 2 . The between-host initial condition for (c, d) is the same as Fig. 5e, f. (a) shows that in general as δ 2 increases, R 1 , R 2 1 increase and R 2 , R 1 2 decrease. (b) shows that there exists some value of δ 2 satisfying both invasion reproduction numbers are less than 1. (c, d) shows that when δ 2 = 13.685, which strain will dominant depends on the initial condition. (e) shows that as δ 2 increases, strain replacement happens. The influence of δ 2 on the prevalence of all infected individuals is small (Color figure online) Fig. 7 a, c Epidemiological basic reproduction numbers and invasion reproduction numbers for the model (3.1) with respect to β 1 and β 2 , respectively. b, d Prevalences of all infected individuals for different values of β 1 and β 2 . All the other parameters are obtained from Tables 1 and 2 . For all panels, the within-host initial condition is the same as Fig. 2 . The between-host initial condition for (b) and (d) is the same as Fig. 4. (a) shows that in general as β 1 increases, R 1 , R 2 1 increase, and R 2 , R 1 2 decrease. (c) shows that in general as β 2 increases, R 1 , R 2 1 decrease, and R 2 , R 1 2 increase. (b) and (d) show that as β 1 increases or β 2 increases, strain replacement happens. The influence of β 1 and β 2 on the prevalence of all infected individuals is small (Color figure online) one. If the invasion reproduction number of strain j is greater than one, then the strain i dominance equilibrium is unstable.

We further define a multi-scale immuno-epidemiological model with two strains and coinfection. We structure each strain and the coinfected class with time-sinceinfection. Individuals in each of the infectious classes exhibit the same with-in host dynamics; however, if strain i is the only strain present, we start the strain j variables in the within-host model from zero. We study the disease-free and endemic equilibria of the immuno-epidemiological model. We further investigate the local asymptotic behavior of the immuno-epidemiological coupled system. We find that the diseasefree equilibrium of the coupled system is locally asymptotically stable if and only if R i < 1 for i = 1, 2. When R i > 1 and the corresponding invasion reproduction number R j i < 1, the strain i dominance equilibrium is locally asymptotically stable; when the invasion reproduction number R j i is greater than 1, the strain i dominance equilibrium is unstable. When all the reproduction numbers and invasion reproduction numbers are greater than 1, we also find that the coupled system has an interior equilibrium, that is, a coexistence equilibrium.

In simulation, we focus on the impact of within-host parameters on the between-host basic reproduction numbers, invasion reproduction numbers as well as the dynamics of the system. As δ 1 or β 2 increases, R 1 and R 2 1 decrease, and R 21 and R 1 2 increase. This indicates that increasing δ 1 or β 2 will depress the transmission of strain one, but is helpful for the transmission of strain two. This is reasonable since δ 1 is the clearance rate for the virions of strain one and β 2 represents infection rate of healthy cells infected by strain two. Similarly, we observe that increasing δ 2 or β 1 will depress the transmission of strain two, but is helpful for the transmission of strain one. As these parameters vary, strain replacement can happen, namely from stain one dominance to strain two dominance or vice versa. In addition, when both invasion reproduction numbers are smaller than one, bistability occurs with one of the strains persisting or the other, depending on initial conditions. However, the influence of within-host parameters on the prevalence of total infected individuals is very small.

Funding This work is supported partially by the National Natural Science Foundation of China (11771017). The research of M.M has been partially supported by NSF Grant DMS-1951595. 

Co-infection and super-infection models in evolutionary epidemiology

Modeling within-host evolution of HIV: mutation, competition and strain replacement

How does within-host dynamics affect population-level dynamics? Insights from an immuno-epidemiological model of malaria

Coinfection is an important factor in epidemiological studies: the first serosurvey of the aoudad (Ammotragus lervia)

Evaluating the importance of within-and between-host selection pressures on the evolution of chronic pathogens

Competitive exclusion in a multi-strain immuno-epidemiological influenza model with environmental transmission

Global asymptotic properties of a heroin epidemic model with treat-age

A model for coupling within-host and between-host dynamics in an infectious disease

Modeling host-parasite coevolution: a nest approach based on mechanistic models

Immunoepidemiology-bridging the gap between immunology and epidemiology

A nested model on HIV/AIDS, antiretroviral therapy and drug resistance

Virus dynamics: a global analysis

Evaluating the importance of within-and between-host selection pressures on the evolution of chronic pathogens

Lyapunov functional and global asymptotic stability for an infection-age model

An immuno-epidemiology model of paratuberculosis

An introduction to mathematical epidemiology

A Lyapunov-Schmidt method for detecting backward bifurcation in agestructured population models

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

An epidemic model structured by host immunity

Linking within-and between-host dynamics in the evolutionary epidemiology of infectious diseases

Modeling cholera dynamics at multiple scales: environmental evolution, between-host transmission, and within-host interaction

Global stability of an infection-age structured HIV-1 model linking within-host and between-host dynamics

Global dynamics and cost-effectiveness analysis of HIV pre-exposure prophylaxis and structured treatment interruptions based on a multi-scale model

Conflict and accord of optimal treatment strategies for HIV infection within and between hosts

Asymptotically autonomous differential equations in the plane

A model of antibiotic-resistant bacterial epidemics in hospitals

Analysis of a multiscale HIV-1 model coupling within-host viral dynamics and between-host transmission dynamics

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations