Modeling and Research on an Immuno-Epidemiological Coupled System with Coinfection.

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.

Introduction

Infectious diseases are caused by parasites (viruses, bacteria, etc.) and can spread from person to person. Infectious diseases have always been a great enemy of human health. Historically, infectious disease epidemics have brought great disasters to human survival, national economies and people’s livelihood. For example, during 1519–1530, 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–1350, 1665–1666, 1720–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; Martcheva 2015; 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; Coombs et al. 2007; Xue and Xiao 2020; Shen et al. 2019, 2019; 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; Coombs et al. 2007). 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 (Ball et al. 2007) 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.

The With-In Host Model with Coinfection

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 . The coinfected cells are denoted by y. The number of free virions of strain i is denoted by .

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 . Productively infected cells produce new virions of strain i. The parameter 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 , and d, respectively. Here we assume and d are greater than . The parameters and 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 denotes the rate at which cells infected by strain 2 are simultaneously infected by the free virus of strain 1, and 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.

Fig. 1

Flow diagram of the model (2.1) and (3.1). The upper panel (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, , 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, , and y, respectively. and 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)

Table 1

Description of variables and parameters in model (2.1)

Symbol	Description	Value	Unit	Source
Virable
x	Number of healthy cells
y_k	Number of infected cells with strain k=1,2
y	Number of co-infected cells with two strains
V_k	Number of virions of strain k=1,2
Parameter
r	Production rate for healthy cells	5×10^7	cells/day	Martcheva and Li (2013)
β_1	Infection rate of healthy cells infected by strain 1	0.15×10^-10	1/(virion∗day)	Martcheva and Li (2013)
β_2	Infection rate of healthy cells infected by strain 2	0.15×10^-10	1/(virion∗day)	Martcheva and Li (2013)
μ	Natural death rate for healthy cells	0.01	day^-1	Martcheva and Li (2013)
d_1	Death rate of infected cells with strain 1	0.5	day^-1	Martcheva and Li (2013)
d_2	Death rate of infected cells with strain 2	1	day^-1	Martcheva and Li (2013)
d	Death rate of co-infected cells	1.2	day^-1	Assumed
γ_1	Number of virions of strain 1 produced by a single infected cell with strain 1	250	virions	Martcheva and Li (2013)
γ_2	Number of virions of strain 2 produced by a single infected cell with strain 2	1000	virions	Martcheva and Li (2013)
α_1	Number of virions of strain 1 produced by a single co-infected cell	250	virions	Assumed
α_2	Number of virions of strain 2 produced by a single co-infected cell	1000	virions	Assumed
δ_1	Clearance rate for virions with strain 1	3	day^-1	Martcheva and Li (2013)
δ_2	Clearance rate for virions with strain 2	23	day^-1	Martcheva and Li (2013)
s_1	Shedding rate for virions with strain 1	0.00008	day^-1	Martcheva and Li (2013)
s_2	Shedding rate for virions with strain 2	0.00015	day^-1	Martcheva and Li (2013)
η_1	Probability of infected cell with stain 2 infected by strain 1	0.5	None	Assumed
η_2	Probability of infected cell with stain 1 infected by strain 2	0.5	None	Assumed

In this within-host model with coinfection, we can determine the infection-free equilibrium . Computing the Jacobian of the model (2.1) at the infection-free equilibrium, we havewhere . Considering , we can expand the matrix in terms of the first column. This will give us an eigenvalue . Then, we can expand the remaining matrix around the last row and obtain another eigenvalue . 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 . Now we can apply the usual conditions that guarantee that the eigenvalues of have negative real parts. We want and . is obvious, and leads to the following condition:For and , the above inequality (2.2) can be written asTherefore, the within-host reproduction number of strain i in model (2.1) is given byWe notice that if , then . If , , which implies there exists a positive eigenvalue. Consequently, if both and hold, the infection-free equilibrium is locally asymptotically stable. If 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 . It satisfies the following equations:The solution isIt can be seen that when , exists.

We compute the invasion reproduction number , 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 , we can obtain the following equationThat isDividing both sides of the equation by , we haveLetThen, we observe that if increases, is decreasing and is increasing. Therefore, ifthen the equation does not have roots with non-negative real parts. We can prove this by contradiction. Suppose with , thenandThus, for all with non-negative real part. In this case, the equilibrium is locally asymptotically stable. On the other hand, if , then the equation has a real positive root. Therefore, the equilibrium is unstable.

We rewrite the inequality as follows:Multiplying both sides of the inequality by , we haveTransposing and simplifying, we getTherefore, the invasion reproduction number is given byWe can see that is equivalent to and is equivalent to . Similar results can be obtained for the strain 2 dominance equilibrium . Its stability depends on the relation between the following invasion reproduction number and one (Martcheva 2015):We summarize the above results in the following proposition:

Proposition 2.1

If , then the boundary equilibrium exists. If , equilibrium is locally asymptotically stable. If , the equilibrium is unstable.

If we do not remove healthy cells infections with and infected cells infections by another strain infected simultaneously with from the class of free virions, that is to say, we ignore the last two terms in the and equation, then the model (2.1) becomesThe within-host reproduction number of this system isAnd the strain i dominance equilibrium has componentsWe can also calculate the invasion reproduction number of strain i in model (2.6). They areandNext 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 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 haveSubstituting into the last two equation of the system (2.10), we getSimplifying and substituting into x, we haveWe 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 and , respectively. Then, we can rewrite the system (2.11) as follows:andWe want to use to express . Since is a decreasing function of , if there is a solution, that solution is unique. There will be a solution if for arbitrary . We haveThen,andTherefore, we only need to prove . From , we can getIf we want to get (2.12), we need to consider the following two cases.

Case 1: Ifthat is,then we havesince is always true for real solutions. Here we use the truth that the death rate for infected cells is bigger than the natural death rate for healthy cells and the probability . Therefore,This indicates that for . So there is a unique solution for all .

Case 2: Ifwe haveThen, there exists a such that for and . Then, there is a unique for any such that .

We defineWhen , is the single strain 2 equilibrium. In fact,soWe haveLet , thenTo show , we assume that ; then, there are two values which satisfy .that is,has two positive real solutions in the interval .

We can rewrite the above equation about as a quadratic equation of one unknown:whereTherefore, , this means that there is only one positive solution for the quadratic equation of one unknown. This conclusion contradicts the hypothesis. Thus, we have and for . For any , there exists a unique such that .

Proposition 2.2

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

The Immuno-epidemiological Model with Coinfection

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 be the density of individuals infected by strain k. Individuals in the class 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 to denote the density of individuals coinfected by the two stains. These individual experience within-host 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.

Table 2

Description of variables and parameters in model (3.1)

Symbol	Description	Value	Unit	Source
S	Number of susceptible human individuals
i_k	Density of individuals infected by strain k=1,2
j	Density of individuals co-infected by the two strains
I_k	Number of individuals infected by strain k=1,2
J	Number of individuals co-infected by the two strains
N	Number of total population
Parameters
Λ	Recruitment rate	275,000	individuals/day	Martcheva and Li (2013)
m_0	Natural death rate of susceptible hosts	1/(365∗70)	day^-1	Martcheva and Li (2013)
m_1	Additional host mortality rate for infected individual with strain 1 due to one virus with strain 1	0.00002	1/(virion∗day)	Martcheva and Li (2013)
m_2	Additional host mortality rate for infected individual with strain 2 due to one virus with strain 2	0.00004	1/(virion∗day)	Martcheva and Li (2013)
m~_1	Additional host mortality rate for co-infected individual due to one virus with strain 1	0.00002	1/(virion∗day)	Assumed
m~_2	Additional host mortality rate for co-infected individual due to one virus with strain 2	0.00004	1/(virion∗day)	Assumed
c_1	Coefficient for transmission rate of strain 1	0.00008	day^-1	Martcheva and Li (2013)
c_2	Coefficient for transmission rate of strain 2	0.00015	day^-1	Martcheva and Li (2013)
q_1	Probability of successful transmission for strain 1	0.5	None	Assumed
q_2	Probability of successful transmission for strain 2	0.5	None	Assumed
θ_1	Recovery rate from infection with strain 1 for co-infected individuals	1/365	day^-1	Assumed
θ_2	Recovery rate from infection with strain 2 for co-infected individuals	1/365	day^-1	Assumed
ρ_12	Parameter to adjust the rate of co-infection for individuals infected with strain 1	0.8	None	Assumed
ρ_21	Parameter to adjust the rate of co-infection for individuals infected with strain 2	0.8	None	Assumed

In model (3.1), is the recruitment rate, and 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 , where is the viral load of the strain k. Hence, we will assume that the mortality of the infected hosts iswhere denotes strain k and is the additional host mortality due to the virus. The mortality rate of hosts coinfected by the two strains is . The transmission coefficient of the disease is also dependent on the within-host viral load. We suppose that is proportional to the viral load at a given age-since-infection , that is,The parameters denote the probability of successful transmission. The total population size is given by the sum of all classes:The parameters 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 rely on the size of the invasion reproduction numbers and . They can be expressed as follows:If and , . The parameters and are not probabilities, but a binary tool to capture two symmetric models in one. Flowchart of the model is given in Fig. 1.

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 bywhereWe assume that the stain 1 dominance equilibrium is . Then satisfies the following equations:From the second and third equations in Eq. (3.3), we obtainReplacing in the equation for , we haveEliminating , we then haveBecausewe haveFrom the first and third equations in (3.3), we obtainDividing both sides of the equation by and substituting and , we can getThe solution iswhere . Therefore, the components of the strain 1 dominance equilibrium arewhereSimilarly, the strain 2 dominance equilibrium iswhereandWe consider the equilibrium . 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 , let . The linearized system isWe look for solutions of the form . Then, from the last two equations of the system (3.4), we have . In addition, we obtainwhere . The solution of the differential equation isReplacing in the equation for , we obtain the following characteristic equation for strain k: , whereTherefore, we have the following result.

Proposition 3.1

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

Proof

AssumeConsider with and for such and each k, we haveHence, , in other words, the solutions of the equations have negative real part. Thus, the disease-free equilibrium is locally asymptotically stable.

Now we assumeIf , and k is fixed, then we have . In addition,Therefore, the equation has a real positive root. That is, the disease-free equilibrium is unstable.

To study the local stability of the strain k dominance equilibrium for a fixed k, we recall that the 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 k. What is more, the local stability of the strain k dominance equilibrium depends on the invasion reproduction numbers and for each of the strains. We can give the expression of the invasion reproduction numbers of the first and the second strains:whereThe invasion number of the strain two at the equilibrium of strain one 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 are summarized as follows

Proposition 3.2

Assume for a fixed k, . If , for , and , then the strain k dominance equilibrium is locally asymptotically stable. If , the strain k dominance equilibrium is unstable.

In order to simplify the proof process, we will assume without loss of generality that , that is, . We study the linearized equations around the strain one dominance equilibrium . We introduce the following notation for the perturbations . The linearized system around the strain one dominance equilibrium becomesWe look for exponential solutions of the form . We are considering the case when the immune system is converging to strain 1 equilibrium, that is, and . Substituting into the differential equations of and w, we can obtain the following linear ordinary differential systemLetBy solving Eq. (3.8), we can getNoteSubstituting into Z and simplifying, we getWe can get rid of Z on both sides of this equation and obtain the following characteristic equationWe denote the right-hand side of the equation with , that is,We can find that is a decreasing function of for is real. Moreover,and . If , then Eq. (3.9) has at least one positive real root. Hence, the strain one dominance equilibrium is unstable. If , for with , then we haveTherefore, all the eigenvalues of the equation have negative real part. In addition, the stability of also depends on the eigenvalues of the following systemSolving the differential equation, we obtainTo facilitate the calculation, we apply the notationFrom the first and third equations in Eq. (3.10), we obtainThe solution isMoreover, linearizing the equation for the total population sizewe havethat is,Substituting (3.11) and (3.12) in the equation for and canceling from both sides of the resulting equation, we can obtain the characteristic equation for : , whereFrom and and integration by parts, we haveThen we can obtainFor , we haveBesides, from the previous calculation, we know thatSubstituting above result in Equation , we getBecause implies , we haveOn the other hand, for with ,Therefore, the characteristic equation only has solutions with negative real parts. That is, the strain one dominance equilibrium 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.

Existence of Coexistence Equilibrium

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 and the two invasion reproduction numbers corresponding to each strain are all greater than one. In this case, . We assume is the coexistence equilibrium of the system (3.1). The coexistence equilibrium satisfies the following equationswhere .

Definethen the system (4.1) becomesWe solve Eq. (4.2) and obtainFirst, substituting in the expression of , we haveNext, we substitute in the expression of and obtainLetthenCanceling , we can getSimilarly, we can also obtainTo determine the expression of , we note , thenSimplifying the above equation, we getIt follows thatThe solution isTherefore, the coexistence equilibrium satisfies the following equationswhereSettingThen, we need to show that the equation and have real solutions. From the above equations, we can observe that decreases with respect to and decreases with respect to . From the expression of and , we haveLet , . Hence,Since , thenWe know that is a decreasing function of . If , we have when and . Suppose for . In this interval, for , that is, the equation have no solution. Let , then we have . DefineThis function is continuous.

To show that for any and for large enough, letting , we obtainWe can find that if . In fact,That is, we have . We can obtain that has a positive solution.

Define . We can have . In this case, is always defined and not equal to 0. Thus,When the limit of tends to infinity, the limit of also tends to infinity. In this case, we can easily conclude thatThat is, we can conclude that the equation has a positive solution. It is easy to get corresponding to the value of .

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 , and . Then, the system (3.1) has at least one coexistence equilibrium.

Simulation

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 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 ()=(100000, 6, 6, 6, 200, 100).

Fig. 2

Prevalences of virions in within-host model (2.1) and its simplified form (2.6). All the parameters are obtained from Table 1. The initial condition is . (a) shows that the dynamics for are almost the same for model (2.1) and (2.6). (b) shows the same result for

First, we choose the clearance rate for virions with strain one as an example. We start from looking at the influence of on within-host model. Figure 3 shows that strain 1 dominant equilibrium is locally asymptotically stable for and strain 2 dominant equilibrium is locally asymptotically stable for . This agrees with the analysis results since correspond to the invasion reproduction numbers , , but leads to , . From Fig. 3a, d, we can see that as increases from 2 to 5, the peak of infected cells with strain 1 and the peak of virus with strain 1 are delayed. Figure 3d also shows that the peak of decreases as increases from 2 to 5.

Fig. 3

Prevalences of infected cells and virions in within-host model (2.6) for different values of clearance rate for virions with strain 1 . All the other parameters are obtained from Table 1. The initial condition is the same as Fig. 2. For , the within-host basic reproduction numbers , , the invasion reproduction numbers , , strain 1 dominant equilibrium is locally asymptotically stable. For , , , , , strain 2 dominant equilibrium is locally asymptotically stable (Color figure online)

Then, we look at the impact of on between-host level. From Fig. 4a, we can see that in general as increases, the between host basic reproduction number and the invasion reproduction number decrease, while and increase. Since , , and are important indicators of disease transmission, this implies increasing will depress the transmission of strain 1, but is helpful for the transmission of strain 2.

Fig. 4

a Epidemiological basic reproduction numbers , and invasion reproduction numbers , for the model (3.1) with respect to . b–e The prevalences of individuals infected with strain 1, strain 2, co-infected individuals and all infected individuals for different values of . All the other parameters are obtained from Tables 1 and 2. The within-host initial condition is the same as Fig. 2. The between-host initial condition is . (a) shows that in general as increases, , decrease, and , increase. (b)–(e) show that as increases from 2 to 5, only tiny difference in all prevalences. These five cases result in strain 1 dominant equilibria. Another case leads to strain 2 dominant equilibrium. As 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)

Figure 4b–e shows that as increases from 2 to 5, that increase causes only tiny difference in all prevalences. These four cases result in strain 1 persisting. Another case leads to strain 2 persisting. This is in accord with our analytical results. To be more specific, correspond to and , correspond to , , , (see Fig. 4a). In both scenarios, strain 1 dominant equilibrium is locally asymptotically stable. However, corresponds to , , , (see Fig. 4a), which means strain 2 dominant equilibrium is locally asymptotically stable. From Fig. 4e, we can see that strain replacement happens as 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 is small. Similar results are observed for , and (see Figs. 6e, 7b, d).

Fig. 6

a, b Epidemiological basic reproduction numbers and invasion reproduction numbers for the model (3.1) with respect to . c, d Prevalences for with different between-host initial conditions. e Prevalences of all infected individuals for different values of . 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 increases, , increase and , decrease. (b) shows that there exists some value of satisfying both invasion reproduction numbers are less than 1. (c, d) shows that when , which strain will dominant depends on the initial condition. (e) shows that as increases, strain replacement happens. The influence of 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 and , respectively. b, d Prevalences of all infected individuals for different values of and . 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 increases, , increase, and , decrease. (c) shows that in general as increases, , decrease, and , increase. (b) and (d) show that as increases or increases, strain replacement happens. The influence of and on the prevalence of all infected individuals is small (Color figure online)

Zooming in Fig. 4b, c, e, we can see that as 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 leads to a high peak for the prevalence of co-infected individuals. One possible explanation is that results in a high peak for co-infected cells in within-host level (see Fig. 3c) and is close to the intersection of and (see Fig. 4a). Hence, the transient dynamic is relatively slow.

Zooming in Fig. 4a, we can see there exist some values for satisfying both invasion reproduction numbers and are less than one (see Fig. 5a). Interestingly, even though both 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 goes to zero eventually in both cases, there are several high peaks in the beginning. Besides, and are less than one, and the order of and changes at some value between and 5.4. Hence, the dominant strain in within-host level and between-host level might be different.

Fig. 5

a Zoom in part of Fig. 4a. b Within-host virus load for and 5.4. c, d Between-host prevalences for and 5.4. e, f Between-host prevalences for 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 . (a) shows that there exist some values of satisfying both invasion reproduction numbers are less then 1. b shows that for both 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 and strain 2 dominant equilibrium is stable when . e, f show that when , which stain will dominant depends on the initial condition (Color figure online)

In general, when both and 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 , 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 increases, , increase, and , decrease. This indicates that increasing will depress the transmission of strain 2, but is helpful for the transmission of strain 1. Similarly, increasing will depress the transmission of strain 2, but is helpful for the transmission of strain 1 (see Fig. 7a). However, increasing will depress the transmission of strain 1, but is helpful for the transmission of strain 2 (see Fig. 7c).

Discussion

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 , . When the basic reproduction number , 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 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-since-infection. 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 disease-free equilibrium of the coupled system is locally asymptotically stable if and only if for . When and the corresponding invasion reproduction number , the strain i dominance equilibrium is locally asymptotically stable; when the invasion reproduction number 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 or increases, and decrease, and and increase. This indicates that increasing or will depress the transmission of strain one, but is helpful for the transmission of strain two. This is reasonable since is the clearance rate for the virions of strain one and represents infection rate of healthy cells infected by strain two. Similarly, we observe that increasing or 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.