Multiscale model within-host and between-host for viral infectious diseases

Multiscale models possess the potential to uncover new insights into infectious diseases. Here, a rigorous stability analysis of a multiscale model within-host and between-host is presented. The within-host model describes viral replication and the respective immune response while disease transmission is represented by a susceptible-infected model. The bridging of scales from within- to between-host considered transmission as a function of the viral load. Consequently, stability and bifurcation analyses were developed coupling the two basic reproduction numbers R0W\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0^{W}$$\end{document} and R0B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0^{B}$$\end{document} for the within- and the between-host subsystems, respectively. Local stability results for each subsystem, including a unique stable equilibrium point, recapitulate classical approaches to infection and epidemic control. Using a Lyapunov function, global stability of the between-host system was obtained. Our main result was the derivation of the R0B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0^{B}$$\end{document} as an increasing function of R0W\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0^{W}$$\end{document}. Numerical analyses reveal that a Michaelis–Menten form based on the virus is more likely to recapitulate the behavior between the scales than a form directly proportional to the virus. Our work contributes basic understandings of the two models and casts light on the potential effects of the coupling function on linking the two scales.


Introduction
Infectious diseases remain a global concern despite the advance of modern medicine and living conditions (Saker et al. 2004). Combating with infectious diseases requires multidisciplinary efforts, among which include holistic understandings of the infection mechanisms (Gates 2015). Infectious disease dynamics, however, are governed by many interconnected scales, from complex within-host infection processes, to between hosts, vectors, and environments. In this context, mathematical modeling is an essential tool to untangle the processes, providing understandings of diseases dynamics and public health policies (Heesterbeek et al. 2015).
Typical infectious diseases models restrict their dynamics to one of the two scales (Feng et al. 2012;Murillo et al. 2013): within-host, focusing on cellular interactions; and between-host, focusing on transmission and infection statuses (Fig. 1). This restriction allows more analytical tractability (Feng et al. 2012) but there are problems where it becomes necessary to consider both scales (Heesterbeek et al. 2015;Handel and Rohani 2015). The last decade has witnessed an emergence of models that bridge within-and between-host disease processes (Gandolfi et al. 2014;Heesterbeek et al. 2015;Heffernan and Keeling 2009;Legros and Bonhoeffer 2016;Longini et al. 2005). This approach has served as a framework to understand the within-host evolution (Alizon et al. 2011;Mideo et al. 2008), epidemic control and prevention (Heffernan and Keeling 2009;Longini et al. 2005;Lukens et al. 2014), and beyond (Cen et al. 2014;Handel and Rohani 2015;Murillo et al. 2013). For example, natural selection of a disease can be investigated by linking pathogen dynamics to population dynamics (Gilchrist and Coombs 2006). In Heffernan and Keeling (2009), the authors investigated the immune responses dynamics in accompany with a SIR model with vaccination. In another study (Legros and Bonhoeffer 2016), the vector-borne transmission was analyzed with an extended SI model and a within-host subsystem.

Fig. 1
A disease model is usually restricted to one of the following scales: within a host (left) and between hosts (right). Our multiscale model considers both scales, by focusing on the interaction of the virus and immune response (e.g., CD8 + T cells) within an infected host, and direct contact between susceptible and infected individuals Applications of multiscale models have broadened the disease modeling landscape and have the potential to accurately describe disease dynamics (Murillo et al. 2013).
In this paper, we considered a nested model assuming a single coupling direction from within-to between-host. The within-host dynamics is taken from the model formulated by Boianelli et al. (2015) whereas between-host transmission dynamics are governed by the susceptible-infected (SI) model (Anderson and May 1992). In addition to the quantitative analysis in Boianelli et al. (2015), our results characterize the withinhost basic reproduction number R W 0 with the existence of a unique positive steady state having viral load V * (i.e., chronic infection). In addition, we obtain the between-host basic reproduction R B 0 assuming the "limiting case" of the between-host subsystem where the viral load remains at the chronic steady state (Gilchrist and Coombs 2006). The novel contribution of our work is an explicit relationship between the two reproduction numbers. Our theoretical analyses are accompanied by numerical examples.

Materials and methods
We considered two models that cover the infection dynamics within a host and the transmission between hosts. Each scale is described below with a two-dimensional ordinary differential equation.

Within-host subsystem
We focused on dynamics of virus (V ) and T-cell populations (E) as follows : (1) In this model, the virus replicates at a self-limiting rate determined by the intrinsic rate p and the carrying capacity K V . The virus is cleared at a per-capita rate proportional to T-cell populations (c V E). In the absence of virus, T cells stay at a homeostatic level given by E 0 := N E /δ E , where N E and δ E denote the replenishment rate and half-life of T cells, respectively. Otherwise, T-cells proliferate following a continuous function G(V ) of the viral load. To allow generality for our within-host subsystem, we imposed the following properties: A special case of G can be the Michaelis-Menten (saturation) form ) which was considered in our analyses. The parameter r represents the maximum rate of T cell proliferation while K E is the half-saturation constant.

Between-host subsystem
We considered the Susceptible-Infected (SI) model (Anderson and May 1992), where the host population is divided into susceptible and infected classes with densities S and I , respectively. This model reads as where β is the transmission rate. Susceptible individuals enter the population at a rate N S , and both subpopulations have distinct half-life parameters δ S and δ I . The host population attains the equilibrium size S 0 := N S /δ S in infection-free situations.

Bridging scales
Considering experimental observations of the impact of viral load on disease transmission (Blaser et al. 2014;Edenborough et al. 2012), we assumed a multiscale model coupling (1) and (4) where β = β(V ) has the following properties: As the exact functional relationship between viral load and the transmission rate can be debatable (Handel and Rohani 2015), we considered three different functional forms of the coupling function, including a linear, logistic, and saturation function that satisfies (5) in our numerical analyses. In particular, . The parameter r W is the rate of transmission due to the viral load. The parameter K W is a threshold of the viral load that a host may need to cross to transmit the infection.

Numerical example
To complement our analytical analyses and illustrate infection-transmission dynamics, we obtained numerical results from the two models using the following values from literature Hernandez-Vargas et al. 2014;Toapanta and Ross 2009): δ E = 0.02, p = 4.4, K V = 5.6234 × 10 5 , c V = 1.24 × 10 −6 , E 0 = 1 × 10 6 , r = 0.33 and K E = 2.7 × 10 3 , whereas S 0 = 1 × 10 3 , δ S = 4.74 × 10 −5 and δ I = 0.5 were chosen judiciously. The value of δ I is based on Brauer and Castillo-Chavez (2012), Communicable Disease Surveillance Centre (1978), where roughly half of the infected population is removed each day. We assumed that the average life expectancy of a healthy individual is 80 years, hence our choice of δ S .

Mathematical analysis of the multiscale model
The basic reproduction number (R 0 ) is a fundamental concept used to establish the cornerstone results for public health initiatives, such as the epidemic and herd immunity thresholds (Anderson and May 1992). The number is defined as the average number of secondary infections caused by a typical infected case in a completely susceptible population over its infectious period (Diekmann et al. 1990). For each of the subsystems, we derived the respective R 0 as a bifurcation parameter; the within-and between-host R 0 are denoted R W 0 and R B 0 , respectively.

Within-host subsystem
The basic reproduction number for the within-host subsystem (1) is given by With the parameter values in Sect. 2.4, R W 0 is approximately 3.5484. To justify our definition, we have for V ≈ 0 and E ≈ E 0 , i.e., at the outset of infection. Thus, the virus is cleared when R W 0 < 1 and persists when R W 0 > 1. According to the differential equation for E, each equilibrium point of (1) takes the form (V, E) where E > 0. Moreover, is the unique infection-free equilibrium for (1). Let x * be an equilibrium point of (1) whose coordinates are positive. To determine x * , we define a function G 0 in the following form Then G 0 has an inverse G −1 0 , given by Furthermore, G −1 0 (0) is a zero of G 0 , i.e., G 0 (G −1 0 (0)) = 0. The following lemma is a major step to obtaining x * .

Lemma 1 There exists a solution V of the following equation
if and only if R W 0 > 1. This solution is unique and lies in the open interval (0, G −1 0 (0)).
Proof According to the conditions in (2), the function G is increasing, from which On the other hand, so that G 0 is decreasing. We argue by defining H = G 0 − G, which is a decreasing function. Then each solution of (9) corresponds to a root of H on (0, K V ). Note that Consequently, H has a root on (0, K V ), and (9) admits a solution, only if R W 0 > 1. Now, suppose that R W 0 > 1, from which H (0) > 0 and 0 < G −1 0 (0) < K V . Observe that Hence, there exists a solution of (9) if and only if R W 0 > 1; this solution is uniquely given by v * ∈ (0, G −1 0 (0)).

Theorem 1 Consider the chronic equilibrium point x
Proof According to the equations of (1), the coordinates of x * satisfy the following: Equation (10) is equivalent to (11), from which 0 < V * < K V . On the other hand, we transform Eq. (12) into Applying Eqs. (6) and (10) to Eq. (13), we get (3), then we can express V * as a function f of R W 0 assuming r > δ E . This requires an equation P(V * , R W 0 ) = 0 from which we solve for V * . To this end, we introduce the following functions: Then Hence, we obtain the quadratic polynomial which has two distinct roots in μ. Let f be the smaller root of P(μ, x) in μ, given by Theorem 2 The function f is increasing (i.e., f > 0), and f (x) > 0 for all x > 1.

Proof
We compute the derivative f as Note that α 1 < 0 and D > 0. Moreover, we obtain The following theorem expresses V * as the function f of R W 0 .

Theorem 3 Assume the function G in the Michaelis-Menten form (3) and r
and note that P(μ) is concave down. Thus, and we must have Remark 1 In principle, given a function G that satisfies the conditions of (2), one can establish a function f that satisfies the properties in Theorem 2 and the equation . This is illustrated in Fig. 2 for three chosen examples of G, where f yields from the intersection of the curves y = G 0 (x) and y = G(x). Note in each example that the qualitative features of f , i.e. monotonicity and boundedness, inherit from G 0 .
Remark 2 Future research may prove the existence of f and its properties for any given G by computing H −1 (0) or implementing a global implicit function theorem.
We now establish the local stability of the within-host model (1).
Each pair of blue and black curves intersects at x = V * (Theorem 1). Moreover, each choice of G determines V * as a monotone bounded function of R W 0 (color figure online) Theorem 4 The following statements hold: Proof Consider the following Jacobian matrix associated with (1): Then J W (x 0 ) is a lower-triangular matrix. Moreover, recalling (2) and (6), we obtain the following eigenvalues of J W (x 0 ): This establishes statement 1. Now, we assume R W 0 > 1 and appeal to Theorem 1. Denote the (k, k)-entry of while Eqs. (7) and (9) Thus, J W (x * ) has a negative trace. On the other hand, G (V * ) > 0 by (2), from which J W (x * ) has a positive determinant. Accordingly, each eigenvalue of J W (x * ) has a negative real part. Therefore, x * is asymptotically stable. Statement 2 is proven.
Theorem 4 indicates a unique equilibrium point x for the system (1) that is either asymptotically stable or is nonhyperbolic; that is, We need the following hypothesis in our forthcoming analysis of the between-host subsystem (4): Thus, we say that x is globally asymptotically stable. A consequence for those solutions in hypothesis (H) is that This reiterates R W 0 as an infection threshold for the system (1): the virus is eradicated when R W 0 < 1, and persists when R W 0 > 1. The phase plane portraits in Fig. 3 supports the global stability of x = ( V , E) in hypothesis (H).

Coupled within-host and between-host model
Consider hypothesis (H), and assume that R W 0 = 1. Then, the system (1) equilibrates to  Fig. 3 Phase plane portraits for the within-host subsystem (1) at two different values of R W 0 . Each curve represents a solution (V (t), E(t)) for t ≥ 0, with initial data represented by the blue endpoint. Equilibrium points and the heteroclinic orbit connecting x 0 to x * are colored red. The following parameter values were used: c V = 1.24 × 10 −6 ; N E = 2 × 10 4 ; δ E = 2 × 10 −2 ; r = 3.3 × 10 −1 ; K E = 2.7 × 10 3 ; K V = 5.6234 × 10 5 ; and p = R W 0 c V N E /δ E (color figure online) according to (15). Thus, we may assume in (4) that V = V , yielding the following system: Let us call (18) the limiting case of the between-host subsystem. If R W 0 ≤ 1, then β( V ) = β(0) = 0 by (17) and the differential equations of (18) yield lim t→∞ (S(t), I (t)) = (S 0 , 0) regardless of initial condition and parameter values. For R W 0 > 1, the equations of (18) are evaluated at V = V * to yield the following system: From now on, we assume that R W 0 > 1 and analyze the limiting case (19). We realize the following between-host basic reproduction number from (19) : Indeed, we have for S ≈ S 0 and I ≈ 0, i.e., in a disease-free population. Thus, we expect the disease to be eradicated when R B 0 < 1, and an endemic to occur when R B 0 > 1. It follows from the differential equation of S that each equilibrium point of the system (19) is of the form (S, I ) where S > 0. Moreover, is the unique disease-free equilibrium point for (19).
Let y * = (S * , I * ) be an endemic equilibrium point of (19), where both S * and I * are positive. Then y * uniquely exists according to the following theorem.
Theorem 5 Assume that R W 0 > 1. Then y * = (S * , I * ) exists, with unique coordinates if and only if R B 0 > 1.
The local stability of the system (19) is given by the following result.
Proof We begin with the following Jacobian matrix associated with (19): Then J B (y 0 ) is an upper-triangular matrix, with eigenvalues −δ S and Thus, statement 1 holds.
To prove statement 2, we assume that R B 0 > 1 and denote the (k, k)-entry of J B (y * ) by J B k (y * ). Note that β(V * ) > 0 by (5), hence J B 1 (y * ) < 0. Meanwhile, by Theorem 5. Thus, the matrix J B (y * ) has a negative trace and a positive determinant, and each eigenvalue of J B (y * ) has a negative real part. Therefore, statement 2 holds.
We obtain from Theorem 6 a unique equilibrium point y = ( S, I ) of (19) that is either asymptotically stable or nonhyperbolic. That is, Denote by W S ( y) the intersection of X with the stable manifold of y. It is enough to establish to prove the theorem. We remark that W S ( y) ⊆ Y due to the forward invariance of the I = 0 plane.
We appeal to the LaSalle invariance principle (Hsu 2005;Wolkowicz and Lu 1992) to establish (25). For w = S, I , we define a function L w : Y → R by Consider the sum L = L S + L I and its time derivative L . Let Then we have after some computations. Observe that Substituting (27)   In addition to Theorem 7, we have lim t→∞ (S(t), I (t)) = (S 0 , 0) = y 0 (regardless of R B 0 ) whenever I (0) = 0. Figure 4 illustrates the phase plane portraits of the system (19), showing the global stability of y = ( S, I ) established in Theorem 7.
Remark 3 A reduction of (19) to a single equation is viable for δ S = δ I = δ, since we would then have lim t→∞ [S(t) + I (t)] = S 0 . However, the difference δ I − δ S is typically assumed a positive quantity to take into account excess mortality by infection (Gilchrist and Coombs 2006).
Remark 4 One may generalize the between-host model by taking δ I to be a bounded positive function on nonnegative real numbers. Then the stability results in Theorems 5, 6 and 7 remain true with δ I = δ I (V * ).

Relating the two basic reproduction numbers
It is useful to describe the dynamics of both systems (1) and (19) with respect to a single basic reproduction number. Assuming that R W 0 > 1, we have by Eqs. (6) and (20). With the equation V * = f (R W 0 ) from Theorem 3, we establish below that R B 0 is an increasing function of R W 0 . Given that β is defined for nonnegative values, we must remark that β • f is defined on the interval [1, ∞) where f is nonnegative.

Theorem 8 Assuming that r > δ E , define a function F on [1, ∞) by
where f is given by Eq. (14).
Proof For x > 1, we have Observe from (5)   In reality, the basic reproduction numbers may change as the disease progresses, from which F becomes a function of time. Considering this possibility, we would like to establish a subhomogenity result.
Given two subhomogeneous functions ϕ and φ with degrees j and k, respectively, it can be shown that the degree of subhomogenity is j for the scalar multiple uϕ (u > 0), j + k for the product ϕ · φ, and jk for the composite ϕ • φ.
We note from Eq. (29) that F is defined on the interval [1, ∞) where f ≥ 0. Thus, we compose F with the following bijective function so that F := F • T is defined on nonnegative real numbers. Recalling Eq. (14), we also denote f := f • T . Proof For all λ ≥ 1 and x > 0, the subhomogenities of β and f yield due to the increasing property of β assumed in (5). Thus, β • f is subhomogeneous of degree jk. Since T is subhomogeneous of degree one and it follows that F is subhomogeneous of degree jk + 1.

Discussion
Multiscale models have shown to be capable of generating hypotheses and adding new insights into infectious disease mechanisms (Handel and Rohani 2015). Linking the two scales has led to the identification of missing pieces in infectious disease research and is the only way to capture the effects of feedback between the scales (Mideo et al. 2008). While there exist many obstacles (Gutierrez et al. 2015;Handel and Rohani 2015;Mideo et al. 2008), multiscale modeling is a promising approach that would improve our understandings of the epidemic control and prevention.
One of the significant bottlenecks for multiscale models of infectious diseases is the lack of appropriate data to quantify the functional link between within-host viral kinetics and between-host transmissibility (Handel and Rohani 2015), leading to results which are difficult to validate and translate into practical policies (Hellriegel 2001).
To close this gap, more tailored experiments are needed to characterize the biological relationship between the viral load and its transmission probability. Before these can be realized, detailed analytical and numerical analyses targeted at the bridge between within-and between-host systems are promptly required to provide initial understandings and help design testable functional hypotheses.
In this paper, we analyzed a multiscale model that comprises a within-host viral infection model nested in a SI transmission model. A rigorous stability analysis and numerical simulations for the multiscale model were performed. We derived stability analyses considering that the stimulation of the immune system by the virus (V ) is represented by any function G(V ) with the properties G(0) = 0 and G (V ) > 0.
The unique non-trivial equilibrium point for the within-host system (1) is established in Theorem 1 and associated with the unique intersection of two curves (Lemma 1). It is known that a linear response of the virus and a bounded functional response of the effector T-cells do not yield periodic solutions for the corresponding model. Figure 3 confirms not only the absence of closed trajectories but also the global stability of the within-host system.
With a specific expression for G(V ), we obtained an expression V * = f (R W 0 ), which consequently provides a relationship between the basic reproduction numbers.
The existence of f for a general G is an open problem, although our numerical examples verify it (Remarks 1 and 2). Meanwhile, the global stability of the betweenhost system (Theorem 7) was established through the construction of a Lyapunov function. We derived two basic reproduction numbers, for the within-host level and between-host, respectively. Viral infections with R W 0 below unity would indicate virion clearance within a host; otherwise, a chronic state could be established (Fig. 3). For viruses with a very high replication rate, such as influenza (Hatta et al. 2010) and Ebola , the virus replication likely outpaces the immune response (Hatta et al. 2010;Nguyen and Hernandez-Vargas 2017). Early interventions to stem down viral replication rate could promote survival and recovery from infection (Hatta et al. 2010;Nguyen and Hernandez-Vargas 2017). Notable examples include the use of neuraminidase inhibitors to contain influenza infection (Ferguson et al. 2005;Longini et al. 2005) and the use of antiretroviral treatment (ART) as a public measure to prevent HIV (Heesterbeek et al. 2015). Given the same disease and similar intervention approaches in place, i.e., fixed values for δ I and β(V * ), the epidemic size then depends on the influx of new susceptible individuals [Eq. (21)]. Approaches, such as childhood vaccinations, diverge the influxes into non-susceptible classes and thus could lead to disease eradication. As per the SI model's assumption of a constant population, Eq. (21) also shows that the epidemic size will saturate regardless of R B 0 . As R B 0 will also saturate when R W 0 crosses a certain value (Fig. 5), it follows that epidemic interventions aimed at controlling viral replication would need to be highly effective to reduce the epidemic size. In other words, if an intervention can largely reduce the R W 0 but this value is still above a certain threshold then the intervention would have a negligible impact on the epidemic size. For our set of parameter values (Sect. 2.4), R W 0 ≈ 1.5. Figure 6 shows the dynamics of the viral load and the number of infected cases for different values of R W 0 .
The above results could still hold with a different coupling function between the two scales was varied (Fig. 5). The correct link function for many diseases is still debatable because empirical data are lacking (Gutierrez et al. 2015;Handel and Rohani 2015;Mideo et al. 2008). For this reason, we tested three functional forms expressing the transmission rate based on viral load including a linear, logistic, and Michaelis-Menten function as discussed elsewhere (Handel and Rohani 2015). Figure 5 shows that while the effect of R W 0 on R B 0 is qualitatively indifferent among the functions, there is a marked difference quantitatively. Linear coupling led to the transmission increased roughly triple compared to the others. Thus, coupling with a linear function for the pathogen may lead to an overestimate in disease transmission. This is illustrated in Fig. 7, where the trajectory corresponding to linear transmission yielded a much faster and stronger epidemic trajectory. For a linear or Michaelis-Menten (saturation) form of G and our examples of β, the relationship between the two reproduction numbers maintained the same qualitative properties (Fig. 5). This can be seen as a consequence of the asymptotic nature of the function G 0 (V ) for increasing values of R W 0 (Fig. 2). Our theoretical results hold with abstract conditions for both G and β, as given by (2) and (5). However, concrete formulations are typically endowed with biological properties like the functional response curves of predator-prey theory (Holling 1959). Nevertheless, the corresponding biological ramifications of our stability analysis apply to a wide range of disease scenarios compatible with our model. While our disease model is tractable enough to yield theoretical results, future works may employ our analytical framework to more general disease problems. Such generalizations can account, for example, logistic growth and carrying capacity from predator-prey dynamics. A particular extension of our work arises in disease-induced mortality, where δ I is a function of the virulence rate, denoted α. If we define δ I by then the stability results of the between-host subsystem (Theorems 5, 6 and 7) remain true where δ I = δ S + α(V * ); see Remark 4. However, the relationship between the reproduction numbers will now depend on both the transmission and virulence rates, whose trade-off connection was revealed in Gilchrist and Coombs (2006). The dependence of the curve R B 0 = F(R W 0 ) on both transmission and virulence is demonstrated in Fig. 8, where F loses its monotonicity in favor of a local maximum near R W 0 = 1 for K α > 100. Further investigation can be done by comparing the dynamics and the variation of reproduction numbers with different forms of virulence.
Another extension of this work may include time-varying parameters centered around vaccination and disease control (Liu and Stechlinski 2012;Ma and Ma 2006;De la Sen et al. 2010). To this end, Proposition 1 serves as an initial step and lends itself to models with subhomogeneous nonlinearities, including classes of monotone, periodic, almost periodic and cooperative systems (Bokharaie et al. 2011;Zhao 2003;Zhao et al. 2017).
In Proposition 1, we established that the subhomogenity of F = F • T , where T is the bijective function in (31), is inherited from the transmission rate β and f = f • T . Our examples of β (Sect. 2.3) are subhomogeneous, suggesting that subhomogenity is a biological property of virion-induced transmission. Further studies should investigate the subhomogenity of f . The analysis of infectious disease dynamics from within-to between-host requires the coupling of systems that are usually studied in isolation. The high level of detail involved in this approach entails considerable effort. As of now, a complete framework to study multiscale infectious diseases has not prevailed, but it is a constant target. This work contributes basic understandings of the within-and between-host models in consideration and casts light on potential effects of the coupling function on linking the two scales. Future experimental studies are needed to verify our results to improve multiscale models.