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

1.   Introduction

A number of mathematical models of infectious diseases (e.g., human immunodeficiency virus (HIV)) have been studied at a single microscale or macroscale ^[1,2,3,4,5]. The single-scale models are proposed to either explore the within-host viral dynamics at the individual level and give guides for treatment strategies, or the between-host transmission dynamics at the population level to predict the future prevalence and suggest effective control measures. However, more evidence indicated that viral load at the individual level affects the progression of infection ^[6,7]. Laith J et al. found that a higher circulating viral load was positively related to a higher rate of host-to-host transmission ^[8]. This means that the dynamics at different scales are not independent, but interrelated. Obviously, this also complicates disease control. For example, implementing measures to change behavior patterns may not only affect the disease transmission at the population level, but even the probability of an individual receiving treatment, and further the within-host virus dynamics. Therefore, formulating a multiscale model that can combine the within-host and between-host scales and study the macro-micro interaction mechanism is significant for determining more effective strategies at different levels.

In recent years, many multiscale systems have been proposed in the area of mathematical biology ^{[9,10,11,12,13]}. Initially, some researchers built nested models to consider the evolution of the host and parasite by evolutionary dynamics ^[14,15]. Further, researchers applied coupled models to study infectious diseases. They integrated the within-host model into the epidemiological model by introducing the viral load-dependent transmission rate or disease-induced mortality rate to explore the potential effect of micro dynamics on the macro dynamics, design the coupled optimal scheme and provide the cost-effectiveness analysis ^[16,17,18]. Later, Feng et al. formulated a multiscale model with an environment compartment to investigate the effects of between-host dynamics on the viral progression within the hosts ^[19,20]. Taking contaminated environment as a coupling link, Wang et al. coupled the age-structured macro model and the micro within-host dynamics model through the bacteria-dependent indirect transmission ^[21], Sun et al. proposed a multiscale model with threshold control strategy to study the effect of threshold-dependent interventions on the spread of infectious disease ^[22], Xiao et al. linked the macroscale to the microscale in a spatiotemporal context to examine effects of an individual movement and spatial control measures on a disease outbreak ^[23].

It should be noted that most of the above-mentioned multiscale models were formulated by linking viral dynamics to between-host transmission dynamics, which can not reflect the influence of transmission of disease on viral loads within the hosts. It is worth noting that for some environmentally-driven infectious disease such as Toxoplasma gondii, there have been models successfully coupling the macro level to the microscale ^[19,20,21]. However, for some infectious diseases such as HIV that the virus can be spread by direct contact with infected individuals, how to construct a bridge which can couple the epidemiological model to the virus dynamics model and study the interactions between macro and micro level remain unclear and fall within the scope of our study.

The main purpose of this study is to propose a multiscale model which can bidirectionally couple the within-host viral dynamics and between-host transmission dynamics, and then analyze the dynamic behaviors of the full system. Note that most existing coupled models provided that between-host transmission rate is bilinear incidence, which is inconsistent with the fact that one can only contact with limited persons in a certain time. Hence, we include the standard incidence in our study. Specifically, in this article, we first formulate a multiscale model with standard incidence in section 2. In section 3 and section 4, we analyze the dynamic behaviors of the fast subsystem and slow subsystem, which can help us understand the dynamics of the full system. In section 5, we discuss the interactions of within-host and between-host dynamics by a series of numerical simulations. Finally, the main conclusions are highlighted.

2.   Model formulation

Here, we take HIV viral dynamics and transmission dynamics as an example to illustrate how we form the links on macro and micro levels. The micro system is embedded into the macro system by incorporating the viral load-dependent transmission and disease-induced mortality. The virus shed by the contacted infections, depending on the contact rate and prevalence, enter into the within-host viral system due to close contact and act as a bridge such that between-host system is coupled into the within-host system, inspired by the ideas by Kostova ^[24], Bhattacharya and Maia Martcheva ^[25]. In order to investigate the impact of macro epidemic dynamics on micro virus dynamics and vice verse, we proposed a multiscale model, motivated by the sexual contact transmission of HIV ^[26]. The model equations are as follows.

where $S(t), I(t), A(t)$ are the number of susceptibles, HIV-positive individuals without clinical manifestation and AIDS patients who have developed one or more opportunistic infections regardless of their CD4 count, $N(t) = S(t)+I(t)+A(t)$. $T(t), T^{*}(t),$ $V(t)$ denote the densities of healthy T cells, infected T cells and viral load, respectively. $\lambda, k, d$ represent the recruitment, per-capita infection rate, per-capital mortality of healthy cells, $\delta$ is the per-capital mortality of infected cells. $\eta$ is the drug efficacy ($0 \leq \eta \leq 1$). For HIV infection, it represents the effectiveness of reverse transcriptase (RT) inhibitors. $N_{1}$ is the virus production rate by an infected cell, $c$ is the clearance rate of virus within host. $c_ {1}, c_ {2}$ represent the average number of contacts with HIV-positive individuals and AIDS patients who have developed clinical symptoms, $p_ {1 }, p_ {2}$ denote the amount of virus released by HIV infections and AIDS patients at each contact. $\Pi, \mu$ denote recruitment rate and natural death rate of hosts. $\xi$ is the transfer rate from HIV-positive stage to the onset of opportunistic infections. $\beta_ {1} (V), \beta_ {2} (V)$ are the transmission probability of HIV-positive individuals and AIDS patients, $\alpha_{1}(V), \alpha_{2}(V)$ are the corresponding disease-related mortality. They are dependent on the viral loads, such as $\beta(V) = aV$, or $\beta(V) = \frac{aV}{1+bV}$, or $\beta(V) = aV^{q}, \ q < 1$ with

Since the system (2.1) contains the dynamics on different time-scales, i.e., within-host dynamics on a fast time-scale and between-host dynamics on a slow time-scale, it is a challenge to directly study the multiscale system. For this problem, some studies analyzed the full system on two time-scales, i.e., fast system and slow system ^[27,28]. The fast system is obtained by assuming that the slow variables are all constants and the slow system is defined by assuming that the fast system will tend to the stable equilibrium very quickly. This methods allows us to derive the analytical results and reveal some meaningful phenomena of the full system from the analyses of the fast and slow subsystems.

3.   Dynamics of the fast subsystem

For the within-host system to have a meaningful couple with the between-host system, we focus on the case $I, A > 0$ in the following analysis. We first introduce a slow time-scale $\tau = \epsilon t$ with $0 < \epsilon\ll1$ and let

Then the system (2.1) can be written as

Setting $\epsilon = 0$ and $\frac{d}{dt} = \cdot$, we get the fast subsystem

where $\hat{I}, \hat{A}$ are considered as constants. Let $E_{f}^* = (\tilde{T}^{*}, \tilde{V}, \tilde{T})$ denote a positive equilibrium of system (3.1). It is obvious that for $\hat{I} > 0$, the fast subsystem (3.1) has no disease-free equilibrium. In addition, we can get the following result.

Theorem 1 For $\hat{I} > 0$, the fast subsystem (3.1) has unique positive equilibrium $E_{f}^*$, which is globally asymptotically stable.

Proof. Let the right sides of system (3.1) be equal to 0 and then we get $\tilde{T}^* = \frac{\lambda-d\tilde{T}}{\delta}, \tilde{V} = \frac{(c_{1}p_{1}\hat{I}/\hat{N}+c_{2}p_{2}\hat{A}/\hat{N})+N_{1}(\lambda-d\tilde{T})}{c}$ and $\tilde{T}$ satisfying

where $a_{0} = k(1-\eta)N_{1}d, a_{1} = -[k(1-\eta)(N_{1}\lambda+c_{1}p_{1}\hat{I}/\hat{N}+c_{2}p_{2}\hat{A}/\hat{N})+cd], a_{2} = c\lambda$.

Solving the above equation, we get

with $a_{1}^2-4a_{0}a_{2} > 0$. Obviously, $\tilde{T}_{1} > \tilde{T}_{2} > 0$. For a positive equilibrium, it should satisfy $\tilde{T}^*, \tilde{V} > 0$, which corresponds to the condition $\tilde{T} < \frac{\lambda}{d}.$ It is easy to verify $\tilde{T}_{2} < \frac{\lambda}{d} < \tilde{T}_{1}$. Hence, the system (3.1) has unique positive equilibrium $E_{f}^* = (\tilde{T}^{*}, \tilde{V}, \tilde{T})$ with $\tilde{T} = \tilde{T}_{2}.$ This completes the proof.

In the following, we verify that $E_{f}^{*}$ is globally asymptotically stable. Formulating the following Lyapunov function

then the derivative of $W(t)$ is

and for all $T, T^*, V > 0$, $\frac{dW}{dt} = 0$ holds only at $E_{f}^{*}$. By LaSalle's invariant principle, $E_{f}^{*}$ is globally asymptotically stable.

4.   Dynamics of the slow subsystem

Writing system (2.1) with respect to the slow time-scale $\tau$, we have

Letting $\epsilon = 0$, $\frac{d}{d\tau} = '$, and dropping the bar for convenience, we obtain the slow subsystem

where $\tilde{V}(S, I, A)$ is the steady state of the fast subsystem. The initial condition for system (4.1) is $S_{0} > 0$, $I_{0}$ or $A_{0} > 0$.

By using the next generation method ^[29], we can calculate the basic reproduction number of the slow system (4.1) as follows

where

with $R_{f} > 1$. Here, $R_{f} = \frac{\lambda N_{1} k(1-\eta)}{cd}$ represents the basic reproduction number of isolated within-host system.

4.1.   Stability of disease-free equilibrium

In the following, we first analyze the stability of disease-free equilibrium for the slow subsystem (4.1).

Theorem 2 The disease-free equilibrium $E_{0}^{s} = (\frac{\Pi}{\mu}, 0, 0)$ of slow subsystem (4.1) is locally asymptomatically stable (LAS) for $R_{s} < 1$ and unstable for $R_{s} > 1$.

Proof. The characteristic equation of Jacobian matrix at $E_{0}^{s}$ is

It is easy obtained

Obviously, for $R_{s} < 1$, $\chi_{2}+\chi_{3} < 0$ and $\chi_{2}\chi_{3} > 0$, that is $\chi_{2} < 0, \chi_{3} < 0$. Thus, the disease-free equilibrium is locally asymptomatically stable. For $R_{s} > 1$, $\chi_{2}\chi_{3} < 0$, i.e., there is a positive solution for above characteristic equation. This demonstrates the disease-free equilibrium is unstable. Hence, we obtain the Theorem 2.

4.2.   Existence and number of positive equilibrium

If there exists the equilibrium for the slow subsystem, we let $E_{s}^{*} = (\tilde{S}, \tilde{I}, \tilde{A})$ denote a positive equilibrium of system (4.1). For the existence and stability of the endemic equilibrium, it is always difficult to be analyzed theoretically with viral load-dependent transmission rate, disease-induced mortality rate and standard incidence, simultaneously. Hence, we provide partial analytic proofs for the case where $\beta_{i}(V)$ satisfies the properties as mentioned in the Section 2 and $\alpha_{i}(V) = \alpha_{i}$ $(i = 1, 2)$. Then we extend the results with numerical simulations.

Let $\tilde{x}_{1} = c_{1}p_{1}\frac{\tilde{I}}{\tilde{N}}+c_{2}p_{2}\frac{\tilde{A}}{\tilde{N}}$, then the equilibrium of the fast subsystem can be written as

Let $\tilde{x}_{2} = c_{1}\beta_{1}(\tilde{V})\frac{\tilde{I}}{\tilde{N}}+c_{2}\beta_{2}(\tilde{V})\frac{\tilde{A}}{\tilde{N}}$, then the equilibrium $E_{s}^{*}$ of the slow subsystem can be written as

Further, $\tilde{x}_{1}$ can be described by $\tilde{x}_{2}$, that is

Thus, $\tilde{T} = \tilde{T}(\tilde{x}_{2}), \tilde{T}^* = \tilde{T}^{*}(\tilde{x}_{2}), \tilde{V} = \tilde{V}(\tilde{x}_{2})$. Using the equation $\tilde{x}_{2} = c_{1}\beta_{1}(\tilde{V})\frac{\tilde{I}}{\tilde{N}}+c_{2}\beta_{2}(\tilde{V})\frac{\tilde{A}}{\tilde{N}}$, we obtain

Let

then $\tilde{{x}}_{2}$ is the solution of equation $H(x_{2}) = 0$. In the following, we only need to identify the number of zeros of $H(x_{2})$ for $x_{2}\geq 0$ by examining the properties of the function $H$. It can be verified that

and $H(\infty) < 0$.

The derivative of $H(x_{2})$ is

where

We find that the sign of $H'(x_{2})$ can be determined by $G(x_{2})$. So, we first analyze the properties of the function $G(x_{2})$.

with

The result indicates that the function $G(x_{2})$ is a monotone decreasing function and $G(\infty) < 0$. Thus, if $G(0) > 0$, then with $x_{2}$ increasing, the function $H'(x_{2}) > 0$ firstly and then $H'(x_{2}) < 0$, which means that the function $H(x_{2})$ increases firstly and then decreases. In this case, if $H(0) < 0$ and $H_{max} > 0$, we can obtain function $H(x_{2})$ has two positive equilibriums; if $H(0) < 0$ and $H_{max} = 0$ or only $H(0) > 0$, there is unique positive equilibrium. If $G(0) < 0$, then the function $H(x_{2})$ always decreases. In this case, if $H(0) < 0$, then there is no positive equilibrium; if $H(0) > 0$, then there is only one positive equilibrium. In summary, we give the following results.

Theorem 3 Let $H = H(x_{2}), G = G(x_{2})$ be the function defined in (4.2) and (4.3), $H_{max} = {\max_{x_{2}}}$ $H(x_{2})$,

(1) For $R_{s} < 1, R_{f} > 1$;

$\; \; \; \; \; \;$(a) If $G(0) > 0$ and $H_{max} > 0$, there exists two positive equilibriums $E^{1*}_{s}, E^{2*}_{s}$;

$\; \; \; \; \; \;$(b) If $G(0) > 0$ and $H_{max} = 0$, there is only one positive equilibrium $E^{*}_{s}$;

$\; \; \; \; \; \;$(c) If $G(0)\leq0$ or $H_{max} < 0$, there is no positive equilibrium;

(2) For $R_{s} > 1, R_{f} > 1$, there is only one positive equilibrium $E^{*}_{s}$;

(3) For $R_{f} < 1$;

$\; \; \; \; \; \;$(a) If $G(0) > 0$ and $H_{max} > 0$, there exists two positive equilibriums $E^{1*}_{s}, E^{2*}_{s}$;

$\; \; \; \; \; \;$(b) If $G(0) > 0$ and $H_{max} = 0$, there is only one positive equilibrium $E^{*}_{s}$;

$\; \; \; \; \; \;$(c) If $G(0)\leq0$ or $H_{max} < 0$, there is no positive equilibrium.

To demonstrate that there are corresponding parameter values which can satisfy all conditions described in Theorem 3, we plot the curve $H(x_{2})$ for different parameter sets to show the three cases in which $H$ has 0, 1 or 2 zeros, that is the slow subsystem may have 0, 1 or 2 interior equilibriums. In Figure 1(A), we fix most parameters and change $R_{s}$ by parameter $c_{2}$. It shows for $0.5492 < R_{s} < 1$, the function $H(x_{2})$ has two zeros meaning there are two interior equilibriums for slow subsystem, but for $R_{s} < 0.5492$, there is no positive equilibrium. This result demonstrates $R_{s} = 0.5492$ is a threshold which determines the number of positive equilibrium. In Figure 1(B), we discuss the case of $R_{f} < 1$. Most parameters have the same values as those in Figure 1(A) except $N_{1}$. We change $R_{f}$ by parameter $N_{1}$. It shows $R_{f} = 0.7995$ is the lower bound such that the slow subsystem has two interior equilibriums.

4.3.   Local stability of the positive equilibrium

In the following analysis, we choose $\beta_{i}(V) = a_{i}V, \alpha_{i}(V) = \alpha_{i}, i = 1, 2$. Then, the Jacobian matrix of the slow subsystem (4.1) at $E_{s}^{*} = (\tilde{S}, \tilde{I}, \tilde{A})$ is

where

Let $\tilde{y} = 1+\frac{k(1-\eta)\tilde{x}_{1}+k(1-\eta)N_{1}\lambda+cd}{\sqrt{(k(1-\eta)\tilde{x}_{1}+k(1-\eta)N_{1}\lambda+cd)^2-4N_{1}dk(1-\eta)}},$ then

Being similar to above description, $\tilde{I}, \tilde{S}, \tilde{N}, \tilde{x}_{1}$ can be written as the functions of $\tilde{x_{2}}$. The characteristic equation of Jacobian matrix at $E_{s}^{*}$ is

where

According to Routh-Hurwitz criterion, we need to identify the sign of $j_{0}$, $j_{1}$ and $j_{1}j_{2}-j_{0}$. Specifically, if $j_{0} < 0$, $E_{s}^{*}$ is unstable, whereas if $j_{0} > 0$ and $j_{1} > 0, j_{1}j_{2}-j_{0} > 0$, $E_{s}^{*}$ is locally stability. It is very difficult to theoretically judge the sign of $j_{i}$ $(i = 1, 2, 3)$, due to the complexity of these functions, so we conduct a large number of numerical simulations. We discover that the Routh-Hurwitz criterion is only satisfied at the larger equilibrium for $R_{s} < 1, R_{f} > 1$ and is always satisfied at the unique equilibrium for $R_{s} > 1, R_{f} > 1$. We plot the equilibria and the corresponding stability for different cases in Figure 2: blue dashed line represents unstable equilibriums and black solid line denotes stable equilibriums and give the following conjectures: (1) for $R_{s} < 1, R_{f} > 1$, if $G(0) > 0$ and $H_{max} > 0$, the larger positive equilibrium is locally asymptotically stable, the other is unstable; (2) for $R_{s} > 1, R_{f} > 1$, the unique positive equilibrium is locally asymptotically stable; (3) for $R_{f} < 1$, if $G(0) > 0$ and $H_{max} > 0$, the larger positive equilibrium is locally asymptotically stable, the other is unstable.

4.4.   Backward bifurcation

We note that the existence of multiple endemic states and their stability may indicate the possibility of a backward bifurcation for our coupled system. In this section, we provide the detailed numerical studies about the parameter regions in which the bifurcation may occur and simulations illustrating the bistability of the disease-free equilibrium and endemic state. We plot the fraction of infected individuals at the positive steady state $\tilde{I}/\tilde{N}$ as a function of $R_{s}$ and $R_{f}$, where the solid and dashed line correspond to the stable and unstable equilibrium, separately. Figure 2(A) shows when $R_{s} < 1$, there is a lower bound $R_{sc}$ for $R_{s}$, above which the system has two positive equilibriums: one is stable, another is unstable and below which the system has no positive equilibria. This denotes the threshold, which governs the eradication of disease, is $R_{s} = R_{sc} = 0.5492$. In other words, the infectious disease can only be controlled until $R_{s} < R_{sc}$. This highlights the challenges in the control of disease. Similar to Figure 2(A), Figure 2(B) shows there exists a stable equilibria and an unstable equilibria for $R_{f}\in (R_{fc}, 1)$, which demonstrates a backward bifurcation can occur for $R_{f} < 1$. This result suggests that contact with infected individuals may cause the virus to persist in the population, although the virus can be cleared in the isolated individuals.

In Section 3 and 4, we analyze the dynamics of fast and slow subsystems, which are useful to obtain the theoretical results and reveal some meaningful phenomena of the full system. To confirm the analytical results, particularly the threshold condition for various dynamic behaviors of the full system, we plot the solution curve of the full system. Figure 3(A)-(B) show the time series of fraction of infected individuals for the case $R_{s} < 1, R_{f} > 1$ and $R_{f} < 1$, separately. We observe there exists bistable attractors for these two cases. The solutions converge to the disease-free equilibrium with lower initial values of $I(0)$, whereas the solutions converge to the interior equilibrium with higher initial values of $I(0)$. Figure 3(C)-(D) plot the time series of fraction of infected individuals and viral load for the case $R_{s} > 1, R_{f} > 1$, in which the system has a unique interior equilibrium. The result shows the two time scales for the fast variable $V$ and slow variable $I/N$. Specifically, the fast variable $V$ quickly approaches the value around its equilibrium, followed by the convergence of the slow variable $I/N$ at a slow rate. What's more, it should be noted that for the case $\alpha_{1}(V) = \alpha_{1}V, \alpha_{2}(V) = \alpha_{2}V$ with bilinear incidence, the system may appear richer dynamics except for the backward bifurcation. Specifically, given suitable parameter values for the case $R_{s} > 1, R_{f} > 1$, we plot the phase trajectories of the between-host system and the within-host system as shown in Figure 4, in which we can observe the appearance of a stable limit cycle for the full system.

5.   The interactions of within-host and between-host systems

In the following, we will discuss the interactions of within-host and between-host dynamics. To explore the influence of macro transmission dynamics on the micro virus dynamics, we plot the viral load at the equilibrium $\tilde{V}$ versus the fraction of infected individuals at equilibrium $\tilde{I}/\tilde{N}$, contact rate $c_{2}$ and viral releasing rate per contact $p_{2}$. Figure 5(A) shows that $\tilde{V}$ is an increasing function of $\tilde{I}/\tilde{N}$, which demonstrates the more infected individuals in the population, the higher the viral loads. In addition, as $R_{f}$ increases, $\tilde{V}$ also increases. Figure 5(B) gives the variation in the number and the size of equilibrium $\tilde{V}$ with contact rate $c_{2}$ increasing. For the case $R_{f} = 0.8582 < 1$, there exists a critical level $\bar{c}_{2} = 13.6$ for contact rate $c_{2}$, above which there are two positive equilibriums: one is stable, the other is unstable and below which there is no positive equilibria. Further, for $c_{2} > \bar{c}_{2}$, the larger the contact rate is, the higher the stable equilibrium level of viral loads is. The results suggest that for the coupled within-host dynamic model, $R_{f} = 1$ is not a threshold to govern the eradication of virus in the hosts. Contact with infected individuals may cause the virus to persist in the population, even if the virus can not be persistent in the isolated individuals. Note that for the case $R_{f} = 1.1013 > 1$, there is no positive equilibrium $\tilde{V}$ for $c_{2} < 5.3$. This result implies that having fewer contacts with infected individuals, the virus will eventually be eliminated at the population level even if there are some individuals who carry virus for long time. Figure 5(C) shows that as viral releasing rate $p_{2}$ increases, the value of stable equilibrium $\tilde{V}$ also increases in the case $R_{f} > 1$. This demonstrates that the final viral loads may tend to high level if contacts occur with infected individuals whose virus replication is in the active period. The results suggest that limiting contact with infected individuals by some protective measures, especially with those whose virus replication is active, is effective to control virus at a low level.

To investigate the influence of micro parameters on the macro disease transmission, we plot the prevalence at equilibrium $\tilde{I}/\tilde{N}$ against drug efficacy $\eta$ under two cases: case 1): $p_1 = 40, p_2 = 50$ and case 2): $p_1 = p_2 = 0$. In particular, the case 2) with $p_{1} = p_{2} = 0$ means decoupling the macro dynamics from micro model. With assumption $\alpha_{i}(V) = \alpha_{i}$, Figure 5(D) shows that increasing treatment efficacy $\eta$ leads to the decline of the prevalence $\tilde{I}/\tilde{N}$. Specifically, $\tilde{I}/\tilde{N}$ declines relatively slow when increasing the drug effectiveness initially until a critical level above which the prevalence begins to decline quickly. This indicates that relatively low treatment efficacy is not enough to effectively control the infectious disease. Compared with case 1) $p_1 = 40, p_2 = 50$, the same drug efficacy leads to a lower disease prevalence for the case 2). Furthermore, the threshold of disease elimination for drug efficacy also reduces from $\eta = 0.55$ to $\eta = 0.47$ in the absence of macro to micro coupling. This demonstrates that ignoring the impact of macro transmission dynamics on micro virus dynamics may underestimate the transmission of diseases in the population and consequently the requirement of drug effectiveness.

Generally, for controlling infectious disease, we should implement measures to reduce the basic reproduction number and consequently the new infections. In this model, we can easily observe that $R_{s}$ decreases with the value of macro parameters $a_{i}$, $c_{i}$ $(i = 1, 2)$ decreasing and the value of micro parameter $\eta$ increasing. However, the existence of backward bifurcation makes the control and elimination of infectious diseases more complicated. To explore the interactions between microscale and macroscale with $R_s\leq R_{sc}$, we plot the trends of contact rate $c_{2}$, transmission coefficient $a_{1}$ and drug effectiveness $\eta$ at surface $R_{s} = R_{sc}$ in Figure 6. If the values of these three parameters are below this surface (i.e., $R_{s} < R_{sc}$), then there is no positive equilibrium and the disease can die out. Otherwise, the disease may be persistent. Furthermore, Figure 6 shows that as $\eta$ increases, the parameter region of $(a_{1}, c_{2})$ in which disease doesn't outbreak enlarges. This demonstrates with high drug efficacy, the control measures with relatively low intensity can still eliminate the disease. However, for small drug effectiveness $\eta$, in order to eliminate the disease, we should not only reduce the infection rate, but also keep the contact rate at a low level.

6.   Discussion

Multiscale systems can simultaneously describe the virus dynamics at the individual level and the transmission dynamics at the population level and have become a research focus in recent years. However, there are many challenges to constitute the linkage from the macroscale to the microscale. In this study, we formulate a multiscale model which can link the slow dynamics for disease transmission to the fast dynamics for viral progression, in order to examine the interactions of between-host dynamics and within-host viral dynamics. In fact, coupling the macro dynamics to micro models has been successfully proposed in environmentally driven infections. But, for some infectious disease such as HIV, the virus can be transmitted by direct contact with infected individuals, which may affect the generation rate of virus in microscopic level. Hence, in our model, the influx of viruses from infected individuals to the within-host system, depending on the contact rate and prevalence, becomes a bridge such that between-host dynamics is coupled to the within-host viral dynamics. What's more, the micro system is embedded into the macro system by introducing the viral load-dependent transmission rate and disease-induced mortality rate.

It is worth noting that the factors such as drug efficacy, individual differences and person-to-person contact influence the within-host dynamics, thereby affecting the value of viral load at the equilibrium, which in turn affects the spread of the disease at the population level. Therefore, it is extremely necessary to formulate a coupled model to reflect this feedback mechanism. Moreover, due to the fact that one can only contact a limited number of individuals per unit time, we choose the standard incidence instead of mass action in the model, which is a advantage compared with most existing coupled models. Of course, this also increases the difficulty of theoretical analysis.

We analyze the formulated model by studying the dynamics of the system on different time-scales, i.e., fast subsystem and slow subsystem. We prove that if $I > 0$, the within-host system globally converges to the unique positive equilibrium, which is independent on the threshold $R_{f}$ obtained from the isolated microscopic model. For the slow subsystem, we find that the system can have zero, one or two interior equilibria, which depends on the magnitudes of $R_{s}$ and $R_{f}$. Specifically, for $R_{s} < 1, R_{f} > 1$ or $R_{f} < 1$, there exists a locally stable positive equilibrium and an unstable positive equilibrium. For $R_{s} > 1, R_{f} > 1$, there is only one locally stable interior equilibria. This demonstrates that the backward bifurcation may occur for the coupled system and $R_{s} = 1$ can not act as the threshold for infection control. Numerical simulations suggest that the full system exhibits the similar dynamic behaviors to the slow subsystem. For $R_{s} < 1, R_{f} > 1$ or $R_{f} < 1$, whether the disease dies out or persists depends on the initial condition.

We numerically examine the synergistic impact between the viral dynamics at the individual level and transmission dynamics at the population level. On the one hand, we observe increasing the drug effectiveness can decrease the prevalence of disease. Previous multiscale models that did not consider the impact of macroscale on microscale may overestimate the control effect of drug treatment on disease transmission. On the other hand, the contact rate may induce the occurrence of backward bifurcation and then affect the number and size of equilibria of viral dynamic model. Contact with infected individuals may cause the virus to persist in the population, even if the virus can not be persistent in the isolated individuals. This emphasizes limiting contact with infected individuals is effective in eliminating virus or controlling virus at a low level at the population level. The existence of backward bifurcation also illustrates the coupled system has more complicated dynamics, and consequently results in the elimination of disease more complex. The suggested control measures may include both the reduction in contact rates or transmission rate at the population level and the increase in drug efficacy at the individual level, and using these control measures together can effectively control the diseases.

In summary, we developed a novel bidirectional coupled model with standard incidence and investigated the synergistic impact between the macro transmission dynamics and micro virus dynamics, which is a advantage compared to most existing multiscale models. We derived some interesting results from this model. The coupled system may appear the backward bifurcation and Hopf bifurcation under certain conditions. Especially, contact rate or drug efficacy is essential to induce these complex dynamic behaviors. The conclusions improve our understanding for prevention and control of infectious disease, which can not be obtained from the pure between-host transmission or within-host viral dynamic models. However, due to the complexity of coupled system, it is challenging to include all progression and transmission routes of infection in the model while providing a complete theoretical analysis.

Although this article is mainly discussed in the context of HIV, the approaches we used are able to be applied more generally in other infectious diseases. Moreover, it should be noted we study the rich dynamics that may occur in our proposed model by theoretical analysis and numerical simulations. However, due to the lack of reliable information on parameter values such as $p_1, p_2,$ our simulations remain more qualitative. For example, the critical level of contact rate $c_{2}$ with two positive equilibria may be larger than the real situation. Hence, it is necessary to further discuss the possible dynamic phenomena and the corresponding parameter region in practical application. We leave it for future work.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (NSFC, 11631012).

Conflict of interest

All authors declare no conflicts of interest in this paper.