On the role of vector modeling in a minimalistic epidemic model

1.   Introduction

The spread of the vector-borne dengue fever in a human population has been studied in many papers: we mention a family of three nested multi-strain compartmental SIR type models analyzed in ^[1,2,3,4]. In ^[3] a realistic description of recorded disease incidents was produced. However, in these models the vector dynamics was taken into account by assuming that the size of the infected vector population is proportional to the infected human-host population.

On the other hand, the analysis of ^[5] shows that an increased degree of personal protection and vector control are effective control measures in dengue-endemic areas. This leads to the conclusion that in order to study the effects of vector-mosquito control an accurate vector population model together with a human-host population model describing the multiple-strain virus outbreaks and spread is required.

In ^[6,7,8] the epidemics of the vectors is modeled explicitly by adding to a model similar to those in ^[1,2,3] extra compartments for the vector epidemics without taking ecological processes into account. The dynamics of the epidemics of the vector is directly related to ecological processes during the very complex life-cycle with four stages: egg (embryo), larva, pupa (embryo) (in the water phase) and imago (adult). Furthermore the numbers of healthy and infected vector individuals will change due to epidemiological processes such as the spread of the disease in the host. The combination of the vector and the host model results then in an eco-epidemiology model where two scientific fields are combined to study ecological and epidemiological issues simultaneously.

The formulation and analyses of such a complicated model and its analysis based almost completely on numerical simulations is beyond the scope of this paper. In the vector-host epidemic model analyzed in this paper we couple a simplified dengue-fever model for the host population based on ^[3] with a simplified model for the mosquito vector population ^[9]. This is done in order to gain insight into the effects of the vector dynamics on the host dynamics, especially when the time-scale of the vector population is much faster than that of the host. In ^[9] it was investigated how far the fast time scale of the mosquito epidemiology can be slaved by the slower human epidemiology, so that for the understanding of human disease data mainly the dynamics of the human time scale is essential because it will be only slightly perturbed by the mosquito dynamics. Useful and realistic parameter values for the dengue fever case from ^[3] were used. It was concluded that owing to the faster vector dynamics it is slaved by the slower dynamics of the hosts.

Mathematically the central manifold approximation in was used ^[9] to improve the quasi-steady-state assumption (QSSA) results, see also ^[10]. Here we use the geometric singular perturbation method where asymptotic expansions give a better approximation when the time-scale difference is moderate, see also ^[10].

The change from the disease-free into the endemic state of the population under variation of the force of infection parameter is studied in the context of bifurcation theory. Instead of the threshold value for the force of infection where the basic reproduction number equals one, $R_0 = 1$ we calculate this value by performing a bifurcation analysis. In this case an eigenvalue of the Jacobian matrix evaluated at an equilibrium of the autonomous system is zero.

Besides time-scale separation, another important issue is the effect of seasonality. In ^[11,12] seasonality in a model similar to the one in ^[9] is studied. In ^[13] it is concluded that assuming parameter stationarity in systems with approximate stationarity in disease incidence is unjustified and may result in missed opportunities to understand the drivers of disease variability. The effect of seasonality is studied where the force of infection changes sinusoidally in time, modeling the annual fluctuations of the presence of the vector. The resulting non-autonomous system is studied in the same way as the autonomous system whereby equilibria are replaced by limit cycles and eigenvalues by Floquet multipliers.

2.   Description of the models

In [9, Sec. 3] an elementary susceptible-infected-removed SIR model for the host where host individuals start off susceptible, at some stage catch the disease, and after a short infectious period become permanently immune, is coupled with also a minimal susceptible-infected SI model for the vector, where mosquito-vectors start off susceptible, at some stage catch the disease and remain infected for the duration of their lifetime. The resulting model is called the SIRUV model.

This simple model is formulated to study the effects of the difference in the time scales of the epidemics in the host population and in the vector population. The total human population is denoted $N = S+I+R$, assumed to be constant. The compartments of the model are formed by susceptible $S$, infected $I$ and recovered $R$, all non-negative. The total vector population is denoted by $M = U + V$ where $U$ and $V$ are the susceptible and infected disease vectors.

The law of mass action is assumed to model random encounters between vector and host. This model can be reduced to an equivalent three-dimensional model:

In order to analyze this model in more detail a simplified model was formulated with a two-dimensional equivalent model [9, Sec. 2]. In this model the SIR model for the host is replaced by a SIS model where host individuals start off susceptible, at some stage catch the disease, and after a short infectious period become susceptible again (no protective immunity), called SISUV model, where $\gamma = 0$: that is, it assumes no recovery compartment and allows the infected hosts after successfully defeating the disease to re-enter immediately into the compartment of the susceptible:

Note that assuming $\gamma = 0$ turns the model (2.1) into (2.2), which is an SI model with no recovery, although formally it appears as if it were indeed an SIS model with the birth/death rate $\mu$ playing the role of the recovery rate.

Alternatively, using the constant host and vector populations, we obtain:

3.   Singular perturbation problem

In this work the geometric singular perturbation technique is used to analyze the SIRUV model (2.1) and the SISUV model (2.3). We will also analyze the quasi-steady state solutions. The reader is referred to ^[14,15] for introductions into perturbation analysis and ^[16] for an illustrative application to the Rosenzweig-MacArthur predator-prey model ^[17,18].

We start with a short overview of singular perturbation technique. Singular perturbation theory deals with systems whose solutions evolve on different time scales whose ratio is characterized by a small parameter $\varepsilon > 0$. When $\varepsilon\ll 1$ the system is a fast-slow system. For instance, in our epidemic models (2.1-2.3) the mosquito-vector dynamics occurs at a much faster time scale than the human-host dynamics. To set (2.3) as an illustration to this approach, the original equations can be rewritten in the form

With $\varepsilon = 0$ we have the fast system also called the layer system:

The infected host population remains constant, so that the trajectories are the vertical lines in Figure 1. With a change of time-scale $\tau = \varepsilon t$ for system (2.3), where the resulting system with $\varepsilon \ll 1$ is named the slow system:

After substitution of $\varepsilon = 0$ we get:

This differential algebraic equation (DAE) is called the reduced system. The trajectory is a part of the $f$-nullcline of the original SISUV-system (2.3) shown in Figure 1. However, in the fast-slow context this is a set of equilibria of the fast layer-system (3.2), called critical manifold, in which $V$ acts as a parameter.

These heuristic results suggest the following approach for dealing with the two different time scales. The first step consists in setting $\varepsilon = 0$, which gives the set of fast equilibria of the fast system (3.2) yielding the algebraic equation (3.3a). This is the critical manifold, namely the set of equilibria on the hyperbola $f(V, I, 0) = 0$.

With a good hypothesis (see below for the details) the hyperbola $f(V, I, 0) = 0$ is equivalent to $I = q(V)$ and we can substitute $I$ by $q(V)$ in equation (3.3b). The result is the slow or reduced-system:

Using (3.3b) and the chain rule

from (3.4) we get formally:

For the case $0 < \varepsilon \ll 1$, we follow the geometric singular perturbation technique. Let us consider system (2.3) again. For $\varepsilon = 0$ the $f$-nullcline, see Figure 1, the set

consists of the critical manifold

${\cal M}$ forms a set of equilibria of the fast system (3.2).

We assume that the critical manifold ${\cal M}$ is normally hyperbolic. Then Fenichel's theorem ^[15,19,20] states that there exists $\varepsilon_0$ such that for $0 < \varepsilon < \varepsilon_0$, there are locally invariant manifolds ${\cal M}_\varepsilon$. Using its invariance, the perturbed manifold ${\cal M}_\varepsilon$ can be approximated by asymptotic expansions in $\varepsilon$. With good hypothesis it can be described as a graph

This manifold is invariant when the following equality holds

which yields with Eq (3.1) and $I = q(V, \varepsilon)$ the invariance equation for the SISUV model:

Later we will furthermore assume that also the inverse of function $q(V) = I$ exists, denoted by $p(I) = q^{-1}(I) = V(I)$. That is, besides the implicit function theorem also the inverse function theorem holds in region of the $(V, I)$ plane we are interested in.

4.   Analysis of the SISUV model

The two-dimensional model is given in (3.1). There is a trivial, disease-free equilibrium $I^0 = 0, V^0 = 0$ and $S^0 = N$ and an interior, endemic equilibrium given by:

The biologically relevant domain for (2.3) is $\Omega = \{(I, V)|0\le I\le N, 0\le V\le M\}$. Following [21, Lemma 2.1], one can demonstrate that $\Omega$ is positively invariant under (2.3). The details are not shown here.

Since the system is two-dimensional we calculate the determinant and trace of the Jacobian evaluated at the equilibria, ^[22]. The trace and the determinant of the $2\times 2$ Jacobian matrix $\mathbf{J}$ read

At the trivial equilibrium we have the following stability properties. Evaluated at this equilibrium the trace and the determinant read:

Note that the $N$ and $M$ do not appear in these expressions, but only the force of infections $\beta$, $\vartheta$ and the birth and death rates $\mu, \nu$ of host and vector respectively are present. The trace $\rm{Tr}\; \mathbf{J} < 0$ is always negative. For the reference values given in Table 1 we have $\beta = 2\mu$ and $\vartheta = 2\nu$. Then the determinant is negative:

Hence, the Routh-Hurwitz criterion shows the instability of the trivial equilibrium.

At the interior point given by Eq (4.1) we have

The determinant has just the opposite value of the trivial equilibrium. It changes sign when $\mu\nu = \beta\vartheta$. At this point both equilibria coincide and the trace $\rm{Tr}\; \mathbf{J} < 0$ is again negative. This means that we have $\det{\mathbf{J}} = 0$ at the transcritical bifurcation TC given by

For values of the parameter $\beta$ above this point we have $\det{\mathbf{J}} > 0$ and the interior equilibrium Eq (4.1) is asymptotically stable.

In the epidemiological literature the basic reproduction number $R_0$ represents the number of secondary cases one case generates on average over the course of its infectious period in an otherwise uninfected population ^[23,24]. $R_0$ equals one at the TC-point.

4.1.   Quasi-steady state assumption

The QSSA equilibrium is calculated, assuming $\varepsilon = 0$, from Eq (2.3b)

This means that the fast variable $V$ changes instantaneously to this value and we substitute this equilibrium value $V^\ast_{\rm QSSA}$ in Eq (2.3a) yielding:

After $I(t)$ is calculated by solving the single ODE Eq (2.3) and thereafter Eq (4.7) gives $V(t)$.

Hence, the time-scale argument gives a further reduction to an equivalent one-dimensional model. When we use the zero-order approximation we got the classical QSSA assumption. Then we finally get using Eq (3.3) the single ODE for the SISUV-system with large time-scale differences between the host and the vector for the infected host population. Observe that in this one-dimensional version the parameters $\nu$, $\vartheta$ from the vector modeling are present in the hyperbolic expression but not the size of the vector population $M$. This hyperbolic expression models a saturated incidence rate. It resembles the Holling type Ⅱ functional response well known in predator-prey models. In ecology the quantity $N\nu/\vartheta$ is called the half-saturation constant. A mechanistic derivation is obtained when besides searching time also handling time of the predator to digest a prey individual is introduced. In biochemistry, Michaelis-Menten kinetics is one of the best-known models of enzyme kinetics. Then the reaction rate is an hyperbolic expression in the substrate concentration. The Michaelis constant is again the substrate concentration at which the reaction rate is half of the maximum rate at abundant substrate. In the epidemiological literature it has also been considered to describe the infection mechanism with a saturated incidence rate ^{[25,26,27,28]}.

Motivated by Holling's derivation it is tempting to follow a similar approach in to obtain the hyperbolic expression in epidemics. In [24, p.155 and p.277–280] it is argued why this is impossible. Namely, it is due to the fact that in the ecological predator-prey setting the two populations are treated independently, whereas in the epidemiology there is one population in which the individuals can reside in two states, namely susceptible and infective.

Interestingly for a vector-borne disease we have here a simple mechanistic formulation where, as in the derivation of the Holling type Ⅱ, a time-scale argument plays a crucial role.

4.2.   Asymptotic expansion for the SISUV model

The following asymptotic expansion in $0 < \varepsilon \ll 1$ is introduced:

hence

Substitution into the invariance equation (3.8) gives

Gathering zero orders of $\varepsilon$ and assuming $V > 0$ results in the following formula, accurate of order $\mathcal{O}(1)$:

and this gives

which are Eq (3.4) and Eq (3.5), respectively.

Gathering the first orders of $\varepsilon$ and assuming $V > 0$ results in the following formula, accurate of order $\mathcal{O}(1)$:

or using Eq (4.12)

or, further

or, finally:

In order to get the slow dynamics along this curve we can solve the ODE for the reduced one-dimensional system

where $p(I) = q^{-1}(I) = V$.

Alternatively the same asymptotic expansion procedure as for $I = q(V, \varepsilon)$ yielding Eq (4.9) can be used for $p(I, \varepsilon)$. The results are not shown here.

5.   Analysis of the SIRUV model

We recall the three-dimensional SIRUV model (2.1)

There is a trivial equilibrium $S^0 = N, I^0 = 0, V^0 = 0$ and an interior equilibrium (whenever $\nu(\mu+\gamma) < \vartheta\beta$) given by:

and furthermore $R^\ast = N-(S^\ast+I^\ast)$ and $U^\ast = M-V^\ast$.

The biologically relevant domain for (2.1) is $\Omega = \{(S, I, V)|0\le S+ I\le N, 0\le S\le N, 0\le I\le N, 0\le V\le M\}$. Following [21,Lemma 2.1], one can demonstrate that $\Omega$ is positively invariant under (2.1). The details are not shown here.

Theorem 1. When the interior equilibrium (5.2) exists, it is locally asymptotically stable and a spiral as long as $\mu$ is sufficiently small.

Proof. We provide the proof for the reference values $\beta = 2 \gamma, \vartheta = 2 \nu$ in Table 1. The Jacobian computed at the interior equilibrium gives a characteristic polynomial

Since $\mu$ is the smallest parameter in the coefficients of $J$, all coefficients are {positive} and Descartes' rule of signs shows that $J$ has no positive real roots. Similarly, because $\gamma$ is the dominant parameter, the coefficient of $\lambda^2$ is the largest with respect to those of $\lambda^3$ and $\lambda$, hence, the local extrema of $J$ are both negative. A result on the geometric interpretation of the complex conjugate roots of a cubic ^[29], shows that these have negative real parts. Therefore, the interior equilibrium $(S^\ast, I^\ast, V^\ast)$ is a spiral.

For the reference parameter values, the slow variables $S, I$ enslave the fast variables $V$ (Figure 8). There is rapid approach to the locus $(V^\ast(S), I^\ast(S))$ shown by the dotted line and fast switching between values near the origin and the interior equilibria.

5.1.   Singular perturbation of SIRUV model

Applying the singular perturbation technique for the SIRUV model (5.1) is now more complicated than for the SISUV model (3.1). The slow manifold is in this case two-dimensional. Using the time-scale argument similar as for Eq (3.3) we get for $\varepsilon = 0$ the two-dimensional $f$-nullspace, the set

which is the critical manifold $\cal{M}$ of the system (2.1) with

In case of the SIRUV model, using the QSSA approach with $\varepsilon = 0$, the two-dimensional slow manifold is

Note that apparently $I$ grows unbounded when $V\rightarrow M^-$. But this is in fact prevented by the constraint $I\le N$, see (5.3). The system (5.4) with $V = \frac{MI}{\nu N/\vartheta+I}$ from (5.5) is the reduced-system. This hyperbolic expression for $V(S, I)$ with saturation with respect to $V$ is the same as that for $V(I)$ in the SISUV model (4.12).

5.2.   Host-only model

It is interesting to compare model (2.1) with the simplest host-only model where the vector dynamics is taken into account implicitly, and where contact rates between the host individuals are described by the law of mass action. To this end we use the linearized version of the inverse of the expression $I = \frac{\nu NV}{\vartheta(M-V)}$, namely $V = \frac{MI}{\nu N/\vartheta+I}$ at its origin, see also the QSSA assumption (4.7). This gives for the hidden infected vector population size $V$:

and the host-only model:

This version can be considered as a simplified version of the model formulated and analyzed in ^[3] when the size of the infected vector population is proportional to the infected host population.

5.3.   Singular perturbation analysis SIRUV model

The algebraic expression (5.5) for the two-dimensional slow manifold ${\cal{M}}$ for the infected host population size $I(V)$ where $\varepsilon = 0$ is the same as given in Eq (3.4) and there are no restrictions for the susceptible population size. Hence, unlike in the SISUV case we can at most express the fast variable $V$ in terms of the slow variables. The statement of Fenichel's theorem ^[15,19,20] under the condition that the critical manifold ${\cal M}$ is normally hyperbolic is as follows. There exists $\varepsilon_0$ such that for $0 < \varepsilon < \varepsilon_0$, there are locally invariant manifolds ${\cal M}_\varepsilon$. Using its invariance, the perturbed manifold ${\cal M}_\varepsilon$ can be approximated by asymptotic expansions in $\varepsilon$.

The invariance equation is similar to (3.8) as follows

where $\frac{d S } {d\tau}$ and $\frac{d I }{d\tau}$ are given in (5.4). With $V(S, I) = p(S, I, \varepsilon)$ we get

5.4.   Asymptotic expansion for the SIRUV model

The following asymptotic expansion in $0 < \varepsilon \ll 1$ is now introduced:

hence

Substituting into the invariance equation (5.8) leads to

Gathering the zero orders of $\varepsilon$ and assuming $V > 0$ gives the following result with $\mathcal{O}(1)$ accuracy:

which in turn yields

and this is the same expression we found for the SISUV model in (4.12). This is also the same as the QSSA result in (5.5).

Gathering the first orders of $\varepsilon$ results in the formula with accuracy $\mathcal{O}(1)$:

or, using (5.13),

Taking only up to the first order terms into account, the reduced model is described by the two-dimensional system (5.4) with $V(S, I)$ given by the following asymptotic expansion up to the first order:

with $p_0(S, I)$ obtained by (5.13) and $p_1(S, I)$ by (5.14). The three-dimensional plot of $p_0, \varepsilon p_1$ for the reference parameter values with $\varepsilon = 1/365$ is shown in Figure 3. The zero-order expression, the tilted surface in the $SIV$-space in the left panel, gives the reduced-system (5.4) where $V(S, I)$ is described by the zero-order asymptotic expansion (5.15). The size of the first-order term in the right panel shows that the the contribution of the term is marginal. Numerical experiments show that the use of the first-order term can be counterproductive already for $\varepsilon > 1/365$ and that spurious equilibria can occur when the trajectory starts not close the the equilibrium.

5.3.   Heuristic analysis of the SIRUV model

Here we present a heuristic analysis of the trajectory of system (2.1) in the two-dimensional slow manifold. The $S, I$ variables for the host population are slow while the variable of the vector population $V$ is fast. However, because the variable $I$ does not appear in the differential equation for $S$, it would be convenient to analyze the behavior of the fast $V$ flow in the $VI$-system, assuming a fixed value of $S$. For a value of the slow variable $S = \bar S$, the $VI$-system reads as

We consider the dynamics of $V, I$ driven by system (5.16) with $\bar S$ acting as a parameter. Depending on the value of $\bar S$, this system has either one or two equilibria: a trivial one $\{(0, 0)\}$ or a pair $\{(0, 0);(V^\ast(\bar S), I^\ast(\bar S))\}$ with the interior equilibrium dependent on $\bar S$ given by

Let $S_c = \frac{\nu(\mu+\gamma)}{\vartheta\beta}N$. It is important to note that for biological relevance we have the constrains $V(t) > 0$ and $I(t) > 0$ and that for $S < S_c$ both are negative in equilibrium. Furthermore, the transcritical bifurcation is degenerate at $S = S_c$ because $V = 0$ and $I = 0$ at the same time.

Theorem 2. The biologically relevant domain of (5.16) is $\Omega = \{(V, I)|0\le V < M, \; 0 < I < (N \beta)/(\mu+\gamma)\}$. Given initial data in $\Omega$, the trivial equilibrium $(0, 0)$ is the single local and global asymptotically stable equilibrium of (5.16) when $\bar S\le S_c$, and otherwise the interior equilibrium $V^\ast(\bar S), I^\ast(\bar S)$ is locally and globally asymptotically stable.

Proof. Denote by $\mathbf{F}$ the flow of the system $(\frac{d}{dt} V = f(I, V), \frac{d}{dt} I = g_2(I, V))$ defined by (5.16). It is straightforward to verify that the boundary of $\Omega$ is impenetrable from the interior under $\mathbf{F}$ so the biologically relevant domain $\Omega$ is positively invariant under $\mathbf{F}$.

Consider first a fixed $0 < \bar S\le S_c$. To show local asymptotic stability of $(0, 0)$, examine the eigenvalues of its respective Jacobian

Its characteristic polynomial is

which cannot have roots with positive real part.

Fix $\bar S > S_c$. To show local asymptotic stability of the interior equilibrium $(V^\ast(\bar S), I^\ast(\bar S))$, examine the eigenvalues of its respective Jacobian

Its characteristic polynomial is

which cannot have roots with positive real part.

To show the global behavior it is enough to exclude the possibility of a closed orbit in $\Omega$. It is straightforward to verify that

in $\Omega$, so Bendixson's negative criterion [30,p. 329] shows that there is no limit cycle in $\Omega$, and all trajectories starting inside $\Omega$ converge to the asymptotically stable equilibrium.

The phase portrait of the system (5.16) is shown in the Figure 4 for changing values of $\bar S$. Observe that the $I$-nullcline changes with $\bar S$, but the $V$-nullcline is independent of $\bar S$.

When $\bar S$ is increasing, the $I$-nullcline rotates counterclockwise and system (5.16) converges to the interior equilibrium (Figure 4, left panel). When $\bar S$ drops beyond the threshold value $S_c$, the $I$-nullcline rotates clockwise, leading to disappearance of the interior equilibrium, and system (5.16) converges to the origin (Figure 4, right panel).

To summarize, system (5.16) undergoes a bifurcation at $\bar S = S_c$. From the trivial branch $B_1 = \{(\bar S, 0, 0)|0\le \bar S\le N\}$ a second branch of interior equilibria $B_2 = \{(\bar S, V^\ast(\bar S), I^\ast(\bar S))|\bar S > S_c\}$ emanates. Note that $B_1$ consists of 2 parts with different stability properties: when $S < S_c$, $B_1$ is locally asymptotically stable, and when $S > S_c$, $B_1$ is a saddle.

Next we characterize the dynamics of the slow variable $S$ in terms of the fast variable $V$. As the dependence of $S$ is explicit on $V$ in (2.1), it is sufficient to consider the dynamics in the $SV$-plane. Consider the locus of points $\{(S, V):\mu(N-S)-\frac{\beta}M V S = 0\}$ representing the $S$ nullcline in the system (2.1) (the dotted curve in Figure 5). First, for $(S, V)$ below this curve, $dS/dt$ is positive and the slow flow is oriented to the right, and for $(S, V)$ above this curve, $dS/dt$ is negative and the slow flow is oriented to the left.

However, when $S$ crosses the critical value $S_c$ the dynamics in the $VI$ plane will rapidly change because of the reappearance or the disappearance of the interior equilibrium. The $V(t)$ value will make a large jump toward the globally asymptotically stable equilibrium of (5.16) (shown as continuous solid line in Figure 5). This would be either an equilibrium in $\{(S, 0, 0):S < S_c\}$ or in $B_2$ (fast $V$-flow shown as long arrows).

Hence, after an initial fast convergence of $(V, I)$ to the locus of equilibria $(V^\ast, I^\ast)$, the slow variable $S$ will evolve monotonically (slow flow shown as short arrows) until the $VI$-system is destabilized and switches between the trivial and interior equilibrium. This switching between the respective trivial and interior equilibria in the $VI$-plane causes the direction of the slow flow to change sign, leading to $S(t)$ spiraling slowly into the equilibrium $(S^\ast, V^\ast, I^\ast)$ which we showed to be locally asymptotically stable.

6.   The role of seasonality

Seasonality is introduced in the non-seasonal, autonomous SIRUV-system by assuming that the size of the vector population changes in a perfect sinusoidal way, this assumption being motivated by dengue fever epidemiology ^[3]. The main motivation lies in being able to reproduce the observed yearly or semiyearly, see ^[12], cycle in dengue incidences. Note that the cause of this yearly periodic changes is outside our modeling approach but includes environmental conditions such as temperature and rainfall.

The expression for the size of the vector population $M(t)$, which is now time-dependent, reads

In ^[3], where no vector dynamics is modeled, the seasonality was introduced as a sinusoidal fluctuation in the host infection rate (based on mass action contact between susceptible hosts $S$ and infected hosts $I$ instead of susceptible hosts $S$ and infected vectors $V$). The expression (6.1) for $M$, with the expression for the infected host based infection rate given in (5.6), is consistent with the expression for the infection rate used in ^[3].

In ^[12], it is shown that the amplitude of the forcing function does not have any substantial influence on the system outcome if a steady state is attained, while simulations indicate that if an oscillatory dynamics arises, the amplitude of the cycles is directly proportional to the one of the forcing function.

7.   Numerical results

In this section we give numerical results of the analysis of the two models: the SISUV model of Eq (2.3) and SIRUV model given by Eq (2.1). To be able to draw conclusions from the obtained results realistic empirical parameter values for the studied eco-epidemiological problem are needed. The reference parameter values used in our study are given in Table 1, with $M_0 = 10000$, and taken from ^[9], which are realistic for the dengue fever case. There is a yearly periodic solution with fixed period equal to 1 year.

In order to apply a time-scale separation technique we have to decouple the three dimensional SIRUV-system (2.1) or (5.1) with state variables $S, I, V$ into one system for the host with state variables $S, I$ and one for the vector with state variable $V$ where the rate of change for the vector is much faster than the one for the host. The ratio of these two rates is defined as $\varepsilon$.

This is most easily done for the SISUV-system (3.1) when we assume $\beta = 2\mu$ and $\vartheta = 2\nu$ and the ratio of the population rates is $\beta/\vartheta$. For the SIRUV-$S, I, V$ system (5.1) the situation is more complex. Here we need to decouple the host $S, I$-system (5.1a) and (5.1b) from the vector $V$-system (5.1c) because $V$ occurs on the right-hand side of system (5.1a) and (5.1b). To that end we can use the QSSA expression (4.7) (effectively the $\varepsilon = 0$ case). The dominant (real part of the) eigenvalues of the Jacobian matrix evaluated at the equilibrium given by (5.2) are the two dominant system rates, given by the return rates to the equilibrium.

With the reference parameter values from Table 1 the calculated time-scale ratio for the host/vector system becomes $\varepsilon = -0.04228/-36.51614 = 0.00116$. In ^[9] and in this paper the ratio of the life expectancies $\mu/\nu$ of the host and vector is used for the SISUV model giving $\varepsilon = 1/365 = 0.00274$. In the case of the SIRUV model we get even $\mu/\nu = 0.000421$. Note that the three values differ by a factor of approximately three and are therefore of the same order. For consistency we used the reference ratio $\varepsilon = 1/365 = 0.00274$ for both models as a natural reference ratio for the host/vector return to equilibrium time-scales.

The small values justify that with the parameters values for Table 1 we can apply a time-scale argument to reduce the full three dimensional $S, I, V$-system (5.1) to the two dimensional $S, I$-system (5.1a) and (5.1b).

7.1.   SISUV model

The trajectory with the reference parameter setting and initial conditions $I = 100, V = 3500$ in shown in Figure 1. From the initial values there is a fast phase where the variable $I$ is almost constant. This fast transient episode ends close to the nullcline $f(V, I, \varepsilon)$ and the trajectory follows this slow manifold and rests at the equilibrium point. Considered as an initial value problem for points far from the slow manifold the first episode is simple, there is an instantaneous reset at the quasi-steady state equilibrium, but this holds only for small $\varepsilon$.

For moderate $\varepsilon$ values we can still do the reduction, but only when the initial condition is close but not at the slow manifold. To that end we use the parameter expansion in $\varepsilon$ explained in the previous section. For larger $\varepsilon$ no reduction is obtained, one has to solve the two-dimensional original system either analytically or numerically.

We study the case with parameter values in Table 1 except for the ratio between the rates of the host and vector which is now $\varepsilon = 60/365$ instead of default value $\varepsilon = 1/365$. The results are shown in Figure 6. The initial point is the same as in Figure 1. The difference in time scale is reduced and $\varepsilon$ is moderately small but not very small as can be seen from the initial part of the trajectory, which differs from a vertical orbit a lot. Nevertheless, when the trajectory approaches the stable equilibrium it is for a large region along the calculated invariant manifold ${\cal M}_\varepsilon$ with first order approximation in the asymptotic expansion in $\varepsilon$ derived here, which is itself an approximation of the stable manifold belonging to the dominant negative eigenvalue of the Jacobian of the stable equilibrium. This approximation is better than that obtained for the zero-order approximation with $\varepsilon = 0$.

7.2.   SIRUV model

7.2.1.   Phase-space analysis

Figure 7 gives a trajectory of SIRUV-system (2.1) in the three-dimensional phase-space. With parameter values in Table 1 there is convergence to the spiral stable interior equilibrium (5.2). The projection of the trajectory on the $I, V$-plane is shown in Figure 2. This shows that the asymptotic dynamics occurs close the the plane $V = q(I)$ given in (4.7) with the saturated incidence rate. Figure 8 gives the results where this relationship is replaced by the linear relationship (5.6) without incidence saturation.

7.2.2.   Bifurcation analysis

The numerical bifurcation analysis was performed using AUTO ^[31] with the reference parameters values from ^[9] given in Table 1 for the SIRUV model. The one-parameter diagram for parameter $\beta$ in Figure 9 shows equilibrium values of the susceptible hosts $S$, infected hosts $I$ and the infected vectors $V$. The solid lines represent the non-seasonal system and the dashed lines denote the maximum and minimum values in the annually fluctuating system.

In both the non-seasonal and seasonal case there is a transcritical bifurcation TC that separates the parameter space into two regimes where the system is disease-free and endemic. In the non-seasonal case the TC occurs where one eigenvalue of the Jacobian matrix evaluated at the disease-free equilibrium is zero and in the seasonal case where one Floquet multiplier of the limit cycle equals one.

The numerical results suggest that both TC's occur at the same parameter threshold value. For the non-seasonal case, below the TC the system is disease-free, that is $S = N$ and $U = M$, and both infected $I = 0$ and recovered $R = 0$ host populations and infected vector $V = 0$ are extinct. For the seasonal case, below the TC the susceptible vector population $U(t)$ oscillates with the same phase and amplitude of the forcing given by Eq (6.1) $U(t) = M(t)$, while the infected population is extinct $V(t) = 0$. For the host population the susceptible population is constant at $S(t) = N = 10000$ while the infected $I(t) = 0$ and the recovered $R(t) = N-S(t)+I(t) = 0$ are extinct.

Above the TC the disease is endemic and all populations oscillate between minimum and maximum values shown in Figure 9 with period of 1 year, equal to that of the forcing function.

In Figure 10 the resulting limit cycles by seasonal forcing are shown for various parameter $\rho$-values where $\beta = 104$.

The two-parameter bifurcation diagram for the two infection rates $\beta$ of the host and $\vartheta$ for the vector is shown in Figure 11. The transcritical bifurcation TC-curve for the non-seasonal system is described by Eq (4.6). The TC-curve for the seasonal case is almost the same; the curves in both cases are indistinguishable in the figure.

8.   Discussion and perspectives

In many dengue epidemic models neither the epidemics nor the ecology of the vector population is modeled explicitly. In the case of dengue fever one vector is the mosquito Aedes aegypti which has a very complex life-cycle with four stages: egg (embryo), larva, pupa (embryo) (in the water phase) and imago (adult). When effects of vector-mosquito control are investigated, an accurate vector population model describing the multiple dengue-strain virus outbreak and spread coupled with a human-host population model is required. This paper is a first step toward an eco-epidemiological approach where we use simplified models for both the host and the vector dynamics.

In ^[9] it was investigated how far the fast time scale of the mosquito epidemiology can be slaved by the slower human epidemiology, so that for the understanding of human disease data, mainly the dynamics of the human time scale is essential and only slightly perturbed by the mosquito dynamics. It was concluded that in a very simple mono-strain SIS host model including vector dynamics, this vector dynamics is slaved by the host dynamics, owing to its slower evolution in time.

In this paper we used the parameter values from ^[9] which are realistic for the dengue fever case. This means that the results from that paper can be directly compared to those obtained here, which in some cases are supplementary. In ^[9] a center-manifold approach is used via eigenvectors of zero eigenvalues at the equilibrium point while here an invariant manifold together with an asymptotic expansion in the perturbation parameter is used. This facilitates the approximation when the differences between time scales of host and vector are moderate and not very large. The results obtained show that in those cases modeling of the vector becomes important. Furthermore, as demonstrated by experimental host-pathogen models ^[32,33,34], pathogenic virus aligns its behavior differently in the target cells of the vector and the host in order to maximize its chance of further transmission. Additional effects such as pest control or personal protection also influence the vectors, factors that speak in favor for an explicit inclusion of the vector dynamics into the epidemic model.

The zero-order assumption, equivalent to the quasi-steady-state, gives a mechanistic derivation of the hyperbolic relationship where the number of contacts saturate with increasing population densities. This resembles the derivation of the Holling type Ⅱ functional response in the case of a predator-prey system.

Another result is that the transcritical bifurcation of the non-seasonal and seasonal forced system by annual abundance cyclic change of the number of mosquitoes that separates the parameter space into two regimes where the system is disease-free and endemic occur at the same parameter value.

In the near future we will use this approach and the obtained results with the analysis of a much more complex eco-epidemiological model for dengue fever by taking explicitly the vector dynamics into account. Examples are the evaluation of vaccines, vector control by application of insecticides, releasing of male mosquitoes infected by a special strain of Wolbachia ^[35] and potential release of genetically modified mosquitoes ^[36].

Acknowledgement

This publication is based upon work from COST Action CA16227 Investigation & Mathematical Analysis of Avant-garde Disease Control via Mosquito Nano-Tech-Repellents, supported by COST (European Cooperation in Science and Technology). Weblink www.cost.eu. Peter Rashkov acknowledges financial support from the Bulgarian National Fund for Scientific Research (FNI) under contract DKOST01/29 and would like to thank the Mathematical Biosciences Institute (funded from the National Science Foundation Division of Mathematical Sciences and supported by The Ohio State University) for the helpful discussions during the Emphasis Semester on Infectious Diseases: Data, Modeling, Decisions-Spring 2018. Maira Aguiar has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 792494. Nico Stollenwerk is supported by National Funding from FCT-Fundação para a Ciência e a Tecnologia, under the project: UID/MAT/04561/2019, and supported as an Investigador FCT.

Conflict of interest

The authors declare there is no conflict of interest.