# Sensitivity analysis in an Immuno-Epidemiological Vector-Host Model

Hayriye Gulbudak† Department of Mathematics, University of Louisiana at
Lafayette 217 Maxim Doucet Hall, P.O. Box 43568, Lafayette, LA
<EMAIL_ADDRESS>,  Zhuolin Qu  Department of Mathematics,
University of Texas at San Antonio, One UTSA Circle, San Antonio, TX 78249
<EMAIL_ADDRESS>,  Fabio Milner  School of Mathematical and Statistical
Sciences, Arizona State University, 825 Wexler Hall, PO Box 871804, Tempe, AZ
85287<EMAIL_ADDRESS>and  Necibe Tuncer  Department of Mathematical
Sciences, Florida Atlantic University, Science Building, Room 234 777 Glades
Road Boca Raton, FL 33431<EMAIL_ADDRESS>
###### Abstract.

Sensitivity Analysis (SA) is a useful tool to measure the impact of changes in
model parameters on the infection dynamics, particularly to quantify the
expected efficacy of disease control strategies. SA has only been applied to
epidemic models at the population level, ignoring the effect of within-host
virus-with-immune-system interactions on the disease spread. Connecting the
scales from individual to population can help inform drug and vaccine
development. Thus the value of understanding the impact of _immunological_
parameters on epidemiological quantities. Here we consider an age-since-
infection structured vector-host model, in which epidemiological parameters
are formulated as functions of within-host virus and antibody densities,
governed by an ODE system. We then use SA for these immuno-epidemiological
models to investigate the impact of immunological parameters on population-
level disease dynamics such as _basic reproduction number, final size of the
epidemic_ or _the infectiousness at different phases of an outbreak_. As a
case study, we consider Rift Valley Fever Disease (RFVD) utilizing parameter
estimations from prior studies. SA indicates that $1\%$ increase in within-
host pathogen growth rate can lead up to $8\%$ increase in $\mathcal{R}_{0},$
up to $1\%$ increase in steady-state infected host abundance, and up to $4\%$
increase in infectiousness of hosts when the reproduction number
$\mathcal{R}_{0}$ is larger than one. These significant increases in
population-scale disease quantities suggest that control strategies that
reduce the within-host pathogen growth can be important in reducing disease
prevalence.

Keywords: immuno-epidemiological model, sensitivity analysis, Rift Valley
Fever, basic reproduction number, multi-scale model

AMS Subject Classification: 92D30, 92D40

†author for correspondence.

##  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 [20]. 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 [23, 13, 10, 2]. 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 [14]. The parameters for this model are
studied in [23] based on Rift Valley Fever Disease (RVFD) immunological data
and human epidemiological data from the $2006-2007$ Kenya Outbreak [6, 21].
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.
[9, 8] considered an ODE vector-host RVF model to assess the effectiveness of
some control interventions on RVFD. Fischer et al. [5] utilized SA to
investigate the effect of temperate climate on the RVFD dynamics. Mpeshe et
al. [20] 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. [24] 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 Section 2, we present an immuno-
epidemiological model, first introduced in [14, 23], and summarize the
stability and persistence conditions for the disease. In Section 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
subsection 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 [14]. 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:

(2.1) $\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.\ $

The initial values are $P(0)=P_{0}>0,M(0)=M_{0}\geq 0,G(0)=G_{0}\geq 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)_ [14, 23]. 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$ [16]. All parameters and state variables of
this within-host model and their definitions are given in Table 2.1 and Table
2.2. Different dynamics of this immune response model can be found in Figure
2.1.

Figure 2.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 [22].

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
[12]. 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 [12]:

######  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})$.

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 Figure 2.1 and as suggested
by the analytical result above (Theorem 1). In general, within-host dynamics
in arbovirus diseases behave as depicted in Figure 2.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 Figure 2.1(right). In surviving monkeys,
the within-host dynamics imitated the third scenario [19]. 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 Figure 2.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 [22]. Our model can capture these
dynamics as well, as can be seen in Figure 2.1 (right).

At the _population scale_ , we consider a time-since-infection structured
vector-host model, (2.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:

(2.2) $\left\\{\begin{aligned}
\frac{dS_{H}}{dt}&=\Lambda-\beta_{V}S_{H}(t)I_{V}(t)-dS_{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{dR_{H}}{dt}&=\displaystyle\displaystyle\int_{0}^{\infty}{\gamma(\tau)i_{H}(\tau,t)d\tau}-dR_{H}(t),\\\
\frac{dS_{V}}{dt}&=\eta-\left(\int_{0}^{\infty}{\beta_{H}(\tau)i_{H}(\tau,t)d\tau}+\mu\right)S_{V}(t),\\\
\frac{dI_{V}}{dt}&=S_{V}(t)\int_{0}^{\infty}{\beta_{H}(\tau)i_{H}(\tau,t)d\tau}-\mu
I_{V}(t).\end{aligned}\right.$

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 [15, 7], as follows:

(2.3) $\left\\{\begin{aligned}
\beta_{H}(\tau)=\frac{C_{\beta}P(\tau)}{C_{0}+P(\tau)},\\\\[10.00002pt]
\alpha(\tau)=\zeta P(\tau),\\\\[10.00002pt] \kappa(\tau)=\xi
M(\tau),\\\\[10.00002pt]
\gamma(\tau)=C_{\gamma}\frac{G(\tau)}{P(\tau)+\epsilon_{0}}.\end{aligned}\right.$

Fig. A.1 in the Appendix shows how some of the corresponding epidemiological
parameters evolve over the course of host infection. Note that Tuncer et
al.[23] considers mass action term, satisfying 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.2.

Figure 2.2. Schematic illustration of the immuno-epidemiological vector-host model(2.1)-(2.3). Table 2.1. Model parameters: Parameter estimates of the immuno-epidemiological model (2.1)-(2.3) fitted [23]. $\theta$ and $\delta$ are the killing efficiency parameters that re-scale the population immune cells. *$q$ value is chosen differently (fit separately)  |  immunological parameters (within-host model (2.1)) | dimension | baseline | range
---|---|---|---|---
$r$ | Unit pathogen growth rate | $\mbox{day}^{-1}$ | $6.3033$ | $1\sim 10$
$K$ | Carrying capacity of pathogens |  TCID50 | $7.57\times 10^{7}$ | $1.5\times 10^{7}\sim 10^{8}$
$a$ | Unit IgM-activation rate per pathogen | $(\mbox{TCID}_{50}\times\mbox{day})^{-1}$ | $3.93\times 10^{-7}$ | $10^{-8}\sim 5\times 10^{-7}$
$b$ | Unit IgG-activation rate per pathogen | $(\mbox{TCID}_{50}\times\mbox{day})^{-1}$ | $9.58\times 10^{-8}$ | $5\times 10^{-9}\sim 10^{-6}$
$\theta$ | Unit IgM pathogen-destruction rate per-pathogen | $(\mbox{Elisa PP}\times\mbox{day})^{-1}$ | $5\times 10^{-7}$ | $10^{-9}\sim 10^{-4}$
$\delta$ | Unit IgG pathogen-destruction rate per pathogen | $(\mbox{Elisa PP}\times\mbox{day})^{-1}$ | $8.3918\times 10^{-8}$ | $1.7\times 10^{-10}\sim 1.7\times 10^{-5}$
$q$ | Unit IgM-to-IgG production switch rate | $\mbox{day}^{-1}$ | $0.4232^{*}$ | $2\sim 20$
$c$ | Unit IgM decay rate | $\mbox{day}^{-1}$ | $0.1155$ | $0.01\sim 1$
|  epidemic parameters (vector-host model (2.2)) |  |  |
$\Lambda=d$ | Host recruitment rate |  host$\times$day-1 | ${1}/{(365\times 10)}$ |
$\eta=\mu$ | Vector recruitment rate |  vector$\times$day-1 | ${1}/{40}$ |
$\beta_{V}$ | Per-capita-per-infected vector transmission rate |  (vector$\times$ day)-1 | $0.2$ |
$\beta_{H}(\tau)$ | Unit age-of-infection-dependent per host transmission rate |  (host$\times$ day)-1 | - |
$\alpha(\tau)$ | Unit pathogen-induced death rate |  day-1 | - |
$\kappa(\tau)$ | Unit immune-response-induced death rate |  day-1 | - |
$\gamma(\tau)$ |  Unit recovery rate $\tau$ days post-infection |  day-1 | - |
$d$ | Unit natural death rate of hosts |  day-1 | ${1}/{(365\times 10)}$ |
$\mu$ | Unit natural death rate of vectors |  day-1 | ${1}/{40}$ |
|  linking function parameters (2.3) |  |  |
$C_{\beta}$ | transmission efficiency of pathogen infection |  (host$\times$ day)-1 | $0.5365$ |
$C_{0}$ | half-saturation constant in transmission rate | $\mbox{TCID}_{50}$ | $3.03\times 10^{3}$ |
$\zeta$ | unit per host pathogen-lethality rate | $(\mbox{TCID}_{50}\times\mbox{day})^{-1}$ | $6.21\times 10^{-8}$ |
$\xi$ | unit per host immune-response-lethality rate | $(\mbox{Elisa PP}\times\mbox{day})^{-1}$ | $8.51\times 10^{-5}$ |
$C_{\gamma}$ | recovery coefficient | $\frac{\mbox{TCID}_{50}}{\mbox{Elisa PP}\times\mbox{day}}$ | $0.7212$ |
$\epsilon_{0}$ | half-saturation constant in recovery rate |  TCID50 | $7.43\times 10^{-4}$ |
Table 2.2. State Variables  | Meaning
---|---
$P(\tau)$ | Pathogen concentration at infection-age $\tau$
$M(\tau)$ | IgM concentration at infection-age $\tau$
$G(\tau)$ | IgG concentration at infection-age $\tau$
$S_{H}(t)$ | Number of susceptible hosts at time $t$
$i_{H}(\cdot,t)$ | Age-of-infection density of infected hosts at time $t$
$R_{H}(t)$ | Cumulative number of recovered hosts at time $t$
$S_{V}(t)$ | Number of susceptible vectors at time $t$
$I_{V}(t)$ | Number of infected vectors at time $t$

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

(2.4)
$\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},$

where $S^{0}_{H}=\displaystyle\frac{\Lambda}{d}$ and
$S_{V}^{0}=\displaystyle\frac{\eta}{\mu}.$ A detailed analysis of our model
(2.1)-(2.3) including also vector age-of-infection was presented in [12]. 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.2) has a unique endemic
equilibrium (EE) $\mathcal{E}^{*}=(S^{*}_{H},\ i^{*}_{H}(\tau),\ R^{*}_{H},\
S^{*}_{V},\ I^{*}_{V})$ (also presented in [14]), where

(2.5) $\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.$

(2.6) $\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.$

with

(2.7)
$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)=e^{-\displaystyle\int_{0}^{\tau}{(\alpha(s)+\kappa(s)+\gamma(s)+d)ds}},$

and the epidemiological parameters $\beta(\tau),\ \gamma(\tau),\
\alpha(\tau),$ and $\kappa(\tau)$ are given by (2.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 [18]
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 (2.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:

$\displaystyle\gamma^{q}_{p}=\frac{\partial q}{\partial p}\times\frac{p}{q}.$

The baseline parameter values and ranges for RVFD are presented in Table 2.1.
The baseline values were fitted to multi-scale data in [23]: 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, 21]. This
simultaneous fitting of immunological and epidemiological model parameters
induces practical identifiability of model parameters [23].

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. [23], 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 (2.4), we have

(3.1) $\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,$

where $\beta_{H}(\tau,p)$ and $\pi_{H}(\tau,p)$ are defined in (2.3) and
(2.7), respectively, and

$\displaystyle\frac{\partial\beta_{H}(\tau,p)}{\partial p}$
$\displaystyle=\frac{C_{\beta}C_{0}}{(C_{0}+P(\tau))^{2}}\frac{\partial
P(\tau,p)}{\partial p},$ $\displaystyle\frac{\partial\pi_{H}(\tau,p)}{\partial
p}$ $\displaystyle=-\pi_{H}(\tau)\int_{0}^{\tau}\left(\zeta\frac{\partial
P(s,p)}{\partial p}+\xi\frac{\partial M(s,p)}{\partial p}\right.$ (3.2)
$\displaystyle\left.\qquad\qquad+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.$

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

(3.3)
$\frac{\partial}{\partial\tau}\left(\frac{\partial\boldsymbol{x}}{\partial
p}\right)=\frac{\partial}{\partial
p}\left(\frac{\partial\boldsymbol{x}}{\partial\tau}\right),$

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

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

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

$\underbrace{\begin{pmatrix}\frac{\partial}{\partial t}\left(\frac{\partial
P}{\partial a}\right)\\\\[6.99997pt] \frac{\partial}{\partial
t}\left(\frac{\partial M}{\partial a}\right)\\\\[6.99997pt]
\frac{\partial}{\partial t}\left(\frac{\partial G}{\partial a}\right)\\\
\end{pmatrix}}_{\frac{\partial}{\partial
t}\left(\frac{\partial\boldsymbol{x}}{\partial
a}\right)}=\underbrace{\begin{pmatrix}r\left(1-\frac{2P}{K}\right)-\theta
M-\delta G&\ -\theta P&\ -\delta P\\\ aM&\ aP-(q+c)&\ 0\\\ bG&\ q&\
bP\end{pmatrix}}_{\mathrm{\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}}}\underbrace{\begin{pmatrix}\frac{\partial
P}{\partial a}\\\\[6.99997pt] \frac{\partial M}{\partial a}\\\\[6.99997pt]
\frac{\partial G}{\partial a}\\\\[6.99997pt]
\end{pmatrix}}_{\mathrm{\frac{\partial\boldsymbol{x}}{\partial
a}}}+\underbrace{\begin{pmatrix}0\\\ MP\\\ 0\\\
\end{pmatrix}}_{\mathrm{\frac{\partial\boldsymbol{f}}{\partial a}}}.$

By solving the extended system (3.4)-(2.1), we obtain the solutions for
$\boldsymbol{x}$ and the derivatives $\partial\boldsymbol{x}/\partial p$ at
all times $\tau$, which are then utilized to estimate the integrands in (3.1)
and (3.2). We then numerically approximate the integrals by using the
trapezoidal rule. Note that (3.4) 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
(2.5) and (2.6), respectively. The corresponding sensitivity matrices are
derived as follows:

$\displaystyle\frac{\partial
I^{*}_{H}}{\partial\boldsymbol{p}}=\displaystyle\int_{0}^{\infty}\\!\bigg{[}\frac{\partial
i^{*}_{H}(0)}{\partial\boldsymbol{p}}\pi_{H}(\tau)+i^{*}_{H}(0)\frac{\partial\pi_{H}(\tau)}{\partial\boldsymbol{p}}\bigg{]}d\tau,\
\ \ \ \frac{\partial
i^{*}_{H}(0)}{\partial\boldsymbol{p}}=\frac{1}{\mathcal{R}_{0}^{2}}\frac{\partial\mathcal{R}_{0}}{\partial\boldsymbol{p}}\bigg{/}\bigg{(}\frac{1}{d}+\frac{1}{\beta_{V}}\bigg{)},$
$\displaystyle\frac{\partial
I^{*}_{V}}{\partial\boldsymbol{p}}=\frac{1}{(\mathcal{R}_{0}S^{*}_{H})^{2}}\bigg{(}\frac{\partial\mathcal{R}_{0}}{\partial\boldsymbol{p}}S^{*}_{H}-\frac{\mathcal{R}_{0}}{d}\frac{\partial
i^{*}_{H}(0)}{\partial\boldsymbol{p}}\bigg{)}.$

Figure 3.1. _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}}_{\boldsymbol{p}}$ (right y-axis, orange curve) are
plotted against the immune parameter values $p$ over a range of values around
the baseline estimates.

Figure 3.2. _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}}_{\boldsymbol{p}},$ are plotted against the immune
parameter values $p$ over a range of values around the baseline estimates.

In Fig.3.1, the immune-response-dependent epidemiological reproduction number,
$\mathcal{R}_{0},$ along with the response curves for the extended SA,
$\gamma^{\mathcal{R}_{0}}_{\boldsymbol{p}},$ are plotted against the immune
parameter values $p$ over a range of values around the baseline estimates. In
Fig.3.1(a), 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 3.3, where the indices are sorted by magnitude. The
numerical results suggest that:

  * 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.1, the reduction in $\mathcal{R}_{0}$ with increasing within-host
pathogen growth rate $r$ is counter-intuitive, but well known in the
literature [1]. The trade-off between pathogen virulence and infectiousness is
first shown in [11], articulated through a within-host model for directly
transmitted diseases. Later it was extended to vector-borne diseases in [14].
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. A.2). 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),
$1\%$ increase in $\mathcal{I}_{H}^{*}$ (Fig.3.2), and $1.5\%$ increase in
$\mathcal{I}_{V}^{*}$ (Fig.A.4), 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.

  * 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}^{*}$.

  * 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.3.1(c)). 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.

  * 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.

  * 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}^{*}$.

  * 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. A.2 (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 the Fig A.3). 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.

Table 3.1. Different linking parameter formulations Definition | Param. | LP1 | Ref. | LP2 | Ref.
---|---|---|---|---|---
Host transmission rate | $\beta_{H}(\tau)$ | $\displaystyle\frac{C_{\beta}P(\tau)}{C_{0}+P(\tau)}$ | [23] | $\displaystyle\frac{C_{\beta}^{*}P^{2}(\tau)}{(C_{0}^{*})^{2}+P^{2}(\tau)}$ | [14]
Pathogen induced death rate | $\alpha(\tau)$ | $\zeta P(\tau)$ | [23] | $r\zeta^{*}P(\tau)$ | [14, 11]
Antibody induced death rate | $\kappa(\tau)$ | $\xi M(\tau)$ | [23] | $a\xi^{*}M(\tau)P(\tau)$ | [14, 11]
Host recovery rate | $\gamma(\tau)$ | $C_{\gamma}\displaystyle\frac{G(\tau)}{P(\tau)+\epsilon_{0}}$ | [23] | $C_{\gamma}^{*}\displaystyle\frac{G(\tau)}{P(\tau)+\epsilon_{0}^{*}}e^{-P(\tau)/P(0)}$ | [14]
Table 3.2. Re-calibrated parameters for the alternative linking functions ($LP_{2}$) in Table 3.1. The parameters are estimated using nonlinear least squares to match the basic reproduction number $\mathcal{R}_{0}$.  | Linking function parameters | Dimension | Estimated | Range
---|---|---|---|---
$C_{\beta}^{*}$ | transmission efficiency $(=C_{\beta})$ | (host$\times$ day)-1 | $0.5365$ | -
$C_{0}^{*}$ | half saturation in transmission rate | $\mbox{TCID}_{50}$ | $3.030\times 10^{3}$ | $10^{3}\sim 10^{4}$
$\zeta^{*}$ | the pathogen cost coefficient | $(\mbox{TCID}_{50})^{-1}$ | $9.840\times 10^{-9}$ | $10^{-9}\sim 10^{-7}$
$\xi^{*}$ | the immune response cost coefficient | $(\mbox{Elisa PP})^{-1}$ | $2.951\times 10^{-6}$ | $10^{-7}\sim 10^{-5}$
$C_{\gamma}^{*}$ | recovery coefficient | $\frac{\mbox{TCID}_{50}}{\mbox{Elisa PP}\times\mbox{day}}$ | $4.860$ | $0\sim 10$
$\epsilon_{0}^{*}$ | half saturation in recovery rate $(=\epsilon_{0})$ |  TCID50 | $7.43\times 10^{-4}$ | -
Table 3.3. Local sensitivity indices $\gamma^{QOI}_{\boldsymbol{p}}$ with respect to model parameters  LP1 | $\gamma^{\mathcal{R}_{0}}_{\boldsymbol{p}}$ | $\gamma^{\mathcal{I}_{H}^{*}}_{\boldsymbol{p}}$ | $\gamma^{\mathcal{I}_{V}^{*}}_{\boldsymbol{p}}$ | LP2 | $\gamma^{\mathcal{R}_{0}}_{\boldsymbol{p}}$ | $\gamma^{\mathcal{I}_{H}^{*}}_{\boldsymbol{p}}$ | $\gamma^{\mathcal{I}_{V}^{*}}_{\boldsymbol{p}}$
---|---|---|---|---|---|---|---
$r$ | $-0.598$ | $-0.817$ | $-0.720$ | $r$ | $-0.7520$ | $-0.925$ | $-0.906$
$K$ | $-0.0789$ | $-0.0595$ | $-0.0950$ | $K$ | $-0.0801$ | $-0.0607$ | $-0.0964$
$a$ | $-0.0194$ | $-0.0146$ | $-0.0234$ | $a$ | $-0.0205$ | $-0.0155$ | $-0.0247$
$q$ | $1.56\times 10^{-3}$ | $1.13\times 10^{-3}$ | $1.88\times 10^{-3}$ | $q$ | $1.64\times 10^{-3}$ | $1.22\times 10^{-3}$ | $1.97\times 10^{-3}$
$c$ | $4.64\times 10^{-4}$ | $3.50\times 10^{-4}$ | $5.59\times 10^{-4}$ | $c$ | $4.61\times 10^{-4}$ | $3.49\times 10^{-4}$ | $5.54\times 10^{-4}$
$\theta$ | $3.19\times 10^{-6}$ | $2.40\times 10^{-6}$ | $3.84\times 10^{-6}$ | $\theta$ | $5.44\times 10^{-6}$ | $4.12\times 10^{-6}$ | $6.55\times 10^{-6}$
$\delta$ | $4.65\times 10^{-7}$ | $4.39\times 10^{-7}$ | $5.60\times 10^{-7}$ | $\delta$ | $6.12\times 10^{-7}$ | $5.50\times 10^{-7}$ | $7.37\times 10^{-7}$
$b$ | $3.10\times 10^{-7}$ | $2.31\times 10^{-7}$ | $3.74\times 10^{-7}$ | $b$ | $4.21\times 10^{-7}$ | $3.18\times 10^{-7}$ | $5.07\times 10^{-7}$
Table 3.4. Local epidemic sensitivity indices $\gamma^{\mathbf{q}}_{\boldsymbol{p}}$ with respect to within-host model parameters. _Fatal_ | $\gamma^{\mathcal{R}_{0}}_{\boldsymbol{p}}$ | _Mild_ | $\gamma^{\mathcal{R}_{0}}_{\boldsymbol{p}}$
---|---|---|---
$r$ | $-0.598$ | $r$ | $0.587$
$K$ | $-0.0789$ | $K$ | $-0.0137$
$a$ | $-0.0194$ | $a$ | $-0.0029$
$q$ | $1.56\times 10^{-3}$ | $q$ | $-0.00193$
$c$ | $4.64\times 10^{-4}$ | $c$ | $5.54\times 10^{-4}$
$\theta$ | $3.19\times 10^{-6}$ | $b$ | $4.91\times 10^{-7}$
$\delta$ | $4.65\times 10^{-7}$ | $\theta$ | $3.88\times 10^{-7}$
$b$ | $3.10\times 10^{-7}$ | $\delta$ | $-3.65\times 10^{-7}$

###  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 Section 3? To study this question, we consider
different linking functions that have been implemented in the literature
(Table (3.1)), and we carry out the corresponding SA using a similar process
to that described in Section 3. To do that we first re-calibrate the
coefficients in the alternative linking functions ($LP_{2}$-columns in Table
3.1) by solving the following optimization problem:

$\min_{C^{*}\in(C_{L},C_{R})}|\mathcal{R}_{0}(C^{*})-\mathcal{R}_{0}|,$

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 3.2, 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 3.2.

We present the resulting sensitivity indices for these new linking functions
in Table 3.3 ($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.[23], does not cause
important differences in $\mathcal{R}_{0}$-values (Fig.A.5). 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 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.

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.

a) b)

Figure 3.3. _Sensitivity of infectability of host (at distinct epidemic phase
at time $t$) to viral replication rate $r$._

###  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)}_{\textbf{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)}_{\textbf{p}}$ quantifies how much the
infectability of cohort of individuals $i_{H}(\tau,t)$ changes when the
parameter 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 [11, 14], 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.3.3 and Fig.A.8(a) 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 [3, 4, 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 and existent virulent
strains.

Furthermore, Fig.3.3 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.3.3(b), 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$
(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 [14, 11, 23], the choice of linking epidemic parameters as
functions of immune variables were different. Therefore to address how
sensitivity indices, $\gamma^{.}_{\textbf{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 [23], 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.

##  5\. Appendix

Figure A.1. _The evolution of epidemiological parameters over the course of
infection w.r.t. days post infection $\tau.$_ a) Probability of the host being
recovered at infection age $\tau,$ _given that infectious hosts can exit the
infected compartment only through recovery._ b) Transmission rate w.r.t. the
infection age $\tau.$ c) Probability of death occurring due to disease by
infection age $\tau.$ Parameter values are given in Table 2.1 and the
corresponding within-host dynamics is displayed in the right subfigure of
Fig.2.1.

Figure A.2. Top row: $\mathcal{R}_{0}$ vs. pathogen growth rate $r.$ As
pathogen growth rate $r$ increases, the infectious period decreases with
larger peak load and the immune response antibodies become more
prevalent.Middle row: $\mathcal{R}_{0}$ vs. IgM activation rate $a$. The
faster the IgM antibodies activate, the faster they clear the pathogen with
smaller peak load and decreasing infectious period. Bottom row:
$\mathcal{R}_{0}$ vs. killing efficiency of IgM $\theta$. As the immune
response triggered by IgM antibodies become more efficient at killing
pathogens, the pathogen load and subsequent immune response get smaller.
Figure A.3. Impact of killing efficacy $\theta$ on the viral dynamics.

###  A.1. SA curves for final infected vector abundance, $\mathcal{I}_{V}^{*}$

See Fig. A.4

Figure A.4. Impact of immune parameters on the steady-state infected vector
abundance, $\mathcal{I}_{V}^{*}.$ Notice that numerical results for SA for
vector component are similar to SA of $\mathcal{R}_{0},$ and the steady-state
host endemic size $\mathcal{I}^{*}_{H}.$

###  A.2. Impact of using different linking functions on SA:
$\mathcal{I}_{H}^{*}$ and $\mathcal{I}_{V}^{*}$

See Fig. A.6 and Fig. A.7.

Figure A.5. Impact of the choice of linking functions on the SA.

Figure A.6. Impact of the choice of linking functions on the SA.

Figure A.7. Impact of the choice of linking functions on the SA.

(a) (b)  (c)  (d)  (e)  (f)

Figure A.8.  _How infection transmissibility of hosts (with different
infection group $\tau$) changes when quantity of interest (QI) such as an
immune parameter $p$ is increased by $1\%.$_ At initial phase of epidemic
($t=40$ days after epidemic started), the percentage change in infection
transmissibility among infected host groups is the largest when the change is
occurred in parasite growth rate $r,$ comparing to the rest of the immune
parameters.

###  A.3. Derivation of infectability Sensitivity Matrix

By utilizing Clairaut’s theorem, we first obtain the derivative of the system
(2.2) with respect to a model parameter $p$:

(A.1) $\displaystyle\frac{\partial}{\partial t}\left(\frac{\partial
S_{H}}{\partial p}\right)$ $\displaystyle=-\beta_{V}\frac{\partial
S_{H}}{\partial p}I_{V}(t)-\beta_{V}S_{H}(t)\frac{\partial I_{V}}{\partial
p}-d\frac{\partial S_{H}}{\partial p}$ $\displaystyle\frac{\partial}{\partial
t}\left(\frac{\partial i_{H}}{\partial p}\right)$
$\displaystyle+\frac{\partial}{\partial\tau}\left(\frac{\partial
i_{H}}{\partial
p}\right)=-\frac{\partial\left((\alpha(\tau)+\kappa(\tau)+\gamma(\tau)+d)\right)}{\partial
p}i_{H}(\tau,t)-(\alpha(\tau)+\kappa(\tau)+\gamma(\tau)+d)\left(\frac{\partial
i_{H}}{\partial p}\right),$ $\displaystyle\frac{\partial i_{H}(0,t)}{\partial
p}$ $\displaystyle=\beta_{V}\frac{\partial S_{H}}{\partial
p}I_{V}(t)+\beta_{V}S_{H}(t)\frac{\partial I_{V}}{\partial p},$
$\displaystyle\frac{\partial}{\partial t}\left(\frac{\partial R_{H}}{\partial
p}\right)$
$\displaystyle=\int_{0}^{\infty}\left(\frac{\partial\gamma(\tau)}{\partial
p}i_{H}(\tau,t)+\gamma(\tau)\frac{\partial i_{H}(\tau,t)}{\partial
p}\right)d\tau-d\frac{\partial R_{H}}{\partial p},$
$\displaystyle\frac{\partial}{\partial t}\left(\frac{\partial S_{V}}{\partial
p}\right)$ $\displaystyle=-\frac{\partial S_{V}}{\partial
p}\int_{0}^{\infty}{\beta_{H}(\tau)i_{H}(\tau,t)d\tau}-S_{V}(t)\int_{0}^{\infty}\left(\frac{\partial\beta_{H}}{\partial
p}i_{H}(\tau,t)+\beta_{H}(\tau)\frac{\partial i_{H}(\tau,t)}{\partial
p}\right)d\tau-\mu\frac{\partial S_{V}}{\partial p},$
$\displaystyle\frac{\partial}{\partial t}\left(\frac{\partial I_{V}}{\partial
p}\right)$ $\displaystyle=\frac{\partial S_{V}}{\partial
p}\int_{0}^{\infty}{\beta_{H}(\tau)i_{H}(\tau,t)d\tau}+S_{V}(t)\int_{0}^{\infty}\left(\frac{\partial\beta_{H}}{\partial
p}i_{H}(\tau,t)+\beta_{H}(\tau)\frac{\partial i_{H}(\tau,t)}{\partial
p}\right)d\tau-\mu\frac{\partial I_{V}}{\partial p}.$

Following the numerical scheme in [23], we discretize the system using the
Backward Euler difference quotient along the characteristic, and obtain the
following fully-discretized finite difference scheme:

(A.2) $\begin{cases}\displaystyle\frac{S_{H}^{n+1}-S_{H}^{n}}{\Delta
t}&=\Lambda-\beta_{v}S_{H}^{n+1}I_{v}^{n}-dS_{H}^{n+1},\
\displaystyle\frac{i_{H}^{k+1,n+1}-i_{H}^{k,n}}{\Delta
t}=-\left(\alpha_{k+1}+\gamma_{k+1}+d\right)i_{H}^{k+1,n+1},\vspace{1.5mm}\\\
i_{H}^{0,n+1}&=\beta_{v}S_{H}^{n+1}I_{v}^{n+1},\vspace{1.5mm}\\\
\displaystyle\frac{S_{v}^{n+1}-S_{v}^{n}}{\Delta t}&=\eta-
S_{v}^{n+1}\displaystyle\sum_{k=1}^{M}\beta_{H}^{k}i_{H}^{k,n+1}-\mu
S_{v}^{n+1},\ \displaystyle\frac{I_{v}^{n+1}-I_{v}^{n}}{\Delta
t}=S_{v}^{n+1}\displaystyle\sum_{k=1}^{M}\beta_{H}^{k}i_{H}^{k,n+1}-\mu
I_{v}^{n+1}\vspace{1.5mm}\\\
\displaystyle\frac{S_{H_{p}}^{n+1}-S_{H_{p}}^{n}}{\Delta
t}&=-\beta_{v}S_{H_{p}}^{n+1}I_{v}^{n}-\beta_{v}S_{H}^{n+1}I_{v_{p}}^{n}-dS_{H_{p}}^{n+1},\vspace{1.5mm}\\\
\displaystyle\frac{i_{H_{p}}^{k+1,n+1}-i_{H_{p}}^{k,n}}{\Delta
t}&=-\displaystyle\frac{\partial}{\partial
p}\left(\alpha_{k+1}+\gamma_{k+1}+d\right)i_{H}^{k+1,n+1}-\left(\alpha_{k+1}+\gamma_{k+1}+d\right)i_{H_{p}}^{k+1,n+1},\vspace{1.5mm}\\\
i_{H_{p}}^{0,n+1}&=\beta_{v}S_{H_{p}}^{n+1}I_{v}^{n}+\beta_{v}S_{H}^{n+1}I_{v_{p}}^{n},\vspace{1.5mm}\\\
\displaystyle\frac{S_{v_{p}}^{n+1}-S_{v_{p}}^{n}}{\Delta
t}&=-S_{v_{p}}^{n+1}\displaystyle\sum_{k=1}^{M}\beta_{H}^{k}i_{H}^{k,n+1}-S_{v}^{n+1}\displaystyle\sum_{k=1}^{M}\displaystyle\frac{\partial\beta_{H}^{k}}{\partial
p}i_{H}^{k,n+1}-S_{v}^{n+1}\displaystyle\sum_{k=1}^{M}\beta_{H}^{k}i_{H_{p}}^{k,n+1}-\mu
S_{v_{p}}^{n+1},\vspace{1.5mm}\\\
\displaystyle\frac{I_{v_{p}}^{n+1}-I_{v_{p}}^{n}}{\Delta
t}&=S_{v_{p}}^{n+1}\displaystyle\sum_{k=1}^{M}\beta_{H}^{k}i_{H}^{k,n+1}+S_{v}^{n+1}\displaystyle\sum_{k=1}^{M}\displaystyle\frac{\partial\beta_{H}^{k}}{\partial
p}i_{H}^{k,n+1}+S_{v}^{n+1}\displaystyle\sum_{k=1}^{M}\beta_{H}^{k}i_{H_{p}}^{k,n+1}-\mu
I_{v_{p}}^{n+1},\vspace{1.5mm}\\\ \end{cases}$

where the system is discretized on the domain $D=\\{(\tau,t):0\leq\tau\leq
A,\;0\leq t\leq T\\},$ the grid sizes for time ($\Delta t$) and age
($\Delta\tau$) are chosen to be the same $\Delta t=\Delta\tau$, and $T=N\Delta
t$ and $A=M\Delta\tau$. We approximate the value at each grid point
$(\tau_{k},t_{n})$, $\tau_{k}=k\Delta t$, $t_{n}=n\Delta
t\,,k=1\ldots,M,\;n=1,\ldots,N\,$, which are denoted by $S(t_{n})\approx
S^{n},\;i_{H}(\tau_{k},t_{n})\approx
i_{H}^{k,n},\;\beta_{H}(\tau_{k})\approx\beta_{H}^{k}\,$ along with
derivatives $\displaystyle\frac{\partial S_{H}(t_{n})}{\partial p}\approx
S_{H_{p}}^{n},\;\displaystyle\frac{\partial i_{H}(\tau_{k},t_{n})}{\partial
p}\approx i_{H_{p}}^{k,n},\;\displaystyle\frac{\partial S_{v}(t_{n})}{\partial
p}\approx S_{v_{p}}^{n},\;\displaystyle\frac{\partial
I_{v}(\tau_{k},t_{n})}{\partial p}\approx I_{v_{p}}^{k,n}.$ Lastly, the linked
parameters $\beta_{H}(\tau_{k}),\;\alpha(\tau_{k}),\;\gamma(\tau_{k})$ and
their derivatives with respect to the model parameters $p$ are pre-calculated
following the numerical method as described in Section 3.

Furthermore recall that whenever $\mathcal{R}_{0}>1,$ independent from (non-
zero) initial conditions, the solutions converge to the EE $\mathcal{E}^{*}.$
In addition to that, in the previous section, we analytically and numerically
derive the sensitivity index, $\gamma^{i_{H}^{*}}_{\boldsymbol{p}}$ from the
explicit expression of $i_{H}^{*}(\tau),$ given in (2.5).

The subfigures (a)-(d) in Fig.A.8 display the sensitivity,
$\dfrac{d\beta\times i_{H}(\tau)}{d\mathbf{p}}$ of the infection
transmissibility (of infected individuals) with respect to immunological
parameters $\mathbf{p}\in\\{r,a,b,\theta,\delta,K\\}$ (at varying host
infection age $\tau$) $t=40$ days after epidemic starts, which we called
initial phase of epidemic (see Fig.3.3(a)). The inserted subfigures in each
figure display $\%$ of change in host infection transmissibility when QI
increased $1\%$ at infection age $\tau;$ i.e. $\gamma^{\beta\times
i_{H}}_{\mathbf{p}}(\tau).$ In the inserted subfigure of Fig.A.8(a), we
observe that the changes in the immune parameters can modulate the infection
transmissibility of distinct infected host groups significantly. For example,
at the initial phase of the epidemic, in particular $t=40$ days after epidemic
started, $1\%$ increase in-host parasite growth rate $r$ leads up to $4\%$
increase in the infectivity of individuals whom infected approximately at
least $1.6$ and at most $3.6$ days ago. Notice that the infectious period for
hosts is approximately $4$ days. Significant reduction in the host infection
transmissibility when host infected more than $3.6$ days ago (before virus
clearance) is due to increase in disease death rate induced by larger parasite
growth, and instant large immune response in this infectious time interval. On
the other hand, the decrease in the transmission infectivity of individuals
whom got infected up to $1.6$ days ago, might be due to stronger immune
response leading containment of virus in-host, or causing death. At initial
phase of epidemic, the percentage change in infection transmissibility among
infected host groups is the largest when the change is occurred in parasite
growth rate $r,$ comparing to the rest of the immune parameters $a$ in
Fig.A.8(b), $b$ in Fig.A.8(c), $\theta$ Fig.A.8(d), $\delta$ in Fig.A.8(e),
$K$ Fig.A.8(f).

c) d)

Figure A.9. _Sensitivity of infection transmissibility at distinct epidemic
time $t$ to immune parameters $a$, and $\theta.$_

## Acknowledgement

The authors thank two anonymous reviewers for their helpful comments and
feedback on the manuscript, and Mac Hyman of Tulane University for his helpful
discussions. Hayriye Gulbudak was supported by NSF grant (DMS-1951759) and
Simons Foundation/SFARI(638193). The research of Necibe Tuncer was partially
supported by NSF grant (DMS-1951626).

## References

  * [1] R. M. Anderson and R. May.  Coevolution of hosts and parasites.  Parasitology, 85(2):411–426, 1982.
  * [2] O. Angulo, F. Milner, and L. Sega.  A sir epidemic model structured by immunological variables.  Journal of Biological Systems, 21(04):1340013, 2013.
  * [3] C. Browne, H. Gulbudak, and G. Webb.  Modeling contact tracing in outbreaks with application to ebola.  Journal of theoretical biology, 384:33–49, 2015.
  * [4] C. J. Browne, H. Gulbudak, and J. C. Macdonald.  Differential impacts of contact tracing and lockdowns on outbreak size in covid-19 model applied to china.  Journal of Theoretical Biology, 532:110919, 2022.
  * [5] E. A. Fischer, G.-J. Boender, G. Nodelijk, A. A. De Koeijer, and H. J. Van Roermund.  The transmission potential of rift valley fever virus among livestock in the netherlands: a modelling study.  Veterinary research, 44(1):58, 2013.
  * [6] C. for Disease Control, P. (CDC, et al.  Rift valley fever outbreak–kenya, november 2006-january 2007.  MMWR. Morbidity and mortality weekly report, pages 73–76, 2007\.
  * [7] C. Fraser, K. Lythgoe, G. E. Leventhal, G. Shirreff, T. D. Hollingsworth, S. Alizon, and S. Bonhoeffer.  Virulence and pathogenesis of hiv-1 infection: an evolutionary perspective.  Science, 343(6177):1243727, 2014.
  * [8] H. Gaff, C. Burgess, J. Jackson, T. Niu, Y. Papelis, and D. Hartley.  Mathematical model to assess the relative effectiveness of rift valley fever countermeasures.  International Journal of Artificial Life Research (IJALR), 2(2):1–18, 2011.
  * [9] H. D. Gaff, D. M. Hartley, and N. P. Leahy.  An epidemiological model of rift valley fever.  Electronic Journal of Differential Equations (EJDE)[electronic only], 2007:Paper–No, 2007.
  * [10] A. Gandolfi, A. Pugliese, and C. Sinisgalli.  Epidemic dynamics and host immune response: a nested approach.  Journal of mathematical biology, 70(3):399–435, 2015.
  * [11] M. A. Gilchrist and A. Sasaki.  Modeling host–parasite coevolution: a nested approach based on mechanistic models.  Journal of Theoretical Biology, 218(3):289–308, 2002.
  * [12] H. Gulbudak.  An immuno-epidemiological vector-host model with within-vector viral kinetics.  Journal of Biological Systems, 28(02):233–275, 2020.
  * [13] H. Gulbudak and C. J. Browne.  Infection severity across scales in multi-strain immuno-epidemiological dengue model structured by host antibody level.  Journal of mathematical biology, 80(6):1803–1843, 2020.
  * [14] H. Gulbudak, V. L. Cannataro, N. Tuncer, and M. Martcheva.  Vector-borne pathogen and host evolution in a structured immuno-epidemiological system.  Bulletin of mathematical biology, 79(2):325–355, 2017.
  * [15] A. Handel and P. Rohani.  Crossing the scale from within-host infection dynamics to between-host transmission fitness: a discussion of current assumptions and knowledge.  Philosophical Transactions of the Royal Society of London B, 370(1675):20140302, 2015.
  * [16] T. Honjo, K. Kinoshita, and M. Muramatsu.  Molecular mechanism of class switch recombination: linkage with somatic hypermutation.  Annual review of immunology, 20(1):165–196, 2002.
  * [17] J. Macdonald, C. Browne, and H. Gulbudak.  Modelling covid-19 outbreaks in usa with distinct testing, lockdown speed and fatigue rates.  Royal Society Open Science, 8(8):210227.
  * [18] P. Magal and C. McCluskey.  Two-group infection age model including an application to nosocomial infection.  SIAM Journal on Applied Mathematics, 73(2):1058–1095, 2013.
  * [19] J. Morrill, F. Knauert, T. Ksiazek, J. Meegan, and C. Peters.  Rift valley fever infection of rhesus monkeys: implications for rapid diagnosis of human disease.  Research in virology, 140:139–146, 1989.
  * [20] S. C. Mpeshe, H. Haario, and J. M. Tchuenche.  A mathematical model of rift valley fever with human host.  Acta biotheoretica, 59(3-4):231–250, 2011.
  * [21] P. Munyua, R. M. Murithi, S. Wainwright, J. Githinji, A. Hightower, D. Mutonga, J. Macharia, P. M. Ithondeka, J. Musaa, R. F. Breiman, et al.  Rift valley fever outbreak in livestock in kenya, 2006–2007.  The American journal of tropical medicine and hygiene, 83(2_Suppl):58–64, 2010.
  * [22] M. Pepin, M. Bouloy, B. H. Bird, A. Kemp, and J. Paweska.  Rift valley fever virus (bunyaviridae: Phlebovirus): an update on pathogenesis, molecular epidemiology, vectors, diagnostics and prevention.  Veterinary research, 41(6):61, 2010.
  * [23] N. Tuncer, H. Gulbudak, V. L. Cannataro, and M. Martcheva.  Structural and practical identifiability issues of immuno-epidemiological vector–host models with application to rift valley fever.  Bulletin of mathematical biology, 78(9):1796–1827, 2016.
  * [24] Y. Xiao, J. C. Beier, R. S. Cantrell, C. Cosner, D. L. DeAngelis, and S. Ruan.  Modelling the effects of seasonality and socioeconomic impact on the transmission of rift valley fever virus.  PLoS neglected tropical diseases, 9(1), 2015.