# Vector-borne diseases with non-stationary vector populations: the case of
growing and decaying populations

Àlex Giménez-Romero  Rosa Flaquer-Galmés  Manuel A. Matias<EMAIL_ADDRESS>csic.es IFISC (CSIC-UIB), Instituto de Física Interdisciplinar y Sistemas
Complejos, Campus UIB
E-07122 Palma de Mallorca, Spain  Grup de Física Estadística, Departament de
Física. Facultat de Ciències, Universitat Autònoma de Barcelona, 08193
Bellaterra (Barcelona), Spain

###### Abstract

Since the last century, deterministic compartmental models have emerged as
powerful tools to predict and control epidemic outbreaks, in many cases
helping to mitigate their impacts. A key quantity for these models is the so-
called Basic Reproduction Number, $R_{0}$, that measures the number of
secondary infections produced by an initial infected individual in a fully
susceptible population. Some methods have been developed to allow the direct
computation of this quantity provided that some conditions are fulfilled, such
that the model has a pre-pandemic disease-free equilibrium state. This
condition is only fulfilled when the populations are stationary. In the case
of vector-borne diseases, this implies that the vector birth and death rates
need to be balanced, what is not fulfilled in many realistic cases in which
the vector population grow or decrease. Here we develop a vector-borne
epidemic model with growing and decaying vector populations and study the
conditions under which the standard methods to compute $R_{0}$ work and
discuss an alternative when they fail. We also show that growing vector
populations produce a delay in the epidemic dynamics when compared to the case
of the stationary vector population. Finally, we discuss the conditions under
which the model can be reduced to the SIR model with fewer compartments and
parameters, which helps in solving the problem of parameter unidentifiability
of many vector-borne epidemic models.

###### keywords:

Epidemics , Vector-borne diseases , Basic Reproduction Number , Compartmental
model , Mathematical model

††journal:

##  1 Introduction

Vector-borne diseases are caused by infectious agents transmitted by living
organisms, called vectors, frequently arthropods. These diseases represent a
significant threat to global human health [Athni et al., 2020], causing
diseases such as malaria, dengue, yellow fever, Zika, trypanosomiasis and
leishmaniasis [Schumacher & Campbell, 2018]. Vector-borne human diseases are
responsible of more than 17% of all human infectious diseases, causing
millions of cases and more than $700\,000$ deaths annually [WHO, ]. Moreover,
crop production and farm profitability are also affected by bacterial [Huang
et al., 2020] and virus [Bragard et al., 2013] vector-borne diseases. Some
examples are the Pierce’s Disease of grapevines, that has resulted in an
annual cost of approximately $100 million in California alone [Tumber et al.,
2014], the olive quick decline syndrome, which could cause about $5 and 17
billion of loss in Italy and Spain over the next 50 years in the absence of
disease control measures [Schneider et al., 2020] and the multiple diseases
caused by viruses [Rybicki, 2015], with diseases like the tobacco mosaic,
tomato spotted wilt, etc. and transmitted by aphids and other vectors.

Compartmental deterministic models, e.g. the well known SIR model [Kermack &
McKendrick, 1927], have been widely used in the modelling of vector-borne
diseases after the work of Ross and Macdonald [Macdonald, 1957], that, opened
the way to controlling malaria outbreaks by acting on the vectors of the
disease (the Anopheles mosquito). These models consider that both host and
vector populations can be divided into different compartments describing
different states of the individuals, such as susceptible, infected or dead
[Brauer, 2008], and the time-evolution of these compartments is expressed as a
system of ordinary differential equations, defining a dynamical system.
Compartmental models provide a mean-field description, that imply well-mixed
(in practice spatially homogeneous) populations. The well mixed approximation
will be valid whenever the mean distance among hosts is smaller than the
mixing length of vectors before they die. In the case of vector-borne diseases
it is also equivalent to every vector effectively interacting with all the
hosts and every host with all the vectors. A mean-field description is not
always valid in spatially extended systems, but still it is normally the first
step before writing a spatially explicit description.

The most relevant piece of information about a disease is whether an epidemic
outbreak will take place. The basic reproduction number, $R_{0}$, measures the
number of secondary infections caused by an initial infected individual in a
fully susceptible population, defining the epidemic threshold [Anderson, 1991,
van den Driessche, 2017], that determines the emergence (or not) of an
outbreak. If $R_{0}>1$ an epidemic outbreak will occur, while there will be no
outbreak otherwise. The standard way of determining $R_{0}$ in deterministic
compartmental models assumes the existence of an initial disease-free (pre-
pandemic) equilibrium, represented by the absence of infected hosts and
vectors [Laukó, 2006, Kamgang & Sallet, 2008]. Some standard methods based on
the linear stability condition of this equilibrium have been developed to
allow the direct computation of $R_{0}$, such as the Next-Generation Matrix
(NGM) method [Diekmann et al., 2010].

In the case of vector-borne diseases models typically assume that populations
(both hosts and vectors) do not change with time (see e.g. [Macdonald, 1957,
Brauer et al., 2016]), assuming equal birth and death rates. This guarantees
the existence of a disease-free equilibrium and the proper use of standard
methods to determine $R_{0}$. However, this assumption could be far from
reality in several pathosystems. In particular, the interaction between
temperature, precipitation variations and other factors may lead to strong
variations in the vector population [Garms et al., 1979, Rocklöv & Dubrow,
2020], implying that the pre-pandemic state may not be an equilibrium state
and that standard methods cannot be applied.

Some authors have explicitly considered more general cases in which the
demographic rates are not identical or even time-dependent with a given
periodicity. In the case of unbalanced birth and death rates one obtains an
asymptotic stationary vector population, but not an initial disease-free
equilibrium, except in the case that the initial vector population coincides
with the asymptotic value. Nevertheless, the basic reproduction number of
these models is often computed by means of the standard methods from the
asymptotic state of the population [Wei et al., 2008, Lashari & Zaman, 2011,
Shah & Jyoti, 2013, Zhao et al., 2020, Esteva & Vargas, 1998]. The use of
these methods is supported by the fact that the asymptotic dynamics of the
system under study can be described by the subsystem in which the vector
population is in its stationary value [Thieme, 1992, Castillo-Chavez & Thieme,
1995]. Indeed in [Esteva & Vargas, 1998] it is explicitly pointed out that “it
is enough to study the asymptotic behaviour of the model”. In the case of
periodic demographic rates it has been shown that the time-averaged basic
reproduction number defines the epidemic threshold under some circumstances
[Wesley & Allen, 2009] and even a generalisation of the NGM method [Diekmann &
Heesterbeek, 2000] has been developed for these cases [Bacaër & Guernaoui,
2006, Bacaër, 2007]

Compartmental models of vector-borne diseases have another feature that may
hinder their practical applicability. Namely the fact that these models have
many compartments, having to describe both hosts and vectors, and as a
consequence a relatively large number of parameters. This may lead to an issue
known as parameter identifiability and uncertainty [Chowell, 2017], depending
on the available data, that is more likely to be found in models with many
compartments and parameters [Roosa & Chowell, 2019]. Usually, parameter
estimation procedures are needed to connect the models with disease data,
mainly using incidence or prevalence over time in the host population.
Unfortunately, under many circumstances the underlying model parameters are
unidentifiable, so that many different sets of parameter values produce the
same model fit [Kao & Eisenberg, 2018]. Moreover, these parameters can be
really difficult to determine from the available experimental data.
Nevertheless, in some cases, mathematical manipulations can be performed to
reduce the model complexity using exact or approximate relations [Giménez-
Romero et al., 2021]. In such cases, the number of parameters of the models
can be usually reduced in terms of new parameters defined as combinations of
some of the original parameters.

The plan of this paper is as follows. In Section 2 we develop a compartmental
model of vector-borne transmitted diseases with constant, but different, birth
and death rates for the vectors, which allow us to describe growing and
decaying vector populations. For simplicity, the model assumes that there is
no horizontal (host to host) direct transmission, which is also a realistic
assumption in some cases, like plant diseases. Section 3 contains the main
results of the study. In particular, we show how and when traditional methods
fail to estimate the $R_{0}$ of the model and provide an alternative way to
compute it. It turns out that the validity of the standard methods depends,
among other things, on some time-scales of the model. Furthermore, we discuss
and apply some approximations that allow to reduce the model in favour of
simpler ones, with both fewer compartments and fewer parameters. In
particular, we show that if some of the parameters fulfil certain conditions,
it is possible to reduce the original model with five compartments and four
parameters to a SIR model, with three compartments and two parameters. It is
expected that model reductions like this one significantly help in solving
possible problems of parameter unidentifiability that plague these models. It
is interesting to note that a model in which hosts do not interact directly,
but only through vectors, in a certain limit becomes described as if hosts
would infect one to each other, what is assumed in some studies without
suitable confirmation. Finally, the main concluding remarks of the study are
presented in Section 4.

##  2 The model

The compartment model for vector-borne diseases that we will use to illustrate
the points to be discussed in this study consists of $5$ compartments, $3$ of
which describe the host population (susceptible, $S_{H}$, infected, $I_{H}$,
and removed, $R_{H}$), while the other $2$ describe the vector population
(susceptible $S_{V}$ and infected vectors, $I_{V}$). Thus, we consider that
the pathogen affects only the hosts and do not consider exposed compartments.
In addition, no horizontal (direct host to host) or vertical (mother to
offspring for vectors) transmission is assumed. The model could be also
generalised to include an exposed host compartment and the above mentioned
transmission modes, which would hinder the theoretical analysis without
altering the qualitative conclusions of the study. Anyhow, the absence of
horizontal transmission would make the model adequate to study vector-borne
diseases in which host-to-host infection is rare, like malaria, and most
phytopathologies, although not for, say, zika. Finally, we do not consider
host recruitment neither infected hosts becoming susceptible again.

The model is defined according to the following processes,

$S_{H}+I_{V}\stackrel{{\scriptstyle\beta}}{{\rightarrow}}I_{H}+I_{V}\quad
I_{H}\stackrel{{\scriptstyle\gamma}}{{\rightarrow}}R_{H}\quad
S_{V}+I_{H}\stackrel{{\scriptstyle\alpha}}{{\rightarrow}}I_{V}+I_{H}\quad
S_{V}\stackrel{{\scriptstyle\mu}}{{\rightarrow}}\varnothing\quad
I_{V}\stackrel{{\scriptstyle\mu}}{{\rightarrow}}\varnothing\ ,$ (1)

which are graphically described in Fig. 1, being the birth of new susceptible
vectors described as a source term. Thus, the host-vector compartmental model
is written as,

$\displaystyle\dot{S}_{H}$ $\displaystyle=-\beta S_{H}I_{v}/N_{H}$ (2)
$\displaystyle\dot{I}_{H}$ $\displaystyle=\beta S_{H}I_{v}/N_{H}-\gamma I_{H}$
$\displaystyle\dot{R}_{H}$ $\displaystyle=\gamma I_{H}$
$\displaystyle\dot{S}_{v}$ $\displaystyle=\delta C-\alpha S_{v}I_{H}/N_{H}-\mu
S_{v}$ $\displaystyle\dot{I}_{v}$ $\displaystyle=\alpha S_{v}I_{H}/N_{H}-\mu
I_{v}\ ,$

where a standard incidence [Martcheva, 2015] has been considered.

The model describes infection of susceptible hosts, $S_{H}$, at a rate $\beta$
through their interaction with infected vectors, $I_{v}$, while susceptible
vectors, $S_{v}$, are infected at a rate $\alpha$ through their interaction
with infected hosts $I_{H}$. Infected hosts exit the infected compartment at
rate $\gamma$, while infected vectors stay infected the rest of their
lifetime, as we consider that the pathogen does not affect them, as it is
customary. Vectors die naturally (or disappear from the population by some
mechanism) at rate $\mu$ and are born (appear) at a constant rate $\delta$
being susceptible. The constant term $C$ sets the scale of the stationary
value of the vector population. Fig. 1 shows an schematic representation of
the model and we refer to [Brauer et al., 2016] for a similar model of vector-
borne diseases. The mentioned model includes exposed compartments and
horizontal transmission, but assumes that the birth and death rate of vectors
are identical, and thus, the population does not change with time and stays as
fixed by the initial condition.

Figure 1: Schematic representation of the model n Eq. 2. Boxes are the
compartments in which the population is divided, solid arrows represent
changes in state (so transitions between compartments), and dashed arrows
depict the crossed interaction between hosts and vectors.

###  2.1 Preliminary analysis of the model

From Eq. 2 it is straightforward to verify that the population of hosts
remains constant over time, $N_{H}=S_{H}+I_{H}+R_{H}$, while the vector
population fulfils,

$\dot{N}_{v}=\dot{S}_{v}+\dot{I}_{v}=-\mu\left(S_{v}+I_{v}\right)+\delta
C=-\mu N_{v}+\delta C\ ,$ (3)

which can be solved to yield,

$N_{v}(t)=\frac{\delta}{\mu}C+\left(N_{v}(0)-\frac{\delta}{\mu}C\right)e^{-\mu
t}\ .$ (4)

From Eq. 4 the stationary value for the vector population, $N_{v}^{*}$, can be
computed as,

$N_{v}^{*}=\lim_{t\to\infty}N_{v}(t)=\frac{\delta}{\mu}C\ .$ (5)

Thus, if the initial population of vectors is below (above) the stationary
value, the vector population will grow (decrease) until it reaches the
stationary value. On the other hand, if $N_{v}(0)=N_{v}^{*}=\delta C/\mu$ the
initial population of vectors is already at the stationary state. The initial
condition for the vector population can be written in terms of its stationary
value Eq. 5, $N_{v}(0)=fN_{v}^{*}$, where both $f<1$ and $f>1$ are possible,
so that one gets,

$N_{v}(t)=N_{v}^{*}\left[1+\left(f-1\right)e^{-\mu t})\right]\ .$ (6)

We note that vector-borne disease models that assume constant vector
populations (e.g.[Brauer et al., 2016]) can be recovered by setting
$\delta=\mu$ and $C=N_{v}(0)$, so that any initial condition for the vector
population is stationary, i.e. $\dot{N}_{v}=0$ in Eq. 3 and
$N_{v}(t)=N_{v}(0)$.

##  3 Results

###  3.1 The effect of non-stationary vector populations into the epidemic
threshold and disease dynamics

Let us start with the case in which the birth $\delta$ and death $\mu$ vector
rates are identical and the initial condition of the vector population is
already at its stationary value, $N_{v}(0)=N_{v}^{*}$, case in which the total
vector population remains unchanged, as discussed in Section 2.1. In such a
case, the initial disease free state of the model, given by
$I_{H}(0)=I_{v}(0)=0$, is a fixed point (equilibrium state) of the dynamical
system Eq. 2 independently of the other initial conditions for the host and
vector populations. This allows the definition of the basic reproduction
number, $R_{0}$, using standard methods such as linear stability analysis or
the Next-Generation Matrix (NGM) method [Diekmann et al., 2010] (see A).

The situation is different when $\delta\neq\mu$, as then the total vector
population will vary with time if the initial condition, $N_{v}(0)$, is not
identical to the asymptotic value at large times, $N_{v}^{*}$, and thus an
initial disease-free state is not an equilibrium (fixed point) of the model.
However, in the literature it is customary to apply the standard techniques,
i.e. NGM, to compute $R_{0}$ using the vector population in the asymptotic
state, that is the post-pandemic disease-free equilibrium [Wei et al., 2008,
Lashari & Zaman, 2011, Shah & Jyoti, 2013, Zhao et al., 2020, Esteva & Vargas,
1998]. The use of these methods is supported by the fact that the asymptotic
dynamics of the model converges to the dynamics of the subsystem where the
vector population is stationary [Thieme, 1992, Castillo-Chavez & Thieme,
1995]. In both cases the basic reproduction number is given by,

$R_{0}=\frac{\beta\alpha}{\mu\gamma}\frac{S_{H}(0)}{N_{H}}N_{v}^{*}\ .$ (7)

As usual, $R_{0}$ accounts for the number of secondary infections produced by
an infected individual in one generation and controls the threshold behaviour
of the model: for $R_{0}<1$ the epidemic dies out and for $R_{0}>1$ an
outbreak occurs. By one generation we refer to the typical time in which new
infections can be produced, being the generation time in our model,

$t_{g}=1/\gamma+1/\mu\ .$ (8)

Figure 2: Numerical verification of the predictive power of the basic
reproduction number relation Eq. 7, by plotting the final size of the
epidemic, $R(\infty)/N_{H}$ as function of $R_{0}$. In panels (a),(b) the
initial vector population is in the stationary value, in panels (c),(d) is
below, $N_{V}^{*}/10$, and in panels (e),(f) above, $10N_{V}^{*}$. Panels
(a),(c),(e) show realisations for different $\gamma$ values with a fixed
$\mu=1$ baseline value. Panels (b),(d),(f) show realisations for different
$\mu$ values with a fixed $\gamma=1$ baseline value.

Here we show that this result can only be used in the more general case
provided that some conditions are fulfilled, but not always. Let us first
illustrate the different regimes in which Eq. 7 is predictive or not about the
onset of the epidemic. In Fig. 2 the final size of the epidemic,
$R_{\infty}/N_{H}$, is plotted as a function of $R_{0}$. Fig. 2(a)-(d) show
that Eq. 7 does indeed regulate the onset of an epidemic when the initial
vector population is in its stationary value or below it. This result is
general and does not depend on the time-scales of the system, $1/\gamma$ and
$1/\mu$, and so all curves in these panels behave similarly. In contrast, Fig.
2(e)-(f) shows how the $R_{0}$ calculated according to Eq. 7 is not predictive
of the onset of epidemic outbreak, this happening when the initial vector
population is larger than the stationary value. Thus, for $R_{0}<1$ (computed
using Eq. 7) severe outbreaks appear, yielding mortalities even above 80% of
the total population. However, one can observe that as $\mu$ is increased, or
$\gamma$ decreased, the predictive power of Eq. 7 is progressively recovered.

The observed behaviour and the regimes of validity of Eq. 7 can be explained
readily. In essence, only if the vector population reaches its stationary
value before infected hosts have produced new infections the onset of an
epidemic can be characterised by Eq. 7, but not otherwise. Given these
different behaviours, let us discuss separately the cases $f>1$ and $f<1$,
with $N_{v}(0)=fN_{v}^{*}$, namely when the initial vector population is above
and below its stationary value, this is, decaying and growing vector
populations towards the asymptotic value.

Let us start with the case of decaying vector populations, $f>1$. From Eq. 6,
the time to approach the stationary value, $t^{*}$, fulfils

$\left(1+\epsilon\right)N_{v}^{*}=N_{v}^{*}\left[1+(f+1)e^{-\mu t^{*}}\right]\
,$ (9)

where $\epsilon\to 0$ is a small parameter controlling the amount by which the
vector population differs from its asymptotic value at time $t^{*}$. Thus, the
time to approach the stationary value, with precision $\epsilon$, is given by

$t^{*}=-\frac{1}{\mu}\ln(\frac{\epsilon}{f-1})=\frac{1}{\mu}\absolutevalue{\ln{\frac{\epsilon}{f-1}}}\
,$ (10)

where the last equality assumes that the small parameter $\epsilon$ satisfies
$\epsilon<(f-1)$, $f>1$.

If the vector population reaches its stationary value before infected hosts
have had time to generate new infections then $R_{0}$ as determined from Eq. 7
is a good prediction of the onset for an epidemic, what is equivalent to the
condition that $t^{*}$ is much smaller than the hosts infectious period,
$t^{*}\ll 1/\gamma$,

$\frac{1}{\gamma}\gg\frac{1}{\mu}\absolutevalue{\ln{\frac{\epsilon}{f-1}}}\quad\textrm{or}\quad\frac{\mu}{\gamma}\gg\absolutevalue{\ln{\frac{\epsilon}{f-1}}}\
.$ (11)

Otherwise, Eq. 7 will not be predictive of the epidemic onset, and as shown in
Fig. 2(e-f) one may have outbreaks with a substantial final size with
$R_{0}<1$.

In the case of growing vector populations, $f<1$, if $R_{0}<1$ an outbreak
cannot occur at all, because $R_{0}$ is calculated with the asymptotic
population, $N_{v}^{*}$, that is larger that the vector population at any
finite time, $N_{v}(t)<N_{v}^{*}\ \forall t$, and so the threshold condition
is never attained. In the $R_{0}>1$ case the behaviour will be richer, and it
will depend on the initial condition, $N_{v}(0)$. One can define an
instantaneous basic reproductive number,

$R_{0}^{(i)}(t)=\frac{\beta\alpha}{\mu\gamma}\frac{S_{H}(0)}{N_{H}}N_{v}(t)=R_{0}\frac{N_{v}(t)}{N_{v}^{*}}\
,$ (12)

using $N_{v}(t)$ instead of $N_{v}^{*}$, with $R_{0}^{(i)}(t)<R_{0}\ \forall
t$ because the vector population grows. In particular, if $R_{0}^{(i)}(0)>1$
there will be an outbreak occurring for short times, and the population of
infected hosts will start growing. If instead, $R_{0}^{i}(0)<1$, and as
$R_{0}>1$ with $R_{0}$ being calculated with the asymptotic state, there must
be an intermediate time, say $t_{D}$, for which $R_{0}^{(i)}(t_{D})=1$. Thus,
from $t>t_{D}$ an outbreak will occur, not initially but after a finite time,
that induces a delay in the outbreak, and the infected host population will
start growing.

The difference between the original and the delayed dynamics stems from the
waiting time to reach $R_{0}^{(i)}=1$, $t_{D}$, plus the non-linear effect
associated to a new initial condition for the epidemic outbreak at $t_{D}$.
Thus, in the case that $R_{0}>1$ and $R_{0}^{(i)}(0)<1$, from Eq. 12 and Eq. 6
we can analytically approximate the delay as the time needed to reach
$R_{0}^{(i)}(t_{D})=1$,

${1+(f-1)e^{-\mu t_{D}}}=\frac{1}{R_{0}}\ ,$ (13)

which yields the relation,

$t_{D}=-\frac{1}{\mu}\ln\left[\frac{1-R_{0}}{\left(f-1\right)R_{0}}\right]\ ,$
(14)

where the argument of the logarithm is always positive because $R_{0}>1$ and
$f<1$. Eq. 14 is only valid if $f<1/R_{0}$, for $R_{0}^{(i)}(0)=fR_{0}<1$, as
if otherwise $R_{0}^{(i)}>1$ the outbreak would already occur initially.

Figure 3: Numerical study of the delay induced by growing vector populations.
(a) Comparison of hosts dynamics for a stationary vector population ($f=1$)
and a growing vector population ($f=10^{-5}$). (b) Time delay as function of
$f$ for different values of the basic reproduction number $R_{0}$. (c) Time
delay as function of the vector natural death rate. (d) Time delay as function
of the basic reproduction number, $R_{0}$, with $f=10^{-5}$.

From Eq. 14 one can see that when the initial vector population is far enough
from its stationary value, $f\rightarrow 0$, the delay saturates to a constant
value, instead of increasing. This is,

$\lim_{f\to
0}t_{D}=-\frac{1}{\mu}\ln(\frac{R_{0}-1}{R_{0}})=\frac{1}{\mu}\ln(\frac{R_{0}}{R_{0}-1})\
.$ (15)

In addition, for increasing values of the basic reproduction number, $R_{0}$,
the delay tends to vanish, and from Eq. 15. This is,

$\lim_{R_{0}\to\infty}t_{D}=\frac{1}{\mu}\ln(1)=0\ ,$ (16)

where the limit $f\rightarrow 0$ is taken simultaneously to guarantee that
$R_{0}^{(i)}(0)=fR_{0}<1$. On the other hand the delay, $t_{D}$, scales with
the vectors lifetime,

$t_{D}\sim\frac{1}{\mu}=\tau_{v}\ .$ (17)

Fig. 3(a) shows an example of the time delay caused in the hosts dynamics when
the vector population grows from an initial condition far from the stationary
value. In Fig. 3(b) we can qualitatively observe that all the predicted
properties of the delay are fulfilled, namely, the time delay saturates for
low $f$ values and decreases with increasing $R_{0}$. Although the analytical
expression (black dashed line) is clearly not exact due to nonlinear effects,
Eq. 14 captures the basic trends of the time delay, $t_{D}$. This is clear
from Fig. 3(c), that shows that the delay scales with $1/\mu$ and in Fig. 3(d)
that shows that the delay tends to $0$ in the limit $R_{0}\rightarrow\infty$,
in agreement with the prediction of Eq. 16.

###  3.2 The basic reproduction number for non-stationary vector populations

As shown in the previous section, traditional methods to compute the basic
reproduction number fail in the case of epidemic models with decaying vector
populations, $f>1$, unless the time scale of vector population fulfils several
conditions, as illustrated in Section 3.1. Here we introduce an effective,
average definition of $R_{0}$, useful to predict the epidemic onset for
vector-borne diseases with decaying vector populations, i.e. the case where
traditional methods fail. It is defined as the average number of infections
produced by an infected individual in one generation Eq. 8,

$\overline{R_{0}}=\left<R_{0}^{i}(t)\right>\Big{|}_{0}^{t_{g}}=R_{0}\left[1-\frac{1}{\tau}\left(f-1\right)\left(e^{-\tau}-1\right)\right]=R_{0}\cdot\mathcal{F}$
(18)

where $\tau=1+\mu/\gamma$ and $\mathcal{F}$ is the expression in brackets
accounting for the effect of the decaying vector population on the stationary
$R_{0}$ (see B for the full derivation of Eq. 18).

Figure 4: Numerical verification of the expression for the basic reproduction
number for vector-borne diseases with decaying vector populations Eq. 18.
Final size of the epidemic as a function of the basic reproduction number in
panels: (a) linear scale; (b) logarithmic scale. Phase space trajectories in
panels: (c) $I_{H}/N_{H}$ vs $I_{v}/N_{v}(0)$ and (d) $I_{H}/N_{H}$ vs
$S_{H}/N_{H}$, where an initial condition
$I_{H}(0)/N_{H}=0.01,S_{H}(0)/N_{H}=0.99$ and $I_{v}(0)/N_{V}(0)=0$ has been
used for the $3$ cases. $\mu=\gamma$ has been used in all the simulations.

A first observation is that $\overline{R_{0}}>R_{0}$ always. This stems from
the fact that $\tau>1$, so that $e^{-\tau}-1<0$, and $f-1>0$, which yields
$\mathcal{F}>1$. This discussion unravels why standard methods fail to predict
the onset of an epidemic under decaying vector populations. Another important
point is that if $\mu/\gamma\gg 1$, which implies $\tau\gg 1$, then
$\overline{R_{0}}\to R_{0}$ as,

$\lim_{\tau\gg 1}\mathcal{F}=\lim_{\tau\gg
1}\left[1-\frac{1}{\tau}\left(f-1\right)\left(e^{-\tau}-1\right)\right]=1+\frac{f-1}{\tau}\
,$ (19)

and if furthermore $\tau\sim\frac{\mu}{\gamma}\gg(f-1)$ then $\mathcal{F}\to
1$ and $\overline{R_{0}}\to R_{0}$. This is in agreement with the discussion
in Section 3.1 showing that the $R_{0}$ computed from standard methods works
if $\mu\gg\gamma$.

Fig. 4(a-b) contrasts numerically the validity of Eq. 18 to predict the final
size of the epidemic as a function of the general basic reproduction number,
$\overline{R_{0}}$, in linear and logarithmic scale, respectively. We observe
that, independently of the initial condition of vectors, the outbreak occurs
for $\overline{R_{0}}>1$. However, we may notice that for large values of the
initial condition of vectors the final size of the epidemic grows more slowly,
so that larger values of $\overline{R_{0}}$ are needed to produce a proper
outbreak. This can be explained by the fact that for $\overline{R_{0}}$
slightly above the threshold, $\overline{R_{0}}=1$, and large values of
$f=N_{v}(0)/N_{v}^{*}$, infections are produced only in the transient period
of the dynamics, as $R_{0}<1$. This is, while the vector population is
decaying to its stationary value, the vectors are able to produce new
infections, but once the vector population reaches the stationary value, the
epidemics stops. This transmission mechanism is radically different to that of
vector-borne diseases with stationary vector populations in which the pre-
pandemic disease-free state is an equilibrium of the system. The phase-space
plots in Fig. 4(c-d) show that the time-averaged basic reproduction number
$\overline{R_{0}}$ is able to accurately predict the conditions under which
the infected host population will grow, in contrast with $R_{0}$ computed in
the post-pandemic fixed point. In essence, for $\overline{R_{0}}>1$ the
infected host population, $I_{H}$, grows before reaching the absorbing state,
$I_{H}=I_{v}=0$, while for $\overline{R_{0}}<1$ the infected host population
is monotonically decreasing. We note that Eq. 18 is similar to the time-
averaged basic reproduction number presented in [Wesley & Allen, 2009] for the
periodic case, which is a first-order approximation to the true basic
reproductive number [Bacaër & Guernaoui, 2006].

###  3.3 Fast-slow approximation

The original $5$-D Eq. 2 model is certainly not amenable to mathematical
analyses due to its high phase-space dimensionality and the fact that it
depends on $4$ parameters. Moreover, in a real-case application, if the
parameters conforming the model are not known the model could suffer from
parameter unidentifiability. However, some approximations can be performed to
reduce the mathematical complexity of the model, as for instance a fast-slow
(or adiabatic) approximation.

If the time-scale of the vector population evolution is faster than that of
the infected hosts, what is expected to be a good approximation in many
practical cases, the vector population will almost instantaneously adapt to
its stationary value. Thus, if $1/\mu\ll 1/\gamma$, or equivalently if
$\mu\gg\gamma$, we can rewrite the time derivative of the vector infected
population as

$\epsilon\dot{I}_{v}=\frac{\alpha}{\mu}S_{v}\frac{I_{H}}{N_{H}}-I_{v}\ ,$ (20)

with $\epsilon=1/\mu$ being a small parameter. Then, $\dot{I_{v}}$ can be
neglected and the infected vector population can be obtained from the
relationship,

$I_{v}\approx\frac{\alpha}{\mu}\frac{S_{v}I_{H}}{N_{H}}\ .$ (21)

Substituting Eq. 21 into the original system Eq. 2 and the identity
$N_{v}(t)=S_{v}(t)+I_{v}(t)$ we obtain the following reduced system,

$\displaystyle\dot{S}_{H}$ $\displaystyle=-\frac{\beta\alpha N_{v}(t)}{\mu
N_{H}+\alpha I_{H}}\frac{S_{H}I_{H}}{N_{H}}$ (22) $\displaystyle\dot{I}_{H}$
$\displaystyle=\frac{\beta\alpha N_{v}(t)}{\mu N_{H}+\alpha
I_{H}}\frac{S_{H}I_{H}}{N_{H}}-\gamma I_{H}$ $\displaystyle\dot{R}_{H}$
$\displaystyle=\gamma I_{H}\ ,$

which can be understood as a SIR model with time-dependent coefficients,

$\displaystyle\dot{S}_{H}$ $\displaystyle=-\xi(t)\frac{S_{H}I_{H}}{N_{H}}$
(23) $\displaystyle\dot{I}_{H}$
$\displaystyle=\xi(t)\frac{S_{H}I_{H}}{N_{H}}-\gamma I_{H}$
$\displaystyle\dot{R}_{H}$ $\displaystyle=\gamma I_{H}\ .$

A SIR model with time-dependent periodic coefficients was discussed in [Bacaër
& Gomes, 2009].

Moreover, if $f\neq 1$ the above mentioned timescales relationship must fulfil
$\displaystyle\frac{\mu}{\gamma}\gg\absolutevalue{\ln{\frac{\epsilon}{f-1}}}$
(cf. Eq. 11) and not only $\displaystyle\frac{\mu}{\gamma}\gg 1$. It is
important to notice that the presence of horizontal transmission would simply
rescale the coefficient $\xi(t)$, and the SIR reduction Eq. 23 would keep its
validity.

In Fig. 5 we numerically verify the validity of the presented fast-slow
approximation. As expected, we observe that the approximation breaks down for
$\mu\sim\gamma$ (Fig. 5(a)), while as $\mu$ becomes larger than $\gamma$ the
approximation improves Fig. 5(b) and it becomes quantitative when
$\mu\gg\gamma$, Fig. 5(c). Finally, we show in Fig. 5(d) a comparison between
the dynamics of the hosts using both the original and the approximated model
using the same parameters than in Fig. 5(c), where the results of both models
are expected to converge.

Figure 5: Numerical verification of the time-scale approximation (Eq. 21) with
$N_{H}=100$, $\alpha=\gamma=1$. $\beta$ is chosen such that $R_{0}=3$. (a)
$\mu=1$, (b) $\mu=10$, (c) $\mu=100$. Panel (d) shows a comparison between the
approximate and original models for the parameters used in (c), where the
approximated models is expected to represent well the original one.

###  3.4 Reduction to a SIR model

The conditions for which the time-scale approximation is valid,
$\mu\gg\gamma$, imply that the vector population will reach its stationary
value almost instantaneously, so that $N_{v}(t)\approx N_{v}^{*}$. Thus, the
system in Eq. 22 can be written approximately as the following SIR-like model,

$\displaystyle\dot{S}_{H}$
$\displaystyle=-\beta^{\prime}\frac{S_{H}I_{H}}{\lambda N_{H}+I_{H}}$ (24)
$\displaystyle\dot{I}_{H}$
$\displaystyle=\beta^{\prime}\frac{S_{H}I_{H}}{\lambda N_{H}+I_{H}}-\gamma
I_{H}$ $\displaystyle\dot{R}_{H}$ $\displaystyle=\gamma I_{H}\ ,$

where $\beta^{\prime}=\beta N_{v}^{*}/N_{H}$ and $\lambda=\mu/\alpha$.

Furthermore, if $\lambda N_{H}\gg I_{H}$ (which is indeed plausible in this
limit) Eq. 24 the model can be written as a SIR model with constant
coefficients,

$\displaystyle\dot{S}_{H}$
$\displaystyle=-\beta_{eff}\frac{S_{H}I_{H}}{N_{H}}$ (25)
$\displaystyle\dot{I}_{H}$
$\displaystyle=\beta_{eff}\frac{S_{H}I_{H}}{N_{H}}-\gamma I_{H}$
$\displaystyle\dot{R}_{H}$ $\displaystyle=\gamma I_{H}\ ,$

where
$\displaystyle\beta_{eff}=\frac{\beta^{\prime}}{\lambda}=\frac{\beta\alpha
N_{v}^{*}}{\mu N_{H}}$.

Figure 6: Comparison between the original model and the reductions, Eq. 24
(SIR-like) and Eq. 25 (SIR) with $N=100$, $\mu/\gamma=10^{3}$ and $f=1$.
$\beta$ was chosen such that $R_{0}=3$. (a) $\lambda=1$, (b)
$\lambda=10^{-3}$, (c) Mean Squared Error between the original model and the
SIR approximations as function of the ratio $\mu/\gamma$ and $f$.

In Fig. 6 we show the validity of the reduced models Eq. 24 and Eq. 25. Fig.
6(a) shows that the SIR-like model (Eq. 24) works when the time-scale
approximation can be performed (as $\mu/\gamma\gg 1$) but the SIR model fails
when the condition $\lambda N_{H}\gg I_{H}$ is not fulfilled. Conversely, in
Fig. 6(b) we show that as $\lambda N_{H}\gg I_{H}$ is fulfilled, then the SIR
model perfectly matches the original model. Finally, Fig. 6(c) shows the
decrease in the mean squared error of the approximation as the condition Eq.
11 is fulfilled for different values of $f$.

##  4 Conclusions

In the present work we have analysed several features of a compartmental
deterministic model for vector-borne diseases with $3$ compartments for hosts
and $2$ for vectors, that does not consider neither horizontal nor vertical
transmission. The focus is to study the behaviour of the model in the case
that the vector population is not stationary. In this case, the pre-pandemic
disease-free state is not a fixed point (equilibrium state) of the dynamical
system, and, in principle, the methods that are customarily used to determine
the basic reproduction number, $R_{0}$ do not work. This is so because these
methods determine the onset of an outbreak by performing a linear stability
analysis of the disease-free state, assuming that it is a fixed point of the
model. A common assumption made in the literature is to determine $R_{0}$ from
the asymptotic state for the vectors (if it is not an extinction state).

We have analysed several initial conditions of the vector population,
characterising different regimes. In the case that the initial condition for
the number of vectors is below the asymptotic state, implying that the vector
population overall grows, then $R_{0}$ as determined from the asymptotic state
correctly predicts the existence (or not) of an epidemic outbreak, but with a
temporal delay in its appearance. This result contrasts with the situation in
which the initial state is above the asymptotic state, with an overall
decrease in the vector population. In this case $R_{0}$ determined from the
asymptotic state may fail badly, predicting no outbreak while a large fraction
of the population might get infected. We present a simple, albeit useful,
generalisation of $R_{0}$ that is able to give a reasonable prediction of the
epidemic threshold for decaying populations, including the case in which
vectors become extinct, a case in which the asymptotic estimation to determine
$R_{0}$ cannot be applied.

Compartmental models of vector-borne diseases usually have many compartments
and parameters, which can lead to a problem of parameter unidentifiability.
The model analysed here is not an exception, and when applied to real-world
cases many different combinations of the parameters could be able to reproduce
the available data. Thus, in order to facilitate the application of the model
to experimental data, we have studied a useful fast-slow (or adiabatic)
approximation that allows to reduce the model if the parameters fulfil certain
conditions. In particular, our study shows that under quite realistic
assumptions (the typical timescale of hosts infection and death is much slower
than vector timescales) it is possible to obtain a reduced SIR model. We
recall that this reduction implies that, under these assumptions, the process
by which hosts (that could be immobile) get infected through the action of
vectors is equivalent to a direct interaction among hosts.

The deterministic compartmental model analysed here, with some modifications,
is a clear candidate to study many vector-borne diseases, in particular
phytopathologies. Furthermore, in case of parameter unidientifiability the
model reductions performed in this work could be useful to solve this issue.
In any case, this description is still idealised, as compartmental models
imply a well-mixed assumption in which space is not explicitly described. This
kind of representations are not always applicable to real-world scenarios
although are useful as a first approximation. Thus, future research should
focus on the integration of space and vector mobility in the model to account
for more realistic situations.

## Acknowledgments

AGR and MAM acknowledge financial support from Grant RTI2018-095441-B-C22
(SuMaEco) funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of
making Europe” and from the María de Maeztu Program for Units of Excellence in
R&D (No. MDM-2017-0711).

##  Appendix A Calculation of $R_{0}$ from standard methods

The standard methods of calculation of $R_{0}$ are based in the linear
stability analysis of the disease-free equilibrium, either directly, through
the linear analysis of the fixed point, that yields the stability condition
from which $R_{0}$ can be obtained, or using the Next Generation Method (NGM)
[Diekmann et al., 2010] that provides directly $R_{0}$ by solving a suitable
linear problem. Customarily these methods are applied to a pre-pandemic
disease-free equilibrium, but as there is no such state in the case of non-
stationary populations, here a similar approach is applied to a post-pandemic
or asymptotic disease-free equilibrium.

### Linear stability analysis

In order to perform the linear stability analysis of the fixed point
($I_{H}=I_{v}=0$) we first need to compute the Jacobian matrix, $J$,

$J=\begin{pmatrix}-\beta\frac{I_{v}}{N_{H}}&0&0&-\beta\frac{S_{H}}{N_{H}}\\\
\beta\frac{I_{v}}{N_{H}}&-\gamma&0&\beta\frac{S_{H}}{N_{H}}\\\
0&-\alpha\frac{S_{v}}{N_{H}}&-\alpha\frac{I_{H}}{N_{H}}-\mu&0\\\
0&\alpha\frac{S_{v}}{N_{H}}&\alpha\frac{I_{H}}{N_{H}}&-\mu\\\ \end{pmatrix}$
(26)

Then, we evaluate the Jacobian at the fixed point (or disease free
equilibrium, DFE), yielding

$J|_{DFE}=\begin{pmatrix}0&0&0&-\beta\\\ 0&-\gamma&0&\beta\\\
0&-\alpha\frac{C}{N_{H}}\frac{\delta}{\mu}&-\mu&0\\\
0&\alpha\frac{C}{N_{H}}\frac{\delta}{\mu}&0&-\mu\\\ \end{pmatrix}$ (27)

where $S_{H}=N_{H}$ has been considered.

Finally, we obtain the eigenvalues of this matrix as

$\begin{split}det(J|_{DFE}-\lambda\mathbb{I})&=-\lambda\bigg{[}-(\mu+\lambda)^{2}(\gamma+\lambda)+(\mu+\lambda)\beta\alpha\frac{C}{N_{H}}\frac{\delta}{\mu}\bigg{]}=0\Rightarrow\\\
\lambda_{0}&=0\\\ \lambda_{\mu}&=-\mu\\\
\lambda_{\pm}&=-\frac{(\gamma+\mu)}{2}\pm\frac{1}{2}\sqrt{(\gamma-\mu)^{2}+4\beta\alpha\frac{C}{N_{H}}\frac{\delta}{\mu}}\end{split}$
(28)

It is straightforward to see that all eigenvalues are real and the stability
of the disease-free equilibrium is determined by the sign of the eigenvalues.
$\lambda_{\mu}=-\mu<0$ as $\mu$ is defined positive, so in order to discuss
the stability of this fixed point, we need to study the $\lambda_{\pm}$
eigenvalues. $\lambda_{-}$ is always negative, but $\lambda_{+}$ changes sign
depending on the values of the parameters. The threshold condition
$\lambda_{+}=0$ leads to:

$\lambda_{+}=0\;\Rightarrow\;\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}=1$
(29)

So, for
$\displaystyle\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}<1\,\Rightarrow\,\lambda_{+}<0$
the fixed point is stable and for
$\displaystyle\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}>1\,\Rightarrow\,\lambda_{+}>0$
a perturbation will grow in the direction of the eigenvector associated to
$\lambda_{+}$. Thus, this threshold defines the basic reproduction number,

$R_{0}=\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}$ (30)

If instead of $S_{H}=N_{H}$ one considers any initial condition of hosts,
$S_{H}(0)$, the basic reproduction number is given by,

$R_{0}=\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}\frac{S_{H}(0)}{N_{H}}$
(31)

### Next Generation Matrix method

The previous result can also be obtained by means of the NGM method, which is
explained in detail in[Diekmann et al., 2010]. Basically the method is based
in decomposing the Jacobian in the form
$\textbf{J}=\mathbf{T}+\mathbf{\Sigma}$, where $\mathbf{T}$ is the
transmission part, that describes the production of new infections, and
$\mathbf{\Sigma}$ the transition part, that describes changes of state
(including death). Then, it can be proved [Diekmann et al., 2010] that the
basic reproduction number $R_{0}$ is given by the spectral radius (i.e. the
largest eigenvalue) of the (next generation) matrix
$\mathbf{K}=-\mathbf{T}\mathbf{\Sigma}^{-1}$.

$\mathbf{K}=-\mathbf{T}\mathbf{\Sigma}^{-1}=\begin{pmatrix}\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}&\frac{\beta}{\mu}\\\
\ ,0&0\\\ \end{pmatrix}$ (32)

with,

$\mathbf{T}=\begin{pmatrix}0&\beta\frac{N_{H}}{N_{H}}\\\
0&0\end{pmatrix}\qquad\mathbf{\Sigma}=\begin{pmatrix}-\gamma&0\\\
\alpha\frac{C}{N_{H}}\frac{\delta}{\mu}&-\mu\end{pmatrix}\;\Rightarrow\;-{\mathbf{\Sigma}}^{-1}=\begin{pmatrix}\frac{1}{\gamma}&0\\\
\frac{\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}&\frac{1}{\mu}\end{pmatrix}$
(33)

The basic reproduction number is the spectral radius of this matrix so:

$det(\mathbf{K}-\sigma\mathbb{I})=0\implies\begin{vmatrix}\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}-\sigma&\frac{\beta}{\mu}\\\
0&-\sigma\\\
\end{vmatrix}=(-\sigma)\bigg{(}\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}-\sigma\bigg{)}=0\implies$
$\sigma=\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}};\quad\sigma=0$ (34)

Therefore, the basic reproduction number is

$R_{0}=\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}$ (35)

If instead of $S_{H}=N_{H}$ one considers any initial condition of hosts,
$S_{H}(0)$, the basic reproduction number is given by,

$R_{0}=\frac{\beta\alpha}{\gamma\mu}\frac{C}{N_{H}}\frac{\delta}{\mu}\frac{S_{H}(0)}{N_{H}}$
(36)

##  Appendix B Calculation of $R_{0}$ for non-stationary vector populations

A simple way of extending the definition of $R_{0}$ in the case of non-
stationary vector populations consists of averaging the number of secondary
infections produced by an infected individual along one generation, that is
equivalent to averaging the instantaneous definition of $R_{0}$, namely
$R_{0}^{i}$, over one generation,

$\overline{R_{0}}=\left<R_{0}^{i}(t)\right>\Big{|}_{0}^{t_{g}}=\frac{R_{0}}{N_{v}^{*}}\left<N_{v}(t)\right>\Big{|}_{0}^{t_{g}}=\frac{R_{0}}{N_{v}^{*}}\frac{1}{t_{g}}\int_{0}^{t_{g}}N_{v}(t)\,\mathrm{d}t\
,$ (37)

where the integral in Eq. 37 is solved as

$\int_{0}^{t_{g}}N_{v}(t)\,\mathrm{d}t=\left[N_{v}^{*}t-\frac{1}{\mu}\left(N_{v}(0)-N_{v}^{*}\right)e^{-\mu
t}\right]_{0}^{t_{g}}=N_{v}^{*}t_{g}-\frac{1}{\mu}\left(N_{v}(0)-N_{v}^{*}\right)\left[e^{-\mu
t_{g}}-1\right]\ .$ (38)

Thus, the basic reproduction number for non-stationary vector populations is
given by

$\overline{R_{0}}=\frac{R_{0}}{N_{v}^{*}}\left\\{N_{v}^{*}-\frac{1}{\mu
t_{g}}\left[N_{v}(0)-N_{v}^{*}\right]\left[e^{-\mu t_{g}}-1\right]\right\\}\
,$ (39)

where the generation time, $t_{g}$, is Eq. 8. Eq. 39 can be rewritten as,

$\overline{R_{0}}=\left<R_{0}^{i}(t)\right>\Big{|}_{0}^{t_{g}}=R_{0}\left[1-\frac{1}{\tau}\left(f-1\right)\left(e^{-\tau}-1\right)\right]=R_{0}\cdot\mathcal{F}\
,$ (40)

where $\tau=1+\mu/\gamma$ and $\mathcal{F}$ is the expression in brackets,
which accounts for the effect of the decaying vector population on the
stationary $R_{0}$.

In our approach, a generation is defined as the time elapsed in the following
sequence of processes: 1) A host individual becomes infected; 2) The infected
host passes the disease to a susceptible vector; 3) The infected vector dies.
Basically, the time elapsed from the first to the last process is the time in
which new infections can be produced, i.e. $t_{g}$ Eq. 8.

## References

  * Anderson [1991] Anderson, R. M. (1991).  Discussion: The Kermack-McKendrick epidemic threshold theorem.  Bulletin of Mathematical Biology, 53, 3–32. doi:10.1016/S0092-8240(05)80039-4.
  * Athni et al. [2020] Athni, T. et al. (2020).  How vector-borne disease shaped the course of human history.  doi:10.22541/au.159135318.82566392.
  * Bacaër [2007] Bacaër, N. (2007).  Approximation of the basic reproduction number r0 for vector-borne diseases with a periodic vector population.  Bulletin of Mathematical Biology, 69, 1067–1091. doi:10.1007/s11538-006-9166-9.
  * Bacaër & Gomes [2009] Bacaër, N., & Gomes, M. G. M. (2009).  On the final size of epidemics with seasonality.  Bulletin of Mathematical Biology, 71, 1954. doi:10.1007/s11538-009-9433-7.
  * Bacaër & Guernaoui [2006] Bacaër, N., & Guernaoui, S. (2006).  The epidemic threshold of vector-borne diseases with seasonality.  Journal of Mathematical Biology, 53, 421–436. doi:10.1007/s00285-006-0015-0.
  * Bragard et al. [2013] Bragard, C., Caciagli, P., Lemaire, O., Lopez-Moya, J., MacFarlane, S., Peters, D., Susi, P., & Torrance, L. (2013).  Status and prospects of plant virus control through interference with vector transmission.  Annual Review of Phytopathology, 51, 177–201. doi:10.1146/annurev-phyto-082712-102346.
  * Brauer [2008] Brauer, F. (2008).  Compartmental models in epidemiology.  In F. Brauer, P. van den Driessche, & J. Wu (Eds.), Mathematical Epidemiology (pp. 19–79).  Berlin, Heidelberg: Springer Berlin Heidelberg.  doi:10.1007/978-3-540-78911-6_2.
  * Brauer et al. [2016] Brauer, F., Castillo-Chavez, C., Mubayi, A., & Towers, S. (2016).  Some models for epidemics of vector-transmitted diseases.  Infectious Disease Modelling, 1, 79–87. doi:10.1016/j.idm.2016.08.001.
  * Castillo-Chavez & Thieme [1995] Castillo-Chavez, C., & Thieme, H. R. (1995).  Asymptotically autonomous epidemic models.  In D. Arino, O. Axelrod, M. Kimmel, & M. Langlais (Eds.), Mathematical Population Dynamics: Analysis of Heterogeneity, Vol. I, Theory of Epidemics (pp. 33–50).  Wuerz publishing (Winnipeg).  URL: https://ecommons.cornell.edu/bitstream/handle/1813/31834/BU-1248-M.pdf.
  * Chowell [2017] Chowell, G. (2017).  Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts.  Infectious Disease Modelling, 2, 379–398. doi:10.1016/j.idm.2017.08.001.
  * Diekmann & Heesterbeek [2000] Diekmann, O., & Heesterbeek, J. A. P. (2000).  Mathematical Epidemiology of Infectious Diseases. Model Building, Analysis and Interpretation.  Chichester (UK): John Wiley and Sons.
  * Diekmann et al. [2010] Diekmann, O., Heesterbeek, J. A. P., & Roberts, M. G. (2010).  The construction of next-generation matrices for compartmental epidemic models.  Journal of The Royal Society Interface, 7, 873–885. doi:10.1098/rsif.2009.0386.
  * Esteva & Vargas [1998] Esteva, L., & Vargas, C. (1998).  Analysis of a dengue disease transmission model.  Mathematical Biosciences, 150, 131–151. doi:10.1016/S0025-5564(98)10003-2.
  * Garms et al. [1979] Garms, R., Walsh, J., & Davies, J. (1979).  Studies on the reinvasion of the onchocerciasis control programme in the volta river basin by simulium damnosum si with emphasis on the south-western areas.  Tropenmedizin und Parasitologie, 30, 345–362.
  * Giménez-Romero et al. [2021] Giménez-Romero, A., Grau, A., Hendriks, I. E., & Matias, M. A. (2021).  Modelling parasite-produced marine diseases: The case of the mass mortality event of pinna nobilis.  Ecological Modelling, 459, 109705. doi:10.1016/j.ecolmodel.2021.109705.
  * Huang et al. [2020] Huang, W., Reyes-Caldas, P., Mann, M., Seifbarghi, S., Kahn, A., Almeida, R. P., Béven, L., Heck, M., Hogenhout, S. A., & Coaker, G. (2020).  Bacterial vector-borne plant diseases: Unanswered questions and future directions.  Molecular Plant, 13, 1379–1393. doi:10.1016/j.molp.2020.08.010.
  * Kamgang & Sallet [2008] Kamgang, J. C., & Sallet, G. (2008).  Computation of threshold conditions for epidemiological models and global stability of the disease-free equilibrium (dfe).  Mathematical Biosciences, 213, 1–12. doi:10.1016/j.mbs.2008.02.005.
  * Kao & Eisenberg [2018] Kao, Y.-H., & Eisenberg, M. C. (2018).  Practical unidentifiability of a simple vector-borne disease model: Implications for parameter estimation and intervention assessment.  Epidemics, 25, 89–100. doi:10.1016/j.epidem.2018.05.010.
  * Kermack & McKendrick [1927] Kermack, W. O., & McKendrick, A. G. (1927).  A contribution to the mathematical theory of epidemics.  Proceedings of the Royal Society of London. Series A, 115, 700–721. doi:10.1098/rspa.1927.0118.
  * Lashari & Zaman [2011] Lashari, A. A., & Zaman, G. (2011).  Global dynamics of vector-borne diseases with horizontal transmission in host population.  Computers & Mathematics with Applications, 61, 745–754. doi:10.1016/j.camwa.2010.12.018.
  * Laukó [2006] Laukó, I. G. (2006).  Stability of disease free sets in epidemic models.  Mathematical and Computer Modelling, 43, 1357–1366. doi:10.1016/j.mcm.2005.06.011.
  * Macdonald [1957] Macdonald, G. (1957).  The epidemiology and control of malaria.  Oxford (UK): Oxford University Press.
  * Martcheva [2015] Martcheva, M. (2015).  An Introduction to Mathematical Epidemiology.  New York: Springer.  doi:10.1007/978-1-4899-7612-3.
  * Rocklöv & Dubrow [2020] Rocklöv, J., & Dubrow, R. (2020).  Climate change: an enduring challenge for vector-borne disease prevention and control.  Nature Immunology, 21, 479–483. doi:10.1038/s41590-020-0648-y.
  * Roosa & Chowell [2019] Roosa, K., & Chowell, G. (2019).  Assessing parameter identifiability in compartmental dynamic models using a computational approach: application to infectious disease transmission models.  Theoretical Biology and Medical Modelling, 16, 1. doi:10.1186/s12976-018-0097-6.
  * Rybicki [2015] Rybicki, E. P. (2015).  A top ten list for economically important plant viruses.  Archives of Virology, 160, 17–20. doi:10.1007/s00705-014-2295-9.
  * Schneider et al. [2020] Schneider, K., van der Werf, W., Cendoya, M., Mourits, M., Navas-Cortés, J. A., Vicent, A., & Oude Lansink, A. (2020).  Impact of xylella fastidiosa subspecies pauca in european olives.  Proceedings of the National Academy of Sciences, 117, 9250–9259. doi:10.1073/pnas.1912206117.
  * Schumacher & Campbell [2018] Schumacher, S. K., & Campbell, J. I. (2018).  Chapter 56 - travel medicine.  In R. P. Olympia, R. M. O’Neill, & M. L. Silvis (Eds.), Urgent Care Medicine Secrets (pp. 352–357).  Elsevier.  doi:10.1016/B978-0-323-46215-0.00056-2.
  * Shah & Jyoti [2013] Shah, N., & Jyoti, G. (2013).  SEIR model and simulation for vector borne diseases.  Applied Mathematics, 4, 13 – 17. doi:10.4236/am.2013.48A003.
  * Thieme [1992] Thieme, H. R. (1992).  Convergence results and a Poincaré-Bendixson trichotomy for asymptotically autonomous differential equations.  Journal of Mathematical Biology, 30, 755–763. doi:10.1007/BF00173267.
  * Tumber et al. [2014] Tumber, K., Alston, J., Fuller, K. et al. (2014).  Pierce’s disease costs California $104 million per year.  California Agriculture, 68, 20–29.
  * van den Driessche [2017] van den Driessche, P. (2017).  Reproduction numbers of infectious disease models.  Infectious Disease Modelling, 2, 288–303. doi:10.1016/j.idm.2017.06.002.
  * Wei et al. [2008] Wei, H.-M., Li, X.-Z., & Martcheva, M. (2008).  An epidemic model of a vector-borne disease with direct transmission and time delay.  Journal of Mathematical Analysis and Applications, 342, 895–908. doi:10.1016/j.jmaa.2007.12.058.
  * Wesley & Allen [2009] Wesley, C. L., & Allen, L. J. (2009).  The basic reproduction number in epidemic models with periodic demographics.  Journal of Biological Dynamics, 3, 116–129. doi:10.1080/17513750802304893.  PMID: 22880824.
  * [35] WHO (2018).  Fact sheets of vector-borne diseases.  URL: https://www.hiobs.org/lib/file/manager/Vector-Borne_Disease_General_Fact_Sheet.pdf.
  * Zhao et al. [2020] Zhao, S., Musa, S. S., Hebert, J. T., Cao, P., Ran, J., Meng, J., He, D., & Qin, J. (2020).  Modelling the effective reproduction number of vector-borne diseases: the yellow fever outbreak in Luanda, Angola 2015-2016 as an example.  PeerJ, 8, e8601–e8601. doi:10.7717/peerj.8601.


# A Universal Railway Obstacle Detection System based on Semi-supervised
Segmentation And Optical Flow

Qiushi Guo
CSRD
<EMAIL_ADDRESS>
###### Abstract

Detecting obstacles in railway scenarios is both crucial and challenging due
to the wide range of obstacle categories and varying ambient conditions such
as weather and light. Given the impossibility of encompassing all obstacle
categories during the training stage, we address this out-of-distribution
(OOD) issue with a semi-supervised segmentation approach guided by optical
flow clues. We reformulate the task as a binary segmentation problem instead
of the traditional object detection approach. To mitigate data shortages, we
generate highly realistic synthetic images using Segment Anything (SAM) and
YOLO, eliminating the need for manual annotation to produce abundant pixel-
level annotations. Additionally, we leverage optical flow as prior knowledge
to train the model effectively. Several experiments are conducted,
demonstrating the feasibility and effectiveness of our approach.

##  1 Introduction

With the rapid advancement of high-speed trains, ensuring the security of
railway systems has emerged as a critical public concern. One of the primary
challenges is obstacle detection, which plays a crucial role in railway
safety. Developing a reliable and scalable obstacle detection system can
empower train operators and dispatchers to take preemptive actions and
mitigate potential accidents.

Deep learning techniques have been widely adopted across various security
domains, including mobile payments[3], disaster detection [9], and fraud
detection [6]. This technology exhibits substantial promise in enhancing
railway safety through sophisticated obstacle detection capabilities.
Significant efforts have recently been devoted to addressing obstacle
detection using deep learning methods. Although these approaches have achieved
some success, they also exhibit notable disadvantages:

Figure 1: The pipeline for synthetic data generation, utilizing SAM and YOLO
to extract target objects from a gallery and superimpose them onto a base
image (specifically, a railway image). Notably, this process does not
necessitate annotations.

  * •

Fragility to complex ambient conditions

  * •

Requirement for extensive manual annotations

  * •

Difficulty in extending to different scenarios

Designing an extendable, annotation-free model with strong generalization
ability remains a significant challenge in both industry and academia.

To address the aforementioned issues, we propose a semi-supervised approach
guided by optical flow. To mitigate the data shortage problem, we employ SAM
[7] and YOLO [8] to generate highly realistic pseudo-images for training.
Instead of manually collecting and annotating images pixel by pixel, we
prepare two image sets: base images (fewer than 100 background images with
only railway areas annotated) and object images. The object images include
categories such as pedestrians, animals, and textures. Using SAM and YOLO, we
obtain masks for the intended objects in these images. These objects are then
pasted onto the base images according to the masks. The entire process are
illustrated as Fig. 1 This process simultaneously generates image and mask
pairs without manual effort.

To address the challenges posed by varying weather conditions, we implement
two complementary strategies. Firstly, we compile a dataset of base images
captured under diverse weather conditions, including rainy, foggy, and clear
(sunny) environments. Secondly, we utilize optical flow to provide positional
information as prior knowledge. For optical flow predictions, we generate
pseudo sequences of obstacles. This involves creating an initial pseudo frame
at point $P_{i}(x,y)$ and subsequently generating a new frame at
$P_{i+1}(x+\delta,y+\delta)$ with the same object superimposed. Experimental
results indicate that our approach yields satisfactory performance across
different weather scenarios.

Our contributions are summarized as follows:

  * •

We reformulate the obstacle detection task as a binary segmentation problem,
distinguishing between railway areas and non-railway areas.

  * •

We introduce a simple yet effective data generation mechanism to synthesize
realistic images using SAM and YOLO.

  * •

Optical flow is leveraged to generate prior knowledge that guides the
segmentation network.

Figure 2: Pipeline of our proposed method.

##  2 Related work

###  2.1 Obstacle Detection in Railyway

Matthias Brucker et al.[2]propose a a shallow netwrok to learn railway
segmentation from normal railway images. They explore the controlled inclusion
of global information by learning to hallucinate obstacle-free images. Zhang
Qiang et al.[12]. conbine segmentation model with the LiDAR in their obstacle
detection system; Amine Boussik et al.[1] propose an unsupervised models based
on a large set of generated convolutional auto-encoder models to detect
obstacles on railway’s track level.

###  2.2 Segmentation with Optical flow

Laura et al. [10]. demonstrate the effectiveness of jointly optimizing optical
flow and video segmentation using an iterative scheme; Volodymyr et al[5].
present an architecture for Video Object Segmentation that combines memory-
based matching with motion-guided propagation resulting in stable long-term
modeling and strong temporal consistency.

##  3 Method

The pipeline of our approach is illustrated in Fig.2. Given a set of base
images $B$ and target images $T$, our objective is to identify potential
obstacles within specific regions $\eta$. Unlike traditional detection methods
that categorically detect each obstacle, we reformulate the problem as a
binary segmentation task. Instead of attempting to detect all potential
obstacles, which is impractical, our emphasis is on segmenting the railway
area, a region that remains consistent over time compared to obstacles.

To simulate these scenarios effectively, we generate highly realistic pseudo-
images using a copy-paste approach. Additionally, to address challenges posed
by extreme weather conditions, which can obscure object segmentation, we
introduce optical flow to provide prior information guiding the segmentation
model. Pseudo images $I_{t}$ and $I_{t+\delta}$ are generated by applying a
small shift $\delta$ to the target object, simulating its movement. The output
of the optical flow model is incorporated along with pseudo images as input to
facilitate accurate predictions. This section will delve into the detailed
methodology employed throughout this process.

###  3.1 Data Acquisition

Base Images are used in our experiments are gathered at our facility in
Chengdu, which features a railway spanning over 60 meters and includes
simulators for fog and rain conditions. To ensure diversity in our dataset, we
capture images under different weather scenarios, specifically rainy, foggy,
and sunny conditions Fig. 3. Due to the fixed position of the camera, only one
mask is required for annotation purposes. Importantly, the railway areas in
the base images are devoid of any potential obstacles. Any obstacles present
are generated using a copy-paste method.

Figure 3: Sample base images depicting various weather conditions. From left
to right, the images illustrate scenes captured under foggy, normal, and rainy
conditions.

Object Image dataset comprises three categories: PennFudanPed, Obj365 (part)
[11], and DTD [4]. To facilitate fully automated application of our
methodology, we proceed under the assumption that no masks are initially
available. We focus on selecting categories likely to occur in our scenario,
such as animals (e.g., deer, horse, cow) and vehicles (e.g., truck, cart).
This ensures our approach is tailored to handle relevant objects effectively.

Figure 4: Sample of generated image-mask pairs. From left to right:
pedestrian, animal and texture

The entire process can be delineated into sequential steps: Initially, object
images are fed into the YOLO model, which returns a list of bounding boxes
identifying detected targets. These bounding boxes serve as inputs for SAM,
which generates segmentation masks to outline the object pixels. Subsequently,
these segmented object pixels are integrated into the base images based on the
segmentation mask guidance. Here, we elaborate on the detailed methodology.

  * •

Object Detection with YOLO: Object images are inputted into the YOLO model,
specifically trained on Obj365, to detect objects belonging to predefined
target categories fitting our scenario.

  * •

Segmentation with SAM: Bounding boxes from YOLO are used as prompts for SAM to
generate segmentation masks. These masks delineate object pixels, facilitating
their extraction from the object images.

  * •

Integration with Base Images: Extracted object pixels are seamlessly
integrated into the corresponding regions of base images, aligning with the
guidance provided by the segmentation masks.

During the SAM stage, while not every segmentation mask achieves perfection,
each contributes to the overall objective of accurately segmenting the railway
area rather than focusing on obstacles. To address challenges related to out-
of-distribution (OOD) scenarios, we introduce random polygon generation with
texture rendering from DTD. Additionally, object resizing and rescaling are
applied to enrich image content and bolster model robustness. The rescale
follow the equation below:

$h=\alpha\ast y+\beta$ (1) $w=\frac{h}{H}\ast W$ (2)

where h,w and H,W are the shapes of the target obj and original obj,
respectively. $\alpha$ and $\beta$ are hyper-parameters to adjust the scale.
In our project, we set $\alpha$ to 0.6 and $\beta$ to 30. The value should be
varied by the camera’s position and it’s parameters. The final generated
samples are demonstrated as Fig .4

###  3.2 Optical-Flow

Optical flow is based on the assumption that the intensity of a point in an
image remains constant as it moves from one frame to the next.

$I(x,y,t)=I(x+\Delta t,y+\Delta t,t+\Delta t)$ (3)

In our scenario, we employ RAFT (Recurrent All-Pairs Field Transforms) as our
chosen model, which demonstrates robust performance across a wide range of
scales from tiny to large. The size of obstacles in our dataset varies,
spanning from hundreds of pixels down to less than 50 pixels in size.
Utilizing the RAFT model requires two consecutive frames for optical flow
estimation. Accordingly, we generate two pseudo images $I_{t}$ and $I_{t+1}$,
where the same target objects are pasted with a slight positional shift
$\eta$.

$Motion=\phi(I_{t},I_{t+1})$ (4)

$I_{t+1}=I_{t}(obj_{x}+\Delta x,obj_{y}+\Delta y)$ (5)

We set $\Delta x$ and $\Delta y$ range between 5-10. The motion prediction
will be leveraged as prior information fused with pseudo image to train the
model.

##  4 Experiments

###  4.1 Dataset and Evaluation Metrics

Dataset Our training dataset is consist of three parts: $obs\\_{person}$,
$obs\\_{animals}$ and $obs\\_{textures}$, namely person obstacles, animal
obstacles and obstacles generated from texture polygons. The details are
described as follow: As for test dataset, we recollect images with various
obstacles under different weather conditions in different distance to the
camera.

Metrics $mIoU$ is used to evaluate the performance of our model. $mIoU$ refers
to the Mean Intersection over union, which is a widely used metric in
segmentation task. It can be calculated as follow:

$IoU_{i}=\frac{TP_{i}}{TP_{i}+FP_{i}+FN_{i}}$ (6)
$mIoU=\frac{1}{n}\sum_{i=0}^{n}IoU_{i}$ (7)

pixel accuracy is also a metric to evaluate the segmentation models.

$Pixel\\_{accuracy}=\frac{N\\_{corr}}{N\\_{total}}$ (8)

where $N\\_{corr}$ is the number of correctly classified pixels, $N\\_{total}$
is the number of total pixels.

Name | Volume | Dis(m) | Category
---|---|---|---
$obs\\_{person}$ | 4000 | 0-70 | person
$obs\\_{animal}$ | 4000 | 0-70 | cow,horse,deer
$obs\\_{texture}$ | 2000 | 0-70 | see DTD
$val\\_{near}$ | 200 | 0-20 | person,rock,board
$val\\_{middle}$ | 200 | 20-50 | person,rock,board
$val\\_{far}$ | 200 | 50-70 | person,rock,board
Table 1: Datasets description.

###  4.2 Implementation Details

Our method is implemented using the PyTorch framework and the model is trained
on an RTX 4070Ti. We select Jaccard loss as the loss function and AdamW as the
optimizer. The batch size is set to 8 and the number of epochs to 20. Data
transformations include horizontal flip, coarse dropout, and random brightness
contrast adjustments.

###  4.3 Results

To validate the performance of our approach, we conduct experiments on our
three self-collected datasets: val_near, val_mid, and val_far. The details are
described in Table 1. The basic training dataset contains 10,000 images
(4,000+4,000+2,000). To fully assess the impact of the number of generated
images, we increase the dataset size by 10% and 50% in rows 4 and 5.

The results are illustrated in Table2, which show that both RAFT and
segmentation-based approaches can effectively segment obstacles in our railway
area experiments. Combining RAFT and pseudo-images enhances model performance.
As more generated images are added to the training dataset, the model’s
performance gradually reaches its limit.

| val_near | val_mid | val_far
---|---|---|---
yolov5 | 0.744 | 0.623 | 0.457
Raft | 0.717 | 0.636 | 0.585
DeepLab | 0.825 | 0.817 | 0.747
DeepLab+Raft | 0.843 | 0.828 | 0.709
DeepLab+Raft+10% | 0.837 | 0.843 | 0.724
DeepLab+Raft+50% | 0.863 | 0.851 | 0.802
Table 2: Experiments results.

###  4.4 Ablation Study

We conduct ablation experiment to validate the effect of different target
objects. The results are demonstrated as Table 3 Comparing the row 1,2,3 with
row 4, we can find that each obs dataset contributes to improving the
robustness and accuracy of the model.

| obs_person | obs_animal | obs_texture | mIoU
---|---|---|---|---
1 | ✓ | ✓ | ✗ | 0.781
2 | ✓ | ✗ | ✓ | 0.817
3 | ✗ | ✓ | ✓ | 0.732
4 | ✓ | ✓ | ✓ | 0.849
Table 3: Ablation study.

##  5 Conclusion

This paper introduces a universal segmentation model based on a semi-
supervised approach. To address out-of-distribution (OOD) challenges, we
generate highly realistic pseudo images instead of relying on manual pixel-
level annotations. Additionally, we enhance performance by incorporating
optical flow techniques. Experimental results demonstrate satisfactory
performance across various potential objects.

## References

  * [1] Amine Boussik, Wael Ben-Messaoud, Smail Niar, and Abdelmalik Taleb-Ahmed.  Railway obstacle detection using unsupervised learning: An exploratory study.  In 2021 IEEE Intelligent Vehicles Symposium (IV), pages 660–667. IEEE, 2021.
  * [2] Matthias Brucker, Andrei Cramariuc, Cornelius Von Einem, Roland Siegwart, and Cesar Cadena.  Local and global information in obstacle detection on railway tracks.  In 2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 9049–9056. IEEE, 2023.
  * [3] Han Cai, Ji Lin, Yujun Lin, Zhijian Liu, Haotian Tang, Hanrui Wang, Ligeng Zhu, and Song Han.  Enable deep learning on mobile devices: Methods, systems, and applications.  ACM Transactions on Design Automation of Electronic Systems (TODAES), 27(3):1–50, 2022.
  * [4] Mircea Cimpoi, Subhransu Maji, Iasonas Kokkinos, Sammy Mohamed, and Andrea Vedaldi.  Describing textures in the wild.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3606–3613, 2014.
  * [5] Volodymyr Fedynyak, Yaroslav Romanus, Bohdan Hlovatskyi, Bohdan Sydor, Oles Dobosevych, Igor Babin, and Roman Riazantsev.  Devos: Flow-guided deformable transformer for video object segmentation.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 240–249, 2024.
  * [6] Qiushi Guo, Yifan Chen, and Shisha Liao.  Enhancing mobile privacy and security: A face skin patch-based anti-spoofing approach.  In 2023 IEEE International Conference on Cloud Computing Technology and Science (CloudCom), pages 52–57. IEEE, 2023.
  * [7] Alexander Kirillov, Eric Mintun, Nikhila Ravi, Hanzi Mao, Chloe Rolland, Laura Gustafson, Tete Xiao, Spencer Whitehead, Alexander C Berg, Wan-Yen Lo, et al.  Segment anything.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 4015–4026, 2023.
  * [8] Joseph Redmon, Santosh Divvala, Ross Girshick, and Ali Farhadi.  You only look once: Unified, real-time object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 779–788, 2016.
  * [9] Cem Sazara, Mecit Cetin, and Khan M Iftekharuddin.  Detecting floodwater on roadways from image data with handcrafted features and deep transfer learning.  In 2019 IEEE intelligent transportation systems conference (ITSC), pages 804–809. IEEE, 2019.
  * [10] Laura Sevilla-Lara, Deqing Sun, Varun Jampani, and Michael J Black.  Optical flow with semantic segmentation and localized layers.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3889–3898, 2016.
  * [11] Shuai Shao, Zeming Li, Tianyuan Zhang, Chao Peng, Gang Yu, Xiangyu Zhang, Jing Li, and Jian Sun.  Objects365: A large-scale, high-quality dataset for object detection.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 8430–8439, 2019.
  * [12] Qiang Zhang, Fei Yan, Weina Song, Rui Wang, and Gen Li.  Automatic obstacle detection method for the train based on deep learning.  Sustainability, 15(2):1184, 2023.