1 Introduction

Sensitivity Analysis (SA) is a common methodological approach for determining the expected impact of control strategies on common disease outbreak quantities such as final size and basic reproduction number Mpeshe et al. (2011). SA has been used extensively for ODE models but has not been extended to age-structured models. In this study, we develop a method that extends SA to immuno-epidemiological models (Tuncer et al. 2016; Gulbudak and Browne 2020; Gandolfi et al. 2015; Angulo et al. 2013). We couple a time-since-infection-structured epidemiological system with an ODE immunological model. The main motivation for this modeling approach is that many recent vector-borne epidemic models share the limitation of exploring only between-host transmission while ignoring the impact of within-host virus-with-immune-response interactions, which may be important to guide drug and vaccine development, for example. Here we demonstrate the approach through a multi-scale model, first introduced in Gulbudak et al. (2017). The parameters for this model are studied in Tuncer et al. (2016) based on Rift Valley Fever Disease (RVFD) immunological data and human epidemiological data from the \(2006-2007\) Kenya Outbreak ([6], Munyua et al. 2010). Using the SA approach developed here, we quantify the impact of within-host parameters on the Rift Valley Fever Disease (RVFD) dynamics. The same multi-scale modeling framework can be adapted to other arbovirus diseases such as Dengue and West Nile Virus (WNV).

RFVD is a viral disease transmitted by mosquitoes, mainly from the Aedes and Culex genera, and causes illness and death in several different mammal species, including livestock (e.g. cattle, buffalo, sheep, goat, and camel), as well as in humans. RVFD has resulted in significant negative socio-economic impacts, for example, due to abortion among RVF-infected livestock and high mortality among younger ones. In 2018, a panel of experts convened by the World Health Organization (WHO) listed RVF among diseases that pose big public-health risks, yet few or no intervention strategies have been developed. SA can provide helpful insights on the impact of possible pharmaceutical interventions for RVFD control.

Sensitivity analysis has been utilized in several ODE models to assess the impact of epidemic parameters on epidemic quantities. For example, Gaff et al. (2007, 2011) considered an ODE vector-host RVF model to assess the effectiveness of some control interventions on RVFD. Fischer et al. (2013) utilized SA to investigate the effect of temperate climate on the RVFD dynamics. Mpeshe et al. (2011) formulated an ODE model of RVF incorporating parameters dependent of human behavior to investigate disease dynamics and explore sensitivity of the model to variation in those parameters. Xiao et al. (2015) recently studied the effect of both seasonality and socioeconomic status in a multi-patch model. To the best of our knowledge, the SA of immuno-epidemiological models has never been carried out, despite the value and usefulness of SA of the underlying immunological model parameters on epidemic variables related quantities related to them.

In this study, we develop a novel approach for SA in immuno-epidemiological models to investigate the impact of immunological parameters on the disease dynamics. In particular, we consider a time-since-infection-structured vector-host model in which epidemiological model parameters are described as functions of within-host virus-antibody densities that are governed by an ODE system. We first define the basic reproduction number, \(\mathcal {R}_0\) that serves as a threshold between extinction and persistence of RVFD. Then we use this SA approach to investigate the impact of changes in immunological parameters on epidemic quantities such as basic reproduction number, and final disease abundance when \(\mathcal R_0>1.\) Interestingly, our analytical and numerical results suggest that immunological parameters such as viral growth rate and immune activation rate can have a large impact on disease outcomes, underscoring the importance of pharmacological intervention strategies.

This paper is organized as follows. In Sect. 2, we present an immuno-epidemiological model, first introduced in Gulbudak et al. (2017), Tuncer et al. (2016), and summarize the stability and persistence conditions for the disease. In Sect. 3, we develop a novel approach for SA in immuno-epidemiological models to assess the impact of the within-host parameters on epidemic quantities. Furthermore in Sect. 3.2, we consider three distinct stages of an outbreak—initial, peak and die-out—and show how infectiousness of hosts at these different stages of infection is altered by slight changes in the immunological parameters through the phases of an outbreak. In the last section, we summarize our results and draw some conclusions.

2 An Immuno-Epidemiological Vector-Borne Disease Model

A simple model that captures vector-borne disease spread on two-scales (namely immunological and the epidemiological) was first introduced in Gulbudak et al. (2017). At the individual scale, we consider the following immune response model with three state variables representing serum density of the pathogen and of two specific antibodies released by B-cell lymphocytes —IgM and IgG— with densities given, respectively, by \(P=P(\tau ),\) \(M=M(\tau )\) and \(G=G(\tau ),\) where \(\tau \) is the time elapsed since infection:

$$\begin{aligned} \frac{dP}{d \tau } =\left( f(P)-\theta M-\delta G\right) P,\ \frac{dM}{d \tau } =\left( aP-(q+c)\right) M,\ \frac{dG}{d \tau } =qM + bGP.\ \end{aligned}$$
(1)

The initial values are \(P(0)=P_0>0, M(0)=M_0\ge 0,G(0)=G_0\ge 0\), and \(M_0+G_0>0\) to ensure that there is a pathogen and an immune response. It is assumed that the pathogen replicates with a logistic per capita growth rate \(f(P)=r\left( 1-\frac{P}{K}\right) \) and that, upon exposure to the virus, the IgM and IgG antibodies get activated at unit rates a and b, respectively. The IgM antibodies are responsible for a a quick immunological response (innate): they kill virus at a unit rate \(\theta \) and decay at a unit rate c. The IgG antibodies kill the pathogen at a unit rate \(\delta \), and they are mainly responsible for long-term immunity (adaptive) (Gulbudak et al. 2017; Tuncer et al. 2016). Mature B-cells activated by antigen stimulation proliferate quickly in lymphoid follicles and undergo genetic alterations resulting in a switch of the immunoglobulin isotype produced from IgM to either IgG, IgE, or IgA. To keep the immune response model lower-dimensional, we model the B-cell population indirectly through the IgM and IgG antibodies it produces. We incorporate this switch in antibody production by the B-cells by “converting” IgM antibodies to IgG antibodies at a unit rate q (Honjo et al. 2002). All parameters and state variables of this within-host model and their definitions are given in Tables 1 and 2. Different dynamics of this immune response model can be found in Fig. 1.

A detailed analysis of the immunological model extended to include vector-to-host inoculum size dependent on age-of-infection of vectors was presented in Gulbudak (2020). In general upon viral progression within an infected mosquito midgut, the amount of pathogen in an infected mosquito’s saliva dynamically changes with respect to its infection-age, determining the vector-to-host inoculum size \(P_0\) that is the amount of pathogen injected to a susceptible host by an infectious vector during the infectious contact. Assuming that \(P_0\) is constant, all analytical results therein also hold for our present immunological model. In particular, the following result is established Gulbudak (2020):

Theorem 1

If \(M_0>0\) (or \(G_0>0\)), then the pathogen eventually clears (\(\displaystyle \lim _{\tau \rightarrow \infty }P(\tau )=0\)), the IgM antibodies decay to zero after viral clearance and, subsequently, the IgG memory antibodies reach a steady-state; i.e. \(\displaystyle \lim _{\tau \rightarrow \infty }M(\tau )=0 \ \text {and} \ \displaystyle \lim _{\tau \rightarrow \infty }G(\tau )=G^+,\) where \(G^+>0\) depends on the initial condition \((P_0,M_0,G_0)\).

Table 1 Model parameters: parameter estimates of the immuno-epidemiological model (1)-(3) fitted Tuncer et al. (2016)
Full size table
Table 2 State variables
Full size table

The innate and adaptive immune responses vary from organism to organism and depend on several factors, including the initial viral dose, the strain of the virus, the age of the host, and the species of the host. However, infections like RVFV in the laboratory generally follow three qualitative scenarios:

  • The viral load within the host grows rapidly to high levels within the host, and the host dies.

  • There is a delayed onset of complications of infections: e.g. viral load appears to be controlled after an initial infection but makes a resurgence after disseminating into new tissue, such as the central nervous system. This often leads to long-term consequences such as blindness and may also be lethal.

  • The viral load within the host is brought under control and annihilated by a robust immune response.

Our immunological model captures the general dynamics of these different infection and immune-response scenarios as seen in Fig. 1 and as suggested by the analytical result above (Theorem 1). In general, within-host dynamics in arbovirus diseases behave as depicted in Fig. 1 (left), corresponding to the third scenario mentioned above. Also, in a lab experiment with monkeys infected with RVFV, it was observed that monkeys died when the viral concentration rebounded as seen in Fig. 1(right). In surviving monkeys, the within-host dynamics resembled the third scenario Morrill et al. (1989). In an ideal immune response, the immune system releases virus-specific antibodies, targeting the virus and bringing it to extinction. These dynamics are captured in Fig. 1 (left). However, sometimes there is a resurgence of virus within the host as the virus spreads into new organs and tissues, such as the central nervous system. This is often associated with delayed onset of symptoms, including blindness, and may result in fatalities Pepin et al. (2010). Our model can capture these dynamics as well, as can be seen in Fig. 1 (right).

Fig. 1
figure 1

Different within-host pathogen-immune response antibody dynamics Left. General within-host dynamics. Right. Resurgence of virus. The experiments done on monkeys suggest that in some cases, viral load appears to be controlled after an initial infection but makes a resurgence after disseminating into new tissue Pepin et al. (2010) (Color figure online)

Full size image

At the population scale, we consider a time-since-infection structured vector-host model, (2), where \(S_H(t)\) and \(R_H(t)\) are the numbers of susceptible and recovered individuals in the host population and \(i_H(\tau ,t)\) is the infected-host density, structured by age-of-infection \(\tau \). The total number of infected individuals is \( I_H(t) = \int _0^\infty i_H(\tau ,t) d\tau \). The vector compartments \(S_V(t)\) and \(I_V(t)\) represent, respectively, the size of susceptible and infected vector populations at time t. The full model is as follows:

$$\begin{aligned} \left\{ \!\! \begin{aligned} \frac{d S_H}{dt}&= \Lambda - \beta _V S_H(t)I_V(t) -d S_H(t),\\ \frac{\partial i_H}{\partial t}+\displaystyle \frac{\partial i_H}{\partial \tau }&= -(\alpha (\tau )+ \kappa (\tau )+ \gamma (\tau )+ d)\; i_H(\tau ,t),\quad i_H(0,t) = \beta _V S_H(t) I_V(t),\\ \frac{d R_H}{dt}&= \displaystyle \displaystyle \int _0^\infty {\gamma (\tau )i_H(\tau ,t)d\tau }-dR_H(t),\\ \frac{d S_V}{dt}&= \eta - \left( \int _0^\infty {\beta _{H}(\tau )i_H(\tau ,t) d \tau }+\mu \right) S_V(t),\\ \frac{d I_V}{dt}&= S_V(t) \int _0^\infty {\beta _{H}(\tau )i_H(\tau ,t) d \tau }-\mu I_V(t). \end{aligned} \right. \end{aligned}$$
(2)

To bridge the scales from individual to population, the unit infectiousness and the disease-induced death and recovery rates are formulated as functions of immunological variables as suggested by data in Handel and Rohani (2015), Fraser et al. (2014), as follows:

$$\begin{aligned} \left\{ \begin{aligned} \beta _H(\tau )=\frac{C_\beta P(\tau )}{C_0+ P(\tau )},\\ \alpha (\tau )= \zeta P(\tau ),\\ \kappa (\tau )=\xi M(\tau ), \\ \gamma (\tau )=C_\gamma \frac{G(\tau )}{P(\tau )+\epsilon _0}. \end{aligned}\right. \end{aligned}$$
(3)

Figure 6 in the Appendix shows how some of the corresponding epidemiological parameters evolve over the course of host infection. Note that Tuncer et al. (2016) considers a mass action term, satisfying the balance equation for describing the interactions between vectors and hosts with following underlying assumptions:

  • \(c_v N_v=\)total number of contacts (bites received) per host. Then total number of contacts (bites received) that hosts’ have: \((c_v N_v) \times (N_H)\).

  • \(c_H N_H=\) total number of contacts (given bites) per vector. Then total number of contacts (given bites) that vectors’ have: \((c_H N_H)\times (N_v)\). Therefore it makes sense when the host population size is small (or not changing much).

  • To satisfy balance equation (conservation law), we assume that \(c_v=c_H\).

  • Here c is proportionality constant. Mass action assumes that per vector (or per host) contact is proportional to total host (or total vector) population. Therefore it has different meaning than the one in standard incidence.

Changes to the functional form of the interaction terms such as considering frequency-depending force of infection rate (standard incidence) will be further explored in later work.

A schematic diagram for the full system is presented next in Fig. 2.

Fig. 2
figure 2

Schematic illustration of the immuno-epidemiological vector-host model(1)–(3) (Color figure online)

Full size image

Next, we highlight the threshold dynamics of the system, including long-term behavior of the solutions, explicit expressions of equilibria and basic reproduction number.

Let us define the immune-response-dependent reproduction number of the epidemic as

$$\begin{aligned} \mathcal R_0 = \displaystyle \frac{\beta _V S^0_H }{\mu }\; \displaystyle \int _0^\infty {S_V^0 \beta _H(\tau ) e^{-\displaystyle \int _0^\tau {(\alpha (s)+\kappa (s)+\gamma (s)+d)ds}}d\tau }, \end{aligned}$$
(4)

where \(S^0_H =\displaystyle \frac{\Lambda }{d}\) and \(S_V^0=\displaystyle \frac{\eta }{\mu }.\) A detailed analysis of our model (1)–(3) including also vector age-of-infection was presented in Gulbudak (2020). The vector compartments here correspond to the ones in that paper integrated over vector-infection-age and, therefore, all analytical results therein (e.g. threshold conditions) also hold for our present model. In particular, the following results are established for our model:

Theorem 2

The disease-free equilibrium (DFE) \(\mathcal E_0=(S^0_H,0,0,S_V^0,0)\) is locally asymptotically stable if \(\mathcal {R}_0<1\) and unstable if \(\mathcal {R}_0>1.\)

This result is actually global for \(\mathcal E_0\):

Theorem 3

\(\mathcal E_0\) is globally asymptotically stable if \(\mathcal {R}_0<1\) and unstable if \(\mathcal {R}_0>1.\)

Furthermore, when \(\mathcal {R}_0>1\), the system (2) has a unique endemic equilibrium (EE) \(\mathcal E^* =(S^*_H,\ i^*_H(\tau ),\ R^*_H,\ S^*_V, \ I^*_V)\) (also presented in Gulbudak et al. (2017)), where

$$\begin{aligned} \left\{ \begin{aligned} S^*_H&= \left( \beta _V \displaystyle \frac{S^*_V}{\mu }\displaystyle \int _0^\infty {\beta _H(\tau )\pi _H(\tau )d\tau }\right) ^{-1},\\ i^*_H(\tau )&= i^*_H(0)\pi _H(\tau ), \\ R^*_H&=\frac{i^*_H(0)}{d} \displaystyle \int _0^\infty {\gamma (\tau )\pi _H(\tau )d \tau }, \end{aligned}\right. \end{aligned}$$
(5)
$$\begin{aligned} \left\{ \begin{aligned} S^*_V&=\eta \left( i^*_H(0)\displaystyle \int _0^\infty {\beta _H(\tau )\pi _H(\tau )d\tau }+\mu \right) ,\\ I^*_V&=\frac{S^*_V}{\mu }\displaystyle \int _0^\infty {\beta _H(\tau )i^*_H(\tau )d\tau }, \end{aligned} \right. \end{aligned}$$
(6)

with

$$\begin{aligned} i^*_H(0)= & {} S_0^H\left( 1-\displaystyle \frac{1}{\mathcal R_0} \right) / \left( \displaystyle \frac{1}{d} + \displaystyle \frac{\mu }{\beta _V\eta }\right) ,\ \pi _H(\tau ) \nonumber \\= & {} e^{-\displaystyle \int _0^\tau {(\alpha (s)+\kappa (s)+\gamma (s)+d)ds}}, \end{aligned}$$
(7)

and the epidemiological parameters \(\beta (\tau ), \ \gamma (\tau ),\ \alpha (\tau ),\) and \(\kappa (\tau )\) are given by (3). For this equilibrium we have the following results.

Theorem 4

The endemic equilibrium \(\mathcal E^* =(S^*_H, i^*_H(\tau ),R^*_H, S^*_V, I^*_V)\) is locally asymptotically stable whenever it exists (i.e. for \(\mathcal {R}_0>1\)).

Theorem 5

If \(\mathcal {R}_0>1,\) then the disease is uniformly weakly endemic.

In fact, the global stability of a unique endemic equilibrium is shown in Magal and McCluskey (2013) under an equivalent threshold condition to \(\mathcal R_0>1.\) These results show that \(\mathcal R_0\) represents both a threshold between extinction and persistence and a strict threshold for disease eradication.

3 Impact of Immune parameters on Epidemic Quantities

SA of immuno-epidemiological models not only provides a measure of the influence of epidemic parameters on the spread and abundance of the disease at the population scale, but also of those at within-host scale. In particular, SA of immunological parameters on epidemic quantities might provide valuable information on the expected impact of control strategies including those targeted toward:(i) the host population such as vaccination, and drug treatment, (ii) the vector population, such as Wolbachia-based control strategies, and (iii) curbing viral production by providing important comparisons of the efficacy of different immune variables. To determine the impact of immunological parameters on the disease dynamics at population scale, we carry out the SA of the epidemiological quantities \(\mathcal {R}_0, I_H^*\) and \(I^*_V\) with respect to the immune model parameters in (1). Notice that the prior studies on SA only focus on the impact of epidemic parameters on epidemic quantities.

The normalized forward sensitivity index of a quantity of interest (QOI) q to a parameter of interest (POI) p can be defined as the ratio of relative change in the variable to the relative change in the parameter:

$$\begin{aligned} \gamma ^ {q} _{p}=\frac{\partial q}{\partial p}\times \frac{p}{q}. \end{aligned}$$

The baseline parameter values and ranges for RVFD are presented in Table 1. The baseline values were fitted to multi-scale data in Tuncer et al. (2016): for immunological parameter estimations, immunological RVF time-series data was obtained from livestock (in laboratory experiments), and for the epidemiological model, incidence data was acquired from the \(2006--2007\) Kenya Outbreak ([6], Munyua et al. 2010). This simultaneous fitting of immunological and epidemiological model parameters induces practical identifiability of model parameters (Tuncer et al. 2016).

We first consider the normalized sensitivity index for the basic reproduction number \(\gamma _p^{\mathcal {R}_0}\). Note that for consistency regarding parameter estimates in Tuncer et al. (2016), we assume that \(\Lambda =d,\) and \(\eta =\mu ,\) which set the total number of susceptible hosts, \(S_H^0,\) equal to 1 at the DFE, \(\mathcal E_0\).

From (4), we have

$$\begin{aligned} \frac{\partial {\mathcal {R}_0}}{\partial p} = \frac{\beta _V}{\mu } \int _0^\infty \left( \frac{\partial \beta _H(\tau ,p)}{\partial p}\pi _H(\tau ,p)+\beta _H(\tau ,p)\frac{\partial \pi _H(\tau ,p)}{\partial p}\right) \,d\tau , \end{aligned}$$
(8)

where \(\beta _H(\tau ,p)\) and \(\pi _H(\tau ,p)\) are defined in (3) and (7), respectively, and

$$\begin{aligned} \frac{\partial \beta _H(\tau , p)}{\partial p}&=\frac{C_{\beta } C_0}{(C_0 + P(\tau ))^2} \frac{\partial P(\tau ,p)}{\partial p},\nonumber \\ \frac{\partial \pi _H(\tau , p)}{\partial p}&=-\pi _H(\tau )\int _0^\tau \left( \zeta \frac{\partial P(s,p)}{\partial p}+\xi \frac{\partial M(s,p)}{\partial p}\right. \nonumber \\&\left. \quad + C_{\gamma }\frac{\frac{\partial G(s,p)}{\partial p}(P(s,p)+\epsilon _0)- G(s,p)\frac{\partial P(s,p)}{\partial p}}{(P(s,p)+\epsilon _0)^2}\right) \,ds. \end{aligned}$$
(9)

Note that we need to compute \(\displaystyle \frac{\partial \varvec{x}}{\partial p}(\tau ,p)\), where \(\varvec{x}(\tau ,p) \!=\! (P(\tau , p), M(\tau , p), G(\tau , p))\), the state variables of the within-host model (1). From Clairaut’s Theorem, we have

$$\begin{aligned} \frac{\partial }{\partial \tau }\left( \frac{\partial \varvec{x}}{\partial p}\right) =\frac{\partial }{\partial p}\left( \frac{\partial \varvec{x}}{\partial \tau }\right) , \end{aligned}$$
(10)

where \(\partial \varvec{x}/\partial \tau \) is the right-hand side of the system (1). Let \(\varvec{f}(\varvec{x}, p) = \partial \varvec{x}/\partial \tau \), and consider the following system for the derivatives \(\partial \varvec{x}/\partial p\):

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial }{\partial \tau } \left( \frac{\partial \varvec{x} }{\partial p}\right) = \frac{\partial \varvec{f} }{\partial \varvec{x}} \frac{\partial \varvec{x} }{\partial p} + \frac{\partial \varvec{f} }{\partial p},\\&\frac{\partial \varvec{x} }{\partial p} (0) = \varvec{0}, \\ \end{aligned} \right. \end{aligned}$$
(11)

where \(\partial \varvec{f} /\partial \varvec{x}\) is the Jacobian matrix of the right-hand side of the system (1). For example, taking parameter \( p = a\), the IgM immune response activation rate, we can derive the corresponding system (11) as follows

$$\begin{aligned} \underbrace{ \begin{pmatrix} \frac{\partial }{\partial t} \left( \frac{\partial P }{\partial a}\right) \\ \frac{\partial }{\partial t} \left( \frac{\partial M }{\partial a}\right) \\ \frac{\partial }{\partial t} \left( \frac{\partial G }{\partial a}\right) \\ \end{pmatrix}}_{\frac{\partial }{\partial t} \left( \frac{\partial \varvec{x} }{\partial a}\right) } =\underbrace{ \begin{pmatrix} r\left( 1- \frac{2P}{K}\right) -\theta M-\delta G &{} \ -\theta P &{}\ -\delta P \\ a M &{} \ a P- (q+c) &{}\ 0 \\ b G &{} \ q &{}\ b P \end{pmatrix}}_{{\frac{\partial {\varvec{f}}}{\partial {\varvec{x}}}}} \underbrace{ \begin{pmatrix} \frac{{\partial P}}{{\partial a}}\\ \frac{{\partial M}}{{\partial a}}\\ \frac{{\partial G}}{{\partial a}}\\ \end{pmatrix}}_{{\frac{\partial {\varvec{x}}}{\partial \mathrm{a}}}} + \underbrace{ \begin{pmatrix} 0\\ M P\\ 0\\ \end{pmatrix}}_{{\frac{\partial {\varvec{f}}}{{\partial \mathrm{a}}}}}. \end{aligned}$$

By solving the extended system (11)–(1), we obtain the solutions for \(\varvec{x}\) and the derivatives \(\partial \varvec{x}/\partial p\) at all times \(\tau \), which are then utilized to estimate the integrands in (8) and (9). We then numerically approximate the integrals by using the trapezoidal rule. Note that (11) is the standard derivation for the variational system of an ODE.

Similarly, we assess the impact of within-host parameters on the final epidemic size, consider the endemic host population size \(I^*_H\) and endemic vector population size \(I_V^*\) as the QOIs, which are given in (5) and (6), respectively. The corresponding sensitivity matrices are derived as follows:

$$\begin{aligned}&\frac{\partial I^*_H}{\partial \varvec{p}}=\displaystyle \int _0^\infty \!\bigg [\frac{\partial i^*_H(0)}{\partial \varvec{p}}\pi _H(\tau )+ i^*_H(0)\frac{\partial \pi _H(\tau )}{\partial \varvec{p}}\bigg ]d\tau ,\ \ \ \ \frac{\partial i^*_H(0)}{\partial \varvec{p}}=\frac{1}{\mathcal {R}_0^2}\frac{\partial \mathcal {R}_0}{\partial \varvec{p}}\bigg / \bigg (\frac{1}{d} + \frac{1}{\beta _V}\bigg ),\\&\frac{\partial I^*_V}{\partial \varvec{p}}=\frac{1}{(\mathcal {R}_0 S^*_H)^2}\bigg (\frac{\partial \mathcal {R}_0}{\partial \varvec{p}} S^*_H -\frac{\mathcal {R}_0}{d}\frac{\partial i^*_H(0)}{\partial \varvec{p}}\bigg ). \end{aligned}$$
Fig. 3
figure 3

Impact of immune parameters on the epidemic threshold, \(\mathcal {R}_0.\) The immune-response-dependent epidemiological reproduction number, \(\mathcal {R}_0\) (left y-axis, blue dotted curve) along with the response curves for the extended SA, \(\gamma ^{\mathcal {R}_0}_{\varvec{p}}\) (right y-axis, orange curve) are plotted against the immune parameter values p over a range of values around the baseline estimates (Color figure online)

Full size image
Fig. 4
figure 4

Impact of immune parameters on the final epidemic size, \(\mathcal I_H^*.\) The final infected host disease abundance, \(\mathcal I^*_H,\) along with the response curves for the extended SA, \(\gamma ^{\mathcal {R}_0}_{\varvec{p}},\) are plotted against the immune parameter values p over a range of values around the baseline estimates (Color figure online)

Full size image

In Fig.3, the immune-response-dependent epidemiological reproduction number, \(\mathcal {R}_0,\) along with the response curves for the extended SA, \(\gamma ^{\mathcal {R}_0}_{\varvec{p}},\) are plotted against the immune parameter values p over a range of values around the baseline estimates. In Fig.3a, the epidemiological threshold \(\mathcal {R}_0\) is a non-monotone function of the within-host pathogen replication rate r and there is a critical pathogen replication rate \(r_c\) at which the reproduction number is maximal. We present the resulting sensitivity indices at the baseline parameters in Table 5, where the indices are sorted by magnitude. The numerical results suggest that:

  1. i)

    At the baseline (estimated) values of within-host immune parameters, the within-host viral growth rate r has the largest impact on \(\mathcal {R}_0:\) \(1\%\) increase in r causes \(0.6\%\) reduction in \(\mathcal {R}_0\), \(0.8\%\) reduction in final infected host abundance \(\mathcal I_H^*\) and \(0.72\%\) reduction in final infected vector abundance \(\mathcal I_V^*\).

In Fig.3, the reduction in \(\mathcal {R}_0\) with increasing within-host pathogen growth rate r is counter-intuitive, but well known in the literature Anderson and May (1982). The trade-off between pathogen virulence and infectiousness is first shown in Gilchrist and Sasaki (2002), articulated through a within-host model for directly transmitted diseases. Later it was extended to vector-borne diseases in Gulbudak et al. (2017). In particular, the reduction in \(\mathcal {R}_0\) for large r occurs due to three mechanisms: (i) death due to “aggressive” immune response, (ii) death due to pathogen exploitation of target cells, and (iii) decreasing infectious period due to large viral clearance. As the pathogen growth increases rapidly, the immune response becomes more robust (See Fig. 7). Both cases can lead to rapid death of the host, shortening the infectious period, and subsequently leading to a reduction in the average number of secondary cases caused by these infectious individuals. In addition, note that a \(1\%\) increase in within-host pathogen growth rate r in the parameter region \(r \in [0.1, 2]\) leads to up to \(8\%\) increase in \(\mathcal {R}_0\) (Fig.3), \(1\%\) increase in \(\mathcal I_H^*\) (Fig.4), and \(1.5 \%\) increase in \(\mathcal I_V^*\) (Fig.9), which are significant increases in population-scale disease quantities, suggesting that control strategies should be targeted toward reducing within-host pathogen growth rate to significantly change the disease outcomes in the long term. For example, within-host pathogen growth rate r can be reduced by immunizing the host population, or through drug treatment during outbreaks, hampering viral growth within-hosts.

  1. ii)

    The within-host pathogen carrying capacity K has the second largest impact on \(\mathcal {R}_0\) at the baseline parameter value \(K=7.57\times 10^{7}\): \(1\%\) increase in K causes \(0.079\%\) reduction in \(\mathcal {R}_0\), \(0.060\%\) reduction in \(\mathcal I_H^*\), and \(0.095\%\) reduction in \(\mathcal I_V^*\).

  2. iii)

    The within-host IgM immune activation rate a has the third largest impact on \(\mathcal {R}_0\) at the baseline parameter value \(a=1.1\times 10^{-7}\): \(1\%\) increase in a causes \(0.019\%\) reduction in \(\mathcal {R}_0\), \(0.015 \%\) reduction in \(\mathcal I_H^*,\) and \(0.023 \%\) in \(\mathcal I_V^*.\)

The faster the IgM antibodies activate, the faster they clear the pathogen, causing reduction in disease transmission (See Fig.3c). A drastic decrease in \(\mathcal {R}_0\) for large values of a suggests that a large improvement in vaccine efficacy is needed for efficient reduction of population-scale transmission.

  1. iv)

    The immune response switching rate q has the fourth largest impact on the epidemic quantities, albeit almost negligible: a \(1\%\) increase in q causes \(1.6 \times 10^{-3}\%\) increase in \(\mathcal {R}_0,\) a \(1.13 \times 10^{-3} \%\) increase in \(\mathcal I_H^*\) and \(1.9 \times 10^{-3} \%\) increase in \(\mathcal I_V^*.\)

These numbers suggest that it might be more favorable for the host population if the loss term due to B-cells switching production of IgM antibodies to IgG antibodies (with a per-capita rate q) is smaller. In other words, because the IgM immune response antibodies are mainly responsible for rapid destruction of virus, for better outcome in disease eradication, rapid destruction of virus within hosts can be more crucial than the life-long immunity (memory) provided by IgG antibodies.

  1. v)

    The IgM immune response decay rate is the next most influential parameter impacting epidemic quantities: a \(1\%\) increase in c causes a negligible \(4.6\times 10^{-4}\%\) increase in \(\mathcal {R}_0\), a \(3.5\times 10^{-4} \%\) increase in \(\mathcal I_H^*\), and a \(5.6 \times 10^{-4} \%\) increase in \(\mathcal I_V^*\).

  2. vi)

    Finally, the killing efficacy parameters \(\theta , \delta \) for both immune responses appear to be least impactful immunological parameters.

It is interesting to note that an increase in the killing efficacy parameters can lead to an increase in epidemic quantities, although just marginally. The biological insight behind this numerical observation is that: because the baseline parameters are in the high disease fatality parameter region, an increase in parameter killing efficacy slows down the viral density growth within-hosts, decreasing disease induced death rates; subsequently, allowing an increase in infectious period, which leads an increase in epidemic quantities (see Fig.7(bottom row subfigures)). In summary, comparing the impact of immune response parameters, the IgG immune activation rate seems to be the least crucial, underscoring the importance of IgM immune parameters in disease eradication over IgG immune response. However this result might be an artifact of how the disease-induced death rate \(\kappa \) is formulated and what the range of the parameter \(\theta \) is. We observed that the value of the IgM killing efficacy parameter largely determines the peak viral load in such a way that, when the IgM killing efficacy is chosen sufficiently large, not only the viral relapse occurs but also there is a larger change in the virus density (see Fig 8). Since the values of \(\theta \) in the chosen range are small, we observe no viral relapse and a smaller impact on the viral dynamics, but a larger reduction in the antibody response dynamics.

3.1 The Impact of the Choice of Linking Parameters on the SA

An important question is, how does the choice of linking functions affect the outcome of the SA performed in Section3? To study this question, we consider different linking functions that have been implemented in the literature (Table 3), and we carry out the corresponding SA using a similar process to that described in Sect. 3. To do that we first re-calibrate the coefficients in the alternative linking functions (\(LP_2\)-columns in Table 3) by solving the following optimization problem:

$$\begin{aligned} \min _{C^*\in (C_L,C_R)}|\mathcal {R}_0(C^*)-\mathcal {R}_0|, \end{aligned}$$

to approximate the basic reproduction number \(\mathcal {R}_0\) at the baseline setting. We have fit the four linking parameters, \(C_0^*,\zeta ^*,\xi ^*, C_\gamma ^*\), one at a time: we vary one linking parameter \(C^*\) over its corresponding range, \((C_L,C_R)\), specified in Table 4, while keeping all the other parameters in the \(LP_1\) forms at their baseline values. This single-variable optimization was repeated for each of the four linking parameters, and it was done using the fminbnd algorithm in MATLAB with TolX = \(10^{-10}\). The resulting estimates are listed in Table 4.

We present the resulting sensitivity indices for these new linking functions in Table 5 (\(LP_2\) columns). Comparing these with the indices for the original linking functions (\(LP_1\) columns), we see the same ordering by magnitude and comparable values for different the paired linking functions. The numerical results suggest that the choice of linking functions, after fitting their linking parameters at the baseline values of the immunological and epidemiological parameters from Tuncer et al. (2016), does not cause important differences in \(\mathcal {R}_0\)-values (Fig.10). We also observe similar patterns in the corresponding curves for final disease abundance \(I^*_H\) and \(I^*_V\) (Appendix A.2).

Furthermore, notice that at the baseline (estimated) parameter values, we have \(\partial \mathcal {R}_0 / \partial r<0\) due to increasing disease induced fatality and decreasing infectious period. Therefore, we classify the diseases according to different ranges of r as: fatal when \(\partial \mathcal {R}_0 / \partial r<0\); mild if \(\partial \mathcal {R}_0 / \partial r>0\). Then, independently of the choice of linking functions, we obtain that when the disease is fatal:

  • The within-host pathogen growth rate r, and pathogen carrying capacity K have the largest impact on the disease reproduction number and final disease abundance size, suggesting that the most efficient control measures are the ones that target reducing within-host viral growth rate, even if virus eradication within hosts cannot succeed.

  • Compared to IgG immune response, the IgM immune response parameters have larger impact on the disease outcomes.

  • Changes in the host immune response parameters have a larger impact on the final infected vector abundance (\(I^*_v\)) than on the final infected host abundance (\(I^*_H\)) because of the indirect transmission route.

Table 3 Different linking parameter formulations
Full size table
Table 4 Re-calibrated parameters for the alternative linking functions (\(LP_2\)) in Table 3
Full size table
Table 5 Local sensitivity indices \(\gamma ^{QOI}_{\varvec{p}}\) with respect to model parameters
Full size table
Table 6 Local epidemic sensitivity indices \(\gamma ^{\mathbf{q}}_{\varvec{p}}\) with respect to within-host model parameters
Full size table

Next we consider the parameter region, where the disease is mild by fixing \(r=1.\) For this case, we observe that SA of QOI to the within-host parameters r, q and \(\delta \) change signs, and the order of the sensitivity indices by magnitude is slightly different, suggesting that:

  • Independently of disease fatality, the per capita parasite growth rate consistently appears to be the most important within-host parameter impacting disease spread. Yet in the mild case (as opposed to fatal case):\(1\%\) increase in r causes \(0.587\%\) increase in \(\mathcal {R}_0\), i.e. increasing the final disease abundance.

  • However \(1\%\) increase in carrying capacity K causes \(0.0137\%\) reduction in \(\mathcal {R}_0.\) This reduction also holds when r is in the “fatal” parameter region. The biological insight behind this finding is that an increase in within-host pathogen carrying capacity causes the host to harbor more pathogens, resulting in increased fatality due to pathogen resource use.

  • An increase in the immune-response activation rates decreases the epidemic quantities. Yet, IgM parameters, compared to IgG, still appear to be more impactful. This suggests that vaccine efficacy can be tested by measuring the increase in IgM and IgG antibody densities comparing data pre- and post-vaccination.

However, in the mild case, other parameters may also see a change in their baseline values corresponding to changing values of r because the parameters in the model may not be independent from each other. For a more rigorous comparison of the mild and the fatal cases, one needs to have another data set representing the mild disease, and fit other parameters as well, which we suggest as valuable future work (Table 6).

3.2 Impact of Viral Replication Rate on Infectability During Distinct Phases of Infection and Epidemic

So far we focused on the impact of immune parameters on the basic reproduction number \(\mathcal R_0\) and the steady state infected host, \(I_H^*\), and infected vector abundance \(I_v^*\). Next, we will show how infectability sensitivity with respect to key immune parameters, \(\gamma ^{\beta (\tau )i_H(\tau ,t)}_\mathbf{p }\) (derived in Appendix A.3), during different outbreak phases and (within-host) infection-times shows heterogeneous impact of virulence-transmission trade-off. Recall that the normalized forward sensitivity index \(\gamma ^{\beta (\tau )i_H(\tau ,t)}_\mathbf{p }\) quantifies how much the infectability of cohort of individuals \(i_H(\tau ,t)\) changes when the parameter \(\mathbf{p} \) is increased by \(1\%.\) In contrast to the reproduction number and steady states, infectability, defined with the epidemic quantity \(\beta (\tau )i_H(\tau ,t)\), is a function of times since outbreak and individual infection start, t and \(\tau \), respectively, allowing us to see within-host effect on the main population level determinant of force of infection which varies with t and \(\tau \). Although sensitivities are calculated for all immune parameters, because the viral replication rate is the fastest evolving parameter determining virulence (Gilchrist and Sasaki 2002; Gulbudak et al. 2017), we highlight here the influence of r on infectability, dependent on both infection-age of cohorts of infected individuals and time since outbreak began on population scale (see Appendix A.3 for other parameter sensitivities and additional figures). In this way we can understand how, in general, a virus might evolve in response to control strategies acting at different epidemic phase and different times during infection course (e.g. isolation of infected cases), or conversely how interventions may be optimally timed given the viral replication rate.

We observe in Fig.5 and 13a that the raw (non-normalized) sensitivity of infectability to viral replication rate, \(\left( \beta \times i_H\right) _{r}(\tau ,t)\), is positive during early infection (small \(\tau \)) and negative after a sufficient amount of time-since infection passes, and the magnitude amplifies at later phases of the epidemic. Indeed, increasing r causes the within-host infection to be more severe and acute, with larger and earlier peaks of both virus and immune response, thereby tending to shift the transmissions caused by an infected individual to occur earlier during infectious period. Certain control strategies, such as contact tracing, self-isolation, case finding and hospitalization, require some time \(\tau \) after infection to act due to delays in symptom onset or tracking, and thus the increased transmissibility during early infection of viral strains with higher r values may make these interventions harder to implement (Browne et al. 2015, 2022, [17]). In this way, a virus may tend to evolve larger replication rates in response to these control strategies, which may increase disease virulence. With this knowledge, public health authorities may prioritize rapid response for active and passive case finding in order to combat evolution of existing virulent strains.

Fig. 5
figure 5

Sensitivity of infectability of host (at distinct epidemic phase at time t) to viral replication rate r (Color figure online)

Full size image

Furthermore, Fig.5 suggests that the change in viral replication rate can impact the infection transmissibility of individuals differently depending on epidemic phase. In particular, in the inset figure in Fig.5b, we observe that when within-host parasite growth rate is increased by \(1\%,\) infection transmissibility increases up to \(4\%\) at epidemic time \(t=40\) (initial epidemic phase) and at \(t=100\) (epidemic die out phase), but only up to \(2.5\%\) during epidemic peak phase (\(t=80\)), with the maximum increase occurring around 1 day post infection and almost 100% decrease in transmissibility during the end of infection because of the more acute infection. This normalized sensitivity and the raw sensitivities show that varying r may dramatically shift infectability to earlier infection ages at the later declining phase of epidemic, making the implications of virulence and control strategies discussed in previous paragraph even more relevant. Note that these results might be due to estimated parameters values for RFVD. More rigorous assessment of immune parameters on the infection transmissibility requires further developed techniques in fitting data from within-host and between-host scales to immuno-epidemiological models.

4 Discussion

There are some good reasons for the use of SA in immuno-epidemiological models. First, infectious diseases operates at both immunological (individual) and epidemiological (population) scales. Therefore, quantifying how the underlying immunological processes affect the dynamics of the epidemic can be naturally approached by SA of multi-scale systems. Secondly, SA of epidemic models is frequently used to assess the impact of control strategies that impact the within-host (or in-vector) virus-and-immune-response dynamics, such as vaccination, drug treatment or bio-control strategies (e.g. Wolbachia-based). A rigorous approach of SA requires the assessment of the sensitivity of epidemic quantities to immunological parameters. Finally, since obtaining population scale data is a challenging task and requires costly surveillance efforts, the parameters estimated from within-host viral and immune-response data that can obtained less costly from laboratory experiments, can be helpful for control strategy development and its assessment. Therefore, to investigate how within-host immune parameters affect disease spread among host population, the SA of the multi-scale models, as the one introduced here, can be informative and give crucial insights on the expected impact of control strategies and their efficacy.

Our results suggest that the within-host per capita parasite growth rate r has the largest impact on disease spread: \(1\%\) increase in r in the parameter region \(r \in [0.1, 2]\) leads to up to \(8\%\) increase in \(\mathcal {R}_0,\) and up to \(1 \%\) increase in \(\mathcal I_H^*,\) and \(1.5 \%\) increase in \(\mathcal I_V^*\), which are significant numbers concerning disease outcomes. Among the rest of the within-host immune parameters none has a significant impact on the epidemic at the population scale. We may still note that, IgM immune response (mainly responsible for the initial rapid destruction of virus within the host) parameters have larger impact compared to those of IgG immune response (responsible for life-long immunity). When the disease is fatal, a \(1\%\) increase in q (\(\hbox {IgM}\rightarrow \) IgG switching rate) leads to an increase in \(\mathcal R_0\) and when it is mild, it decreases the epidemic quantities. This suggests that when the disease is mild, it is more favorable for the host population to have more IgG antibodies that provide life-long immunity. On the other hand, when the disease is fatal, it emphasizes the benefit of a smaller switching rate. Interestingly, increases in the IgM or IgG virus killing efficacy parameters \(\theta ,\) and \( \sigma \) (respectively), both lead to increases in \(\mathcal R_0\), albeit insignificant ones.

In previous studies (Gulbudak et al. 2017; Gilchrist and Sasaki 2002; Tuncer et al. 2016), the choice of linking epidemic parameters as functions of immune variables were different. Therefore to address how sensitivity indices, \(\gamma ^{.}_\mathbf{p },\) change with respect to linking functions chosen, we first non-dimensionalize all linking functions, and then fit the unknown parameters in them using data from prior studies. By doing so, we show that the different linking functions considered do not cause differences in the main outcomes of SA of immunological parameters on epidemic quantities. An interesting numerical result is that even though the unknown parameters in linking functions are fitted at the baseline (estimated) parameters, the consistency on the overlapping parameter values continues over a range of parameter values given as nontrivial intervals around those baseline values.

SA is a common methodological approach to determine the impact of model parameters on the within- and between-host population dynamics quantities. However, these two scales are typically treated separately. Since infectious diseases operate at both scales, in order to understand the crucial mechanisms behind disease dynamics and the impact of targeted control strategies, we need the formulation and analysis of unified models describing disease progression on both scales. One of the biggest challenges for using immuno-epidemiological models is the lack of multi-scale parameter estimations. The sensitivity analysis approach described here is local in nature. Generally, the sensitivity near one choice of parameter values maybe very different from that near a different choice of parameter values, so that a separate sensitivity analysis would need to be performed for every different choice of baseline parameters. However, for the system at hand, we based our parameter estimates on Tuncer et al. (2016), where it was shown that the model is globally structurally identifiable; that is, a unique set of parameters produces the best fit to the data. Therefore, the fitted parameter estimates we use are reliable.

Finally, the model introduced here can be reformulated to incorporate seasonality and temperature-dependent vector growth cycle and vertical transmission in the mosquito population. We plan to extend the present model in the future by incorporating these mechanisms to investigate the impact of seasonality and temperature-dependent vector growth cycle along together with within-host immunological parameters on the disease dynamics at population scale. This might enable us to forecast the impact of implemented control strategies and inform changes to increase their efficacy to help reduce or eradicate the disease. Even though not a vector-host disease, we may also extrapolate implications to the COVID-19 pandemic. In the presence of unprecedented control efforts which included intensive case finding, COVID-19 evolved more virulent variants, such as the Delta variant, which rapidly spread in some countries that were in declining phases of epidemic waves. This evolution of increased viral replication rate contrasts some previous theory that viruses might evolve to be less virulent, and suggest that closer consideration of multi-scale infection-age models may be important to understand virulence-transmission trade-offs and evolution.