![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
Abstract
We introduce a compartmental host-vector model for dengue with two viral strains, temporary cross-immunity for the hosts, and possible secondary infections. We study the conditions on existence of endemic equilibria where one strain displaces the other or the two virus strains co-exist. Since the host and vector epidemiology follow different time scales, the model is described as a slow-fast system. We use the geometric singular perturbation technique to reduce the model dimension. We compare the behaviour of the full model with that of the model with a quasi-steady approximation for the vector dynamics. We also perform numerical bifurcation analysis with parameter values from the literature and compare the bifurcation structure to that of previous two-strain host-only models.
1. Introduction
Dengue fever is a tropical disease caused by an arbovirus and spread by female mosquitos of the genus Aedes. Its epidemiology is characterized by the possibility of multiple reinfections of an individual due to the presence of four viral strains which co-circulate in the tropical regions of the Earth [Citation34]. Individuals who have experienced a primary infection with one dengue strain develop life-long immunity against this strain, but after a period of temporary cross-immunity against the other strains become again susceptible to a different dengue strain. Secondary infections may result in health complications such as dengue haemorrhagic fever and pose a health concern in dengue-endemic countries. Mathematical study of dengue epidemiology is motivated by its increasing global incidence in the past decades and recent explosive outbreaks [Citation35].
Epidemiological models often have to find a balance between accurate description of the disease dynamics, the different scales of modelling (within-host/micro scale or population/macro scale), and the associated levels of complexity. By using different dimension reduction arguments to construct mathematically tractable models, it is possible to focus on different aspects of the disease dynamics. While some host-only models include vector dynamics implicitly in the model parameters [Citation1,Citation2,Citation4,Citation18,Citation22,Citation28], etc., studies of effects such as pest control or personal protection via repellents, which influence the mosquito dynamics, speak in favour of an explicit inclusion into the model. In this way, control measures could be tailored to the vector's ecology to achieve an optimum effect as they alter diverse vector characteristics, parametrized in the model as birth and death rates, biting rate, infection rate, etc.
Models of multi-strain dengue epidemiology are based on a susceptible-infected-recovered (SIR) model for the host with the following modification: after recovery from a primary infection with one of the strains, the human hosts gain temporary cross-immunity against both strains, but afterwards become susceptible to possible secondary infection with the other strain. Only after recovery from the secondary infection does a host obtain immunity. Dengue models with two strains and their variations are proposed and studied in [Citation1,Citation2,Citation4,Citation22,Citation23]. These compartmental models are host-only of SIR type. The infected classes distinguish between the viral strains in both primary and secondary infections but do not incorporate explicitly the vector dynamics of mosquitos.
Female mosquitoes, however, play a role in the transmission of the disease and are principal target for most control strategies. Dengue models which include mosquito dynamics explicitly [Citation5,Citation14,Citation15,Citation31,Citation36] add compartments for susceptible and dengue-infected mosquitos. The vector population is governed by a susceptible-infected model because the mosquitos do not recover from the virus. Control strategies for models of single-strain dengue with seasonality have been studied in [Citation5,Citation31] and for a two-strain model without seasonality in [Citation36].
In this study, we consider a variation and extension of these models by introducing susceptible and infected vector classes into host-only models which include susceptible hosts with a background of dengue infection as in [Citation1,Citation2,Citation22,Citation23]. Further, we introduce into the model a parameter that quantifies the differences in likelihood for host-to-vector transmission between infected individuals in a primary and a secondary dengue infection, unlike [Citation36]. Such a parameter serves as a measure whether hosts with a secondary infection have on average a lower or a higher contribution to the overall force of infection upon mosquitos than hosts with a primary infection. This requires besides the new equations also a reinterpretation of the model parameters from [Citation1,Citation4,Citation22,Citation23] which must reflect the per-vector rather than per-host contribution to force of infection. In this way, we incorporate micro-scale factors (immune response) besides macro-scale factors (climate effect on the vector). This allows us to examine parameter regimes corresponding to different scenarios (Section 5). Our model is able to capture phenomena such as seasonal variation and irregular patterns in the incidence of the dengue observed in regions such as Thailand, as well as the co-circulation of dengue strains.
We analyse properties of the full model and study its dynamics in relation to the simplified host-only models, as well as to reduced versions based on quasi-steady state assumptions and geometric singular perturbation analysis as done previously [Citation27]. We also compare the qualitative behaviour of our model that explicitly includes a vector population to those host-only models by performing numerical bifurcation analysis with parameter values from the literature [Citation1,Citation2,Citation22,Citation23]. We also study the bifurcation structure of the autonomous model to a non-autonomous model where seasonal variation in the vector population exists. Our analysis leads us to the conclusion that in order to remain in agreement with seasonally varying data on dengue fever incidence, to reflect disease characteristics of primary vs. secondary infections, and to understand the epidemic complexity, it is essential to include vector ecology in the model. In this way, long-term control strategies for a multi-strain vector-borne disease such as dengue fever could become realistically tractable based on such models.
2. Mathematical model
The host population is subdivided into classes with the following transitions (Figure , left): susceptible individuals without a previous dengue infection , primary infected with ith strain
, and recovered
from a primary infection respectively with strain
. After a period of a temporary cross-immunity, recovered individuals in class
transition to the class
of susceptible individuals with a history of a primary dengue infection with ith strain i = 1, 2. These individuals can become infected again with the other virus strain: hosts from
transition to class
and from
to class
. Finally individuals from the secondary infected classes
recover and obtain life-long immunity against both strains (class R).
Figure 1. Flow chart of the transitions between the host and vector compartments. The birth and death of hosts and vectors are not explicitly shown for the sake of clarity.
![Figure 1. Flow chart of the transitions between the host and vector compartments. The birth and death of hosts and vectors are not explicitly shown for the sake of clarity.](/cms/asset/e3505fa1-b75c-4f93-8a51-1c894d11c907/tjbd_a_1864038_f0001_ob.jpg)
To model the vector dynamics, we follow [Citation15] in assuming that once a mosquito is infected, it never recovers and cannot be reinfected with a different strain because of its short lifespan. We label the susceptible vector population by U and the vector populations carrying i-th dengue strain by (i = 1, 2). Reinfections with the dengue virus of a different strain, therefore, may take place only in the host. We also assume that upon biting an infected individual, a mosquito is infected with the strain that infected individual currently carries (Figure , right).
Hence, transmission of the first dengue strain occurs between vector class to host classes
and between host class
to vector class U. Similarly, that transmission of the second strain occurs between vectors in class
to host classes
and between hosts in classes
to vectors in class U.
The cumulative disease transmission parameter (j is the number of consecutive infection a human experiences) is defined as the product of the biting rate r and the per-bite infection probability
from mosquito vector to human host:
. We assume that
, in other words, the infection probability in a secondary dengue infection is equal to or greater than that in a primary infection due to the phenomenon of antibody-dependent enhancement (ADE) [Citation20] leading to
. Thus, in contrast to [Citation23] the index j of B in our model does not refer to the infectivity of a specific dengue viral strain. As we do not model such strain-specific infectivity, we assume that both strains have equal infectivity upon a susceptible host regardless of his infection background.
The cumulative disease transmission parameter ϑ is the product of the biting rate r with the per-bite infection probability from human host to mosquito vector:
.
The shape of the force of infection term for mosquitos upon humans and humans upon mosquitos follows [Citation11]. As in [Citation1] individuals recover from any dengue infection after a period of length years and their temporary cross-immunity lasts on average
years. The hosts' population demography (human birth and death rates) is parametrized by μ. For the mosquito vector, the demography (mosquito birth and death rates) is parametrized by ν. We assume
to account for the significantly longer life span of the human hosts.
It is important for such a model to capture some characteristics of the host-to-vector transmission as well. The factors behind a successful transmission of dengue from an infected human to a biting mosquito are poorly understood [Citation26]. No studies have investigated human host-seeking ability in dengue-infected mosquitos [Citation6]. However, factors such as personal protection measures like use of repellents, protective clothing by asymptomatic patients, or effective isolation of infected patients from mosquitos could influence the dynamics of the disease transmission.
Hence, we exploit this opportunity to study mathematically these characteristics. Unlike Zheng and Nie [Citation36] who assume that individuals with primary and secondary infection have equal contribution to the force of infection upon susceptible mosquitos, we introduce a parameter ϕ that quantifies the difference in transmission likelihood from hosts with primary or secondary infections to vectors.
When ϕ is larger/smaller than 1, individuals in their secondary infection are more/less likely to transmit the disease to mosquitos than individuals in their first infection. This could be due to factors such as changes in biting rates or virus transmission probabilities. This differs from the interpretation of the parameter ϕ in models such as [Citation1,Citation4,Citation22] where ϕ quantifies the effect of ADE, but rather relates to its interpretation as disease transmissibility during a secondary infection as defined by [Citation28]. We note the models in [Citation1,Citation4,Citation22,Citation28] do not include explicitly any vector population.
A value corresponds to a scenario where hosts with a secondary infection have on average a lower contribution to the overall infectivity of mosquitos upon humans than hosts with a primary infection. A reason behind this could be that secondary infection is clinically more apparent (due to disease severity) and thus more people seek medical care and need hospitalization [Citation7]. Hence, such individuals have on average a lower chance of exposure to mosquitos on search for a blood meal, and the cumulative rate ϑ is then lowered by a factor of ϕ.
A value corresponds to a scenario where hosts with a secondary infection have on average a higher contribution to the overall infectivity of mosquitos upon humans than hosts with a primary infection. In particular, the transmission likelihood from host to vector could be greater if during the course of a secondary infection, the infected host is more attractive to mosquitos on search for a blood meal, because, e.g. of physiological symptoms such as high tympanic temperature or higher viral titres in the host's blood [Citation26]. In particular, evidence suggests a higher mean maximum body temperature in patients is usually associated with a secondary infection [Citation7].
2.1. Model equations
Initially we assume both the host population size N and the vector population size M are constant in time. Hence, the dimension of system can be reduced by eliminating two of the variables by setting and
. We shall work with the following autonomous system of ordinary differential equations:
(1a)
(1a)
(1b)
(1b)
(1c)
(1c)
(1d)
(1d)
(1e)
(1e)
(1f)
(1f)
(1g)
(1g)
(1h)
(1h)
(1i)
(1i)
(1j)
(1j)
(1k)
(1k) where
denotes the time derivative.
Observe that for the host infection rate the symbol B is used instead of β to avoid a misinterpretation. In host-only models [Citation2,Citation4,Citation18,Citation22] infection rate per infected host β was used and in host-vector models [Citation27,Citation29] as well as here infection rate per infected vector.
Implicitly in [Citation2,Citation22] and other earlier models described in [Citation4,Citation18], the vector dynamics was taken into account by assuming that the size of the infected vector population is proportional to that of the infected human-host population. The two-strain host-only model is described in the Appendix (A1).
In [Citation27], the following equation was derived for the relationship of the infection rates in the host-vector model and the host-only model
:
(2)
(2) Using the values from [Citation1,Citation22] this yields
, and makes it possible to compare the results of the host-vector model (1) with those from for the host-only model (A1) in [Citation22].
The biologically relevant domain for (1) is
Following [Citation24, Lemma 2.1], one can demonstrate that Ω is positively invariant under (1). The details are not shown here.
We also consider some ecological facts about the mosquito vector, namely, the case of a seasonally changing mosquito population due to climate factors such as temperature or rainfall. The number of mosquito vectors M is modelled as seasonal forcing and is given explicitly by a cosine function
(3)
(3) In (Equation3
(3)
(3) ),
is the mean population, η describes the magnitude of the seasonal fluctuation and φ is the phase, which we assume to be zero. The frequency
models yearly variations in the mosquito population. In this case, the system of governing ordinary differential Equations (1) with (Equation3
(3)
(3) ) becomes non-autonomous.
2.2. Basic reproduction number
In the epidemiological literature, the basic reproduction number denoted by represents the number of secondary cases; one infected case generates on average over the course of its infectious period, in a fully susceptible population [Citation9,Citation10]. There are several approaches to compute
. We compute the parameter values at which an endemic equilibrium appears (e.g. through a transcritical bifurcation TC, see Table ), that is, the disease can persist in the population and associate it with the value
.
3. Equilibria and local stability analysis
We perform analysis of the equilibria of the model (1). For sake of shortness, whenever necessary, we denote by
Of interest for the dynamics of the disease is first, the disease-free equilibrium, where all individuals and mosquitos are susceptible,
and all other classes are zero, and second, any endemic equilibria, where some classes of infected individuals and mosquitos are non-zero.
For the sake of clarity in the rest of the paper, we shall denote .
3.1. Symmetry
The system (1) possesses -symmetry, namely, if we denote the symmetry transformation matrix by
and by F the vector of the right-hand side of (1), we have
We recall that an equilibrium/limit cycle
is called fixed when
. Then the state variables in their own classes are equal, for instance,
for all t. If the solution is a fixed limit cycle, then the variables
with i, j = 1, 2 are in phase.
Two equilibria/limit cycles where
, are called S-conjugate if their corresponding solutions satisfy
(and because
also
). There is another type of periodic solution with period T that is not fixed but called symmetric when
In this case, it holds that
etc. (see Figure ). For a detailed discussion of the
-symmetry in dengue epidemiology, the reader is referred to [Citation22].
3.2. Disease-free equilibrium
The trivial, disease-free equilibrium of (1) is .
Proposition 3.1
Let
(4)
(4) The disease-free equilibrium equilibrium
is locally asymptotically stable if
, and locally asymptotically unstable if
.
The proof is given in the Appendix.
At , the local stability of the disease-free equilibrium
changes. Because two virus strains are present in the epidemiological setting, several endemic equilibria with one or both strains present are possible. We consider first the boundary equilibria in Ω with a single strain. From an ecological viewpoint, this is the case of competitive exclusion of one of the strains.
3.3. One-strain, competitive exclusion equilibria
Let denote the equilibrium where the first dengue strain is absent, so
with
, and let
denote the equilibrium where the second dengue strain is absent, so
with
. It is straightforward to show using [Citation27] the following
Lemma 3.2
The hyperplanes and
are positively invariant sets for (1). Furthermore, any trajectory starting in
converges to
whenever μ is sufficiently small.
Because the model is symmetric in the two strains, the stability analysis of is analogous, and we shall analyse the linear stability of
.
In , the equilibrium values of the nonzero classes (see Appendix for the derivation) are
(5a)
(5a)
(5b)
(5b)
(5c)
(5c)
(5d)
(5d)
(5e)
(5e) These values are positive (and they belong to the biologically relevant set Ω) if and only if
. Summarizing the above, we have
Proposition 3.3
The equilibria exist in Ω when
.
Next we show that under a mild assumption, both competitive exclusion equilibria , are locally asymptotically unstable for all
. This is due to the presence of temporary cross-immunity for the recovered individuals in
modelled by the parameter α. Remember that the average duration of the temporary immunity for these individuals equals
years, e.g.
corresponds to a duration of 6 months.
Proposition 3.4
The one-strain competitive exclusion equilibria of (1) are locally asymptotically unstable for all
if
Proof.
Consider exclusion equilibrium . For
, we have
. Denote the Jacobian matrix evaluated at the equilibrium point by
by
(see Appendix).
The product of the eigenvalues of is the same as
. Because
is an odd-sized matrix, if
, it must have a positive real eigenvalue. We compute
Positivity of
together with the condition
implies
. Therefore, if the condition of Proposition 3.4 is met, the matrix
has a positive real eigenvalue.
The analysis of the determinant is analogous for the equilibrium .
Thus, both boundary equilibria (in ecological sense, equilibria of competitive exclusion of one dengue virus strain) are S-conjugate and are always saddle points if the susceptibility to a secondary infection is equal to the susceptibility to a primary infection.
Corollary 3.5
Whenever , the condition from Proposition 3.4 holds automatically, so both boundary equilibria are saddle points.
The stable manifold to each is given by the hyperplane
listed in Lemma 3.2.
Note, however, if the condition from Proposition 3.4 is not met, we cannot conclude the opposite, namely, that the equilibrium is locally asymptotically stable, without a numerical validation.
To sum up, unlike [Citation15], the model (1) allows for both boundary equilibria to be locally asymptotically unstable. The model in [Citation4] arrives to the same conclusion as Corollary 3.5, namely, due to the equal epidemiological parameters for the dengue strains neither boundary equilibrium will ever be stable.
3.4. Interior, two-strain equilibrium
Next, we consider equilibria where both virus strains are present in the population. In these endemic equilibria, the classes of host individuals and mosquito vectors carrying either strain are strictly positive.
Proposition 3.6
Let . The interior endemic equilibrium is given by
where x is the root of
(6)
(6) with coefficients given by
This equilibrium is unique if
. There are two interior endemic equilibria if
and
.
Note that due to the size of the Jacobian we are not able to make conclusions about stability of this equilibrium without resorting to numerical approximation.
When the class of secondary infected individuals is not able to transmit the virus to mosquitos (), the interior equilibrium is a curve:
Proposition 3.7
Let . The set of interior equilibria
is a curve consisting of the points
parametrized by
. Further,
is a locally stable.
An illustrative example of the parameter ranges for existence of single or pair of interior equilibria is provided in Figure .
Figure 2. The range corresponding to 0, 1 or 2 coexistence equilibria (the dashed vertical line is the set with
such that
. The shaded area shows the pairs
with a locally asymptotically stable interior equilibrium (eigenvalues computed numerically). The parameter values are given in Table and
.
![Figure 2. The range (B1,B2) corresponding to 0, 1 or 2 coexistence equilibria (the dashed vertical line is the set with B1 such that R0=1. The shaded area shows the pairs (B1,B2) with a locally asymptotically stable interior equilibrium (eigenvalues computed numerically). The parameter values are given in Table 1 and ϕ=0.5.](/cms/asset/649cac7d-919b-4787-9f5d-b390fd4f06a7/tjbd_a_1864038_f0002_ob.jpg)
Table 1. Parameter values, most after [Citation1,Citation22] for the host-only model (A1) and [Citation29] for the host-vector model.
Details on the proofs of Propositions 3.6 and 3.7 are provided in the Appendix.
3.5. Special case ![](//:0)
![](//:0)
Whenever the per-bite infection probability from mosquito vector to human host in the secondary infection equals that in the first infection, , we have
, and we can simplify the conditions on signs from Proposition 3.6. Denote
Whenever
,
, and otherwise,
if
. Since
is a small parameter, observe that
Hence, if additionally the discriminant of (Equation6
(6)
(6) )
is positive, the quadratic has two positive roots.
Hence, the number of positive solutions of (Equation6(6)
(6) ) for
can be summarized as follows.
Proposition 3.8
Let and
.
Case A. Let . Whenever
, there are no positive solutions
to (Equation6
(6)
(6) ) and hence, no interior equilibria
. If
, there is exactly one positive solution
to (Equation6
(6)
(6) ) and one interior equilibrium
.
Case B. Let . Whenever
, there are no positive solutions
to (Equation6
(6)
(6) ) and hence, no interior equilibria
to (1). If
and
, there are exactly two positive solutions
to (Equation6
(6)
(6) ); in turn, there are two interior equilibria
to (1). If
and
, there are no interior equilibria
to (1). If
, there is exactly one positive solution
to (Equation6
(6)
(6) ), and one interior equilibrium
to (1).
Note: One cannot a priori guarantee that the parameter values always fulfil the relations required by Case A. Observe that whenever ,
, and hence, if ϕ is sufficiently large, it is possible that
. However, the relation
holds for the reference parameter values used in recent dengue models [Citation22,Citation29,Citation30].
In the case , we can make the following estimate for the root x and the endemic equilibrium values of
in the case A in Proposition 3.8.
Lemma 3.9
Let and
. In the endemic equilibrium
.
The proof is given in the Appendix.
4. Dimension reduction
Following the approach in [Citation27], we seek to reduce the model dimension and complexity by the means of a quasi-steady state approximation, and the techniques of geometric singular perturbation analysis [Citation21].
Singular perturbation theory studies systems whose solutions evolve on multiple time scales whose ratio is characterized by a small parameter . When the ratio
, the system is a fast-slow system. For instance, in the epidemiological models considered in [Citation27,Citation29] and in model (1), the vector dynamics occurs at a much faster time scale than the host dynamics.
The equations for the host dynamics from (1) can be recast into the following terms
(7a)
(7a)
(7b)
(7b)
(7c)
(7c)
(7d)
(7d)
(7e)
(7e)
(7f)
(7f)
(7g)
(7g)
(7h)
(7h)
(7i)
(7i) With a change of time scale
for system (7), we obtain equations for the host populations on the slow time scale
and for the vector variables (Equation1j
(1j)
(1j) )–(Equation1k
(1k)
(1k) ) on the fast time scale
(8a)
(8a)
(8b)
(8b)
4.1. Quasi-steady state approximation model
For , Equation (8) become a system of algebraic equations. For this so-called reduced system, we solve a matrix equation for the fast variables
Using the inverse matrix
which exists (whenever
) for all
, the fast variables for the quasi-steady state approximation (QSSA) model are found:
(9a)
(9a)
(9b)
(9b)
Observe that are constant on the lines
. The expressions (Equation9
(9a)
(9a) ) are hyperbolic relationships similar to those of the Holling type II functional response in ecology, but with competitive inhibition from the host populations infected with the other virus strain. The equations for the QSSA model, thus, take the form (Equation1a
(1a)
(1a) )–(Equation1e
(1e)
(1e) ) with the values
given by (Equation9
(9a)
(9a) ).
4.2. Asymptotic expansion using the invariance equation
The set of critical points (the equilibria of the fast system (8), with ) is normally hyperbolic and contains a critical manifold
(10)
(10) We refer to the Appendix A.4 the proof of normal hyperbolicity of the set of critical points.
The statement of Fenichel's theorem [Citation16,Citation17,Citation21] under the condition that the critical manifold (Equation10
(10)
(10) ) is compact and normally hyperbolic is as follows: there exists
such that for
, there exist locally invariant manifolds
that are
-close, diffeomorphic to
and locally invariant under the flow (7).
The invariance equation for can be recast in terms of
as
(11)
(11) where
are given in the slow time scale. Using its invariance, the perturbed manifold
can be represented as a graph of a function with
and approximated by asymptotic expansions in ε.
With this notation, we introduce the following asymptotic expansion in :
(12a)
(12a)
(12b)
(12b) Substitution of these expressions into (Equation11
(11)
(11) ) and gathering the same powers of ε allows us to solve for the coefficients
. Solving for the terms
turns, in fact, to be equivalent to solving a recursive matrix equation
(13)
(13) with the elements of vector
determined by
and their partial derivatives. This recursive procedure can be used to find the terms
up to the desired order.
Substituting and simplifying (see details in the Appendix) we obtain the coefficients of ε-terms in (Equation12(12a)
(12a) ). The coefficients
found in (EquationA10
(A10a)
(A10a) ) in are exactly the QSSA values for
which were computed in (Equation9
(9a)
(9a) ).
Observe that the first-order coefficients (EquationA11
(A11a)
(A11a) ) do not contain all slow variables, so it might be necessary to include higher order terms in order to increase accuracy of the power series approximation (Equation12
(12a)
(12a) ). Otherwise, convergence to spurious equilibria or limit cycles, or perhaps divergence could occur. Second order coefficients are algebraically more involved and their derivation is given in the Appendix.
5. Numerical results and discussion
Here we present and discuss results for the model obtained by numerical simulations. We recall that integration of the system in time gives transient dynamics starting from given initial conditions, while sensitivity analysis gives dependence of model outcomes by variation of parameter values. Bifurcation theory gives dependence of qualitative (long-term dynamics) model outcomes by continuation of the parameter value. The numerical bifurcation analysis results are obtained by the software package AUTO [Citation12] for continuation and bifurcation problems in ordinary differential equations. A list of the bifurcations discussed in this paper is given in Table .
Table 2. List of bifurcation types.
Examples of long-term dynamics in the model are convergence to a stable equilibrium or limit cycles but also chaotic attractors. The bifurcation diagrams depict the dependency of the long-term dynamics on parameter values. The regions with different dynamics are separated by bifurcations points (one-parameter) or curves (two-parameter) which are calculated by continuation.
The different modelling approaches of host-vector or host-only for the application to two-strain dengue fever are assessed by comparing the resulting bifurcation structures. We provide one- and two-parameter bifurcation diagrams in the model parameters and compare the diagrams for the host-only model (A1) of [Citation1,Citation4,Citation22] and the host-vector model (1) considered in this study. We recall that ϕ in the host-only model stands for the ADE ratio describing the contribution of the secondary infection to the force of infection, while ϕ in the host-vector model (1) quantifies the differences in likelihood for host-to-vector transmission between infected individuals with a primary and a secondary dengue infection. In the figures, we denote the total number of infected individuals by
for the sake of shortness.
5.1. Results for autonomous systems
Here we discuss the model under the assumption of a constant vector population in time. All results for the autonomous system (1) use reference parameter values given in Table from [Citation1,Citation22,Citation29]. The results are compared with those for earlier models in the literature: [Citation1,Citation4,Citation18,Citation22]. Important is to recall that in [Citation2] a realistic description of recorded disease incidents was produced for empirical data sets of dengue fever in Thailand.
In Figure , the two parameter diagrams are shown for the host-only model (A1) from [Citation22] (panel a) and host-vector model (1) (panel b).
Figure 3. Two-parameter bifurcation diagram for the range
and
. One Hopf bifurcation converges for
to the ϕ parameter value [†] for the model of [Citation4] for (a) the host-only model (A1) from [Citation1] and (b) the host-vector model (1). [*] marks the value
in [Citation18].
![Figure 3. Two-parameter (ϕ,α) bifurcation diagram for the range ϕ∈[0,12] and α∈[0,250]. One Hopf bifurcation converges for α→∞ to the ϕ parameter value [†] for the model of [Citation4] for (a) the host-only model (A1) from [Citation1] and (b) the host-vector model (1). [*] marks the value ϕ=2.5 in [Citation18].](/cms/asset/c1ed4faf-304e-454e-9139-0b84398c3c68/tjbd_a_1864038_f0003_oc.jpg)
The bifurcation pattern for the host-only model in Figure (a) is dominated by two Hopf bifurcations: one supercritical and the other
. Furthermore, there is one Torus bifurcation TR also called a Neimark-Sacker bifurcation. These results have been discussed in [Citation22] and the reader is referred to this paper for a detailed discussion. The Hopf bifurcations occur also for the host-vector model Figure (b) but there is no torus bifurcation TR.
Most important is the parameter region for and
. The bifurcation pattern shown in Figure is much more complex in this region for the host-only model (A1) in Figure (a). These results have been discussed in [Citation2] and the reader is referred to this paper for a detailed discussion.
Figure 4. Two-parameter bifurcation diagram for the range
and
(a) for the host-only model (A1) and (b) for the host-vector model (1).
![Figure 4. Two-parameter (ϕ,α) bifurcation diagram for the range ϕ∈[0,1.3] and α∈[1,7] (a) for the host-only model (A1) and (b) for the host-vector model (1).](/cms/asset/72044b7b-6799-40aa-86b3-a867d8fb11a5/tjbd_a_1864038_f0004_oc.jpg)
The resulting pattern for the host-vector model (1) in Figure (b) is much simpler with only the supercritical Hopf bifurcation and a pitchfork bifurcation P and again not a torus bifurcation. This is important because the presence of a torus bifurcation makes chaotic dynamics via a Ruelle-Takens-Newhouse route [Citation32] possible.
For the host-only model with , the one-parameter bifurcation diagram is shown in Figure (a). There is a rich bifurcation scenario that leads to chaos in the interval
. These results have been discussed in [Citation2] and the reader is referred to this paper for a detailed discussion.
Figure 5. One-parameter diagram for ϕ where and B = 52. Equilibria or maximum and minimum values for limit cycles for total infected I (a) for the host-only model (A1) and (b) for the host-vector model (1). Red indicates stable and blue unstable solutions.
![Figure 5. One-parameter diagram for ϕ where α=2 and B = 52. Equilibria or maximum and minimum values for limit cycles for total infected I (a) for the host-only model (A1) and (b) for the host-vector model (1). Red indicates stable and blue unstable solutions.](/cms/asset/971c1f3c-93c6-4eaf-a8c3-07eb5e292983/tjbd_a_1864038_f0005_oc.jpg)
Figure (b) gives the one-parameter bifurcation diagram for the host-vector model (1) for . For
, the set of interior coexistence equilibria is given by the curve
, parametrized by
. We have shown in the Appendix that
is a locally attracting manifold. Increasing the values
the unique interior equilibrium
becomes unstable at the Hopf bifurcation
. Above this Hopf bifurcation the limit cycle is stable and the maximum and minimum values are shown in Figure (b). Further continuation gives a stable limit cycle which becomes unstable at a pitchfork bifurcation P. At this point, two new stable limit cycles originate. In reverse order, it occurs at a second pitchfork bifurcation P. At
the unstable limit cycle is symmetric (shown in Figure (a)). The pairs of two stable S-conjugate limit cycles are shown in Figure (b).
Figure 6. (a) Infected classes (red),
(blue) unstable symmetric solutions for the limit cycle with parameter values
,
, B = 52 and normalized time for one period. (b) Infected classes
(red),
(blue) and
(green),
(pink) solutions for the two stable S-conjugate limit cycles.
![Figure 6. (a) Infected classes I1(t) (red), I2(t) (blue) unstable symmetric solutions for the limit cycle with parameter values α=2, ϕ=0.8, B = 52 and normalized time for one period. (b) Infected classes I11(t1) (red), I21(t1) (blue) and I12(t2) (green), I22(t2) (pink) solutions for the two stable S-conjugate limit cycles.](/cms/asset/42b66145-6672-4018-985f-b58af85d234b/tjbd_a_1864038_f0006_oc.jpg)
No chaotic attractors exist for the host-vector model in the relevant region in the parameter space studied for . This results from the non-existence of a torus bifurcation such as in the host-only model (A1) with temporary cross-immunity for the second strain infection [Citation2].
Under the assumption of equal host infectivity in a primary and secondary infection, we employ the parameter as a bifurcation parameter. In Figures , we compare the bifurcation structure of the full model (1) with that of the QSSA model for some values
. We plot the combinations of pairs
such that (a)
, (b) a Hopf bifurcation or (c) a pitchfork bifurcation occurs in (1) and its QSSA counterpart in Figures –. These bifurcation diagrams find application, in particular, for evaluation of control strategies which act on the vector-to-host and host-to-vector transmission of the virus, based potentially on application of anti-mosquito repellents [Citation3,Citation36].
Figure 7. Two-parameter bifurcation diagram with parameter
of (a) the system (1), and (b) the QSSA system where the fast variables
are assumed to be in steady state (Equation9
(9a)
(9a) ). TC curve is for values
. Left of TC there is stable disease-free equilibrium where I = 0. Between the TC and H curve, there is a stable endemic equilibrium, and beyond H, there is a stable limit cycle. No pitchfork bifurcation P is observed in this parameter region.
![Figure 7. Two-parameter (B,ϑ) bifurcation diagram with parameter ϕ=0.2 of (a) the system (1), and (b) the QSSA system where the fast variables V1,V2 are assumed to be in steady state (Equation9(9a) V1=ϑM(I1+ϕI21)ϑ(I1+I2+ϕ(I21+I12))+νN,(9a) ). TC curve is for values R0=1. Left of TC there is stable disease-free equilibrium where I = 0. Between the TC and H curve, there is a stable endemic equilibrium, and beyond H, there is a stable limit cycle. No pitchfork bifurcation P is observed in this parameter region.](/cms/asset/6a707a87-9d90-44f7-90d5-538ce43de78a/tjbd_a_1864038_f0007_oc.jpg)
Figure 8. Two-parameter bifurcation diagram with parameter
of (a) the system (1), and (b) the QSSA system where the fast variables
are assumed to be in steady state (Equation9
(9a)
(9a) ). TC curve is for values
. Left of TC there is stable disease-free equilibrium where I = 0. Between the TC and H curve there is a stable endemic equilibrium, and between H and P, there is a stable limit cycle. Inside the region bordered by P, there is bistablity: a pair of S-symmetric limit cycles.
![Figure 8. Two-parameter (B,ϑ) bifurcation diagram with parameter ϕ=0.5 of (a) the system (1), and (b) the QSSA system where the fast variables V1,V2 are assumed to be in steady state (Equation9(9a) V1=ϑM(I1+ϕI21)ϑ(I1+I2+ϕ(I21+I12))+νN,(9a) ). TC curve is for values R0=1. Left of TC there is stable disease-free equilibrium where I = 0. Between the TC and H curve there is a stable endemic equilibrium, and between H and P, there is a stable limit cycle. Inside the region bordered by P, there is bistablity: a pair of S-symmetric limit cycles.](/cms/asset/1f25fd21-b9c1-4970-9e86-f647980cb284/tjbd_a_1864038_f0008_oc.jpg)
Figure 9. Two-parameter bifurcation diagram with parameter
of (a) the system (1), and (b) the QSSA system where the fast variables
are assumed to be in steady state (Equation9
(9a)
(9a) ). TC curve is for values
. Left of TC there is stable disease-free equilibrium where I = 0. Between the TC and H curve, there is a stable endemic equilibrium, and between H and P, there is a stable limit cycle. Inside the region bordered by P, there is bistablity: a pair of S-symmetric limit cycles.
![Figure 9. Two-parameter (B,ϑ) bifurcation diagram with parameter ϕ=0.8 of (a) the system (1), and (b) the QSSA system where the fast variables V1,V2 are assumed to be in steady state (Equation9(9a) V1=ϑM(I1+ϕI21)ϑ(I1+I2+ϕ(I21+I12))+νN,(9a) ). TC curve is for values R0=1. Left of TC there is stable disease-free equilibrium where I = 0. Between the TC and H curve, there is a stable endemic equilibrium, and between H and P, there is a stable limit cycle. Inside the region bordered by P, there is bistablity: a pair of S-symmetric limit cycles.](/cms/asset/aca310dc-cca0-4352-9605-2a2d54441ec4/tjbd_a_1864038_f0009_oc.jpg)
The curve for where an interior equilibrium appears is computed from
using (Equation4
(4)
(4) ). A Hopf bifurcation leads to the appearance of a symmetric limit cycle. The pairs
where the Hopf bifurcation occurs are computed numerically by finding pairs of complex conjugate eigenvalues with zero real part of the Jacobian evaluated at the interior, two-strain equilibrium
.
We observe that the QSSA model undergoes a Hopf bifurcation for values which are closer to the values for the transcritical bifurcation TC than the full model. Also the band of stability in
-space of the resulting limit cycle in the QSSA model is much narrower for a given ϕ and beyond that it becomes unstable due to a pitchfork bifurcation.
Observe that as , the border of the Hopf region in the host-vector model (1) approaches the TC point (Figure (a)). As ϕ increases, the limit cycle in the QSSA model rather quickly becomes unstable (Figure (b)), leading to a pair of S-conjugate limit cycles (compare Figure (b) versus (b)).
Figure 10. One-parameter bifurcation diagram for B for the case of equal infection rates in the host-vector model (1) (column (a)), and the QSSA model (column (b)). Solid lines indicate stable and dashed lines unstable solutions. Above the Hopf bifurcation point the minima and maxima for limit cycles are shown for total infected I for
.
![Figure 10. One-parameter bifurcation diagram for B for the case of equal infection rates B1=B2=B in the host-vector model (1) (column (a)), and the QSSA model (column (b)). Solid lines indicate stable and dashed lines unstable solutions. Above the Hopf bifurcation point the minima and maxima for limit cycles are shown for total infected I for ϕ=0.5.](/cms/asset/f9912cf6-2c59-42d3-8498-ae1686c79c9a/tjbd_a_1864038_f0010_oc.jpg)
Figure 11. One-parameter bifurcation diagram for B for the case of equal infection rates in the host-vector model (1) (column (a)) and the QSSA model (column (b)). Solid lines indicate stable and dashed lines unstable solutions. Above the Hopf bifurcation point the minima and maxima for limit cycles are shown for total infected I for
.
![Figure 11. One-parameter bifurcation diagram for B for the case of equal infection rates B1=B2=B in the host-vector model (1) (column (a)) and the QSSA model (column (b)). Solid lines indicate stable and dashed lines unstable solutions. Above the Hopf bifurcation point the minima and maxima for limit cycles are shown for total infected I for ϕ=0.8.](/cms/asset/43711f0c-5ee3-4365-8517-5accd1dde599/tjbd_a_1864038_f0011_oc.jpg)
The bifurcation structure of the solutions for the host-vector model and its QSSA counterpart is examined in more detail in Figures – by varying the bifurcation parameter B and for two values of ϕ. Beyond the Hopf point H the steady state becomes unstable and a limit cycle emerges. Because of the -symmetry of these two systems, S-conjugate limit cycles are possible (for more details see the two-strain model of [Citation23]). After H with increasing B the limit cycle loses stability via a pitchfork bifurcation P. A pair of S-conjugate limit cycles emerges beyond P. Details on the minima and maxima for the two classes of primary infected
are shown in online Appendix (Figure F2 and F3).
Introduction of control measures such as personal protection by means of repellents or repellent-treated clothing changes the parameter B. The numerical bifurcation analysis shown in Figures – reveal that changing B could lead to a switch of regime from periodic oscillations to stable equilibrium depending on the efficacy of the control. Thus, for the chosen parameter values, peak values of infected individuals may decrease by an order of magnitude
(Figure (a)). Sensitivity analysis is a key tool in the study of model responses to changes in conditions, see [Citation36]. The starting point of a sensitivity analysis is the calculation of the local sensitivities of all paramters involved to quantify their contribution to changes in the model output. This gives here a method to assess which parameters gives the most effective control mechanism here to reach a situation where .
For larger values of ϕ which satisfy the conditions of Proposition 3.8, case B, namely, and
, we observe from numerical experiments that the transcritical bifurcation at
is catastrophic.
Figure (a) demonstrates the stability of the two-strain endemic equilibria for as B varies. A pair of two locally asymptotically unstable equilibria for
appears. Observe that in this range the host-vector model (1) displays traits of an excitable system (Figure (b)). In other words, introduction of a small number of infected individuals in the population leads to a long trajectory in the phase space, which represents a transient epidemic, converging to extinction (orange trajectory in Figure (b)). If the initial condition is in the vicinity of the unstable limit cycle, the system displays an oscillatory behaviour before converging to the disease-free equilibrium (red trajectory in Figure (b)).
Figure 12. (a) One-parameter bifurcation diagram in B for the case when
. Observe that in the range where two interior equilibria occur, both are locally asymptotically unstable. The dot ° marks the value of B where the two-strain equilibrium undergoes a subcritical Hopf bifurcation. (b) Phase portrait for B = 25.6 for two trajectories starting near the two unstable equilibria, marked by a dot °. The orange curve shows an epidemic with a single peak, later converging to extinction. Observe the oscillatory behaviour in the vicinity of the unstable limit cycle (red curve). c Phase portrait for B = 39.979. The trajectory follows an irregular regime.
![Figure 12. (a) One-parameter bifurcation diagram in B for the case B1=B2=B when ϕ=3. Observe that in the range where two interior equilibria occur, both are locally asymptotically unstable. The dot ° marks the value of B where the two-strain equilibrium undergoes a subcritical Hopf bifurcation. (b) Phase portrait for B = 25.6 for two trajectories starting near the two unstable equilibria, marked by a dot °. The orange curve shows an epidemic with a single peak, later converging to extinction. Observe the oscillatory behaviour in the vicinity of the unstable limit cycle (red curve). c Phase portrait for B = 39.979. The trajectory follows an irregular regime.](/cms/asset/19215fec-1faa-4061-8670-1ac787db3d48/tjbd_a_1864038_f0012_oc.jpg)
Note that whenever , the trivial disease-free equilibrium is unstable and the border equilibria
are saddle points as well (Proposition 3.4). The numerical stability analysis shows the endemic equilibrium undergoes a subcritical Hopf bifurcation due to a pair of conjugate complex eigenvalues (marked by the dot in Figure (a)). The arising limit cycle is locally asymptotically unstable, as well as the coexistence equilibrium. Since any trajectory starting within Ω remains within Ω, periodic or chaotic behaviour is possible for larger values of ϕ and values of B such that
. Example of chaotic behaviour is shown in Figure (c).
We recall that the dengue models in [Citation4,Citation18] include multi-strain interactions via ADE but without a temporary cross immunity period. They show deterministic chaos when strong infectivity on secondary infection was assumed: namely, for . In [Citation22] and in the host-vector model (1) considered here, we see a similar behaviour for large ϕ values.
We look next into more detail at the effect of ϕ. In Figure (a), the bifurcation diagram with is shown. These results for the host-only model has been discussed in [Citation22]. Note that the value
corresponds to an average period of temporary cross-immunity of approximately 1 week.
Figure 13. One-parameter diagram for ϕ where and B = 52. Equilibria or maximum and minimum values for limit cycles for total infected I for (a) the host-only model [Citation22] and (b) the host-vector model (1). Red indicates stable and blue unstable solutions.
![Figure 13. One-parameter diagram for ϕ where α=52 and B = 52. Equilibria or maximum and minimum values for limit cycles for total infected I for (a) the host-only model [Citation22] and (b) the host-vector model (1). Red indicates stable and blue unstable solutions.](/cms/asset/23f5d10b-baaf-4d5b-885f-0f6196077754/tjbd_a_1864038_f0013_oc.jpg)
The results for plotted in Figure (b) show that the host-vector model (1) exhibits chaos around
. Here the chaotic region starts at the supercritical Hopf bifurcation
. We remark that the first Lyapunov coefficient at this bifurcation is almost zero because it is close to the Bautin bifurcation GH, Table . This results in the explosion of the amplitude of the resulting dynamics stable which directly governs the chaotic dynamics of the autonomous system near the Hopf bifurcation.
Finally, we present some simulations comparing the reduced models obtained via the geometric singular perturbation theory outlined in Section 4. We explore various parameter regimes for the full model (1), the QSSA model (1) with (Equation9(9a)
(9a) ), and the model (1) with asymptotic expansion of the fast variables
to order ε given by (Equation12
(12a)
(12a) ).
The simulations shown in Figure are for illustrative purposes. The solutions in the full model for the chosen parameter values are symmetric limit cycles (Figure F1 in online Appendix). The limit cycle arises as a result of a Hopf bifurcation from the unique locally asymptotically unstable interior equilibrium . We compare the cycle trajectory of full model to that in the model with QSSA assumption and the approximation using asymptotic expansion for
in Figure . We observe that the limit cycle for the QSSA model has a much larger amplitude in I than the other limit cycles. Further, the QSSA limit cycle in Figure (a–c) is S-conjugate, whereas the limit cycle in the full model is not. The second-order asymptotic expansion in ε (Equation12
(12a)
(12a) ) for
(Figure (a–b)) gives a somewhat better approximation of the limit cycle orbit than the QSSA model, which tends to overshoot significantly in amplitude (compare the time series for
from the model simulations in Figure F1 in online Appendix).
Thus, the simulations demonstrate a trade-off between the level of complexity retained in the QSSA model and the level of accuracy in the approximation of limit cycle orbits. In particular, for the power series approximation of higher-order terms in ε have to be included because even in the simple epidemiological model [Citation27] the first-order asymptotic expansion did not always give a high accuracy. Of course, these numerical experiments serve for illustration of the different approaches to dimension reduction.
5.2. Results for non-autonomous systems
Here we discuss the model under the assumption of a seasonally varying vector population in time which is due to seasonally changing climate factors such as temperature or rainfall that affect the vector's breeding pattern. We model the changing population as periodic forcing (Equation3(3)
(3) ) which adds essentially one new dimension to the system. This means that equilibria become limit cycles, limit cycles dynamics on an invariant torus, either limit cycles, periodic solution or quasiperiodic solutions [Citation25].
Figure (and the enlarged portion near the origin in Figure (b)) for the non-autonomous system (1) and (Equation3(3)
(3) ) look similar to Figure (b) for the autonomous system (1), however, note the torus bifurcation instead of the Hopf bifurcation because of the seasonal forcing of the mosquito population.
Figure 15. Two-strain non-autonomous host-vector model (1) and (Equation3(3)
(3) ) with B = 52. Torus bifurcation TR bifurcation exist at the same position as the Hopf H for the autonomous system (1), shown in Figure (b). L: limit cycle, C: complex dynamics. Details for the range
are shown enlarged in b.
![Figure 15. Two-strain non-autonomous host-vector model (1) and (Equation3(3) M=M(t)=M0(1+ηcosω(t+φ)).(3) ) with B = 52. Torus bifurcation TR bifurcation exist at the same position as the Hopf H for the autonomous system (1), shown in Figure 3(b). L: limit cycle, C: complex dynamics. Details for the range (ϕ,α)∈{(0,1.3)×(0,7)} are shown enlarged in b.](/cms/asset/b64d7d39-3946-4261-9e62-e3119121be82/tjbd_a_1864038_f0015_oc.jpg)
Figure is the one-parameter ϕ diagram where . There is chaos above the torus bifurcation TR. Remarkably the shape of the chaotic attractor changes when the value of the Pitchfork bifurcation P from Figure of the autonomous system is crossed. For the limit cycle case, there the limit cycle goes from a symmetric limit cycle in Figure (a) to S-conjugate limit cycle in Figure (b). There are no periodic solutions but obviously the chaotic solution shows the same change of type of solution of a
-symmetric system.
Figure 16. Two-strain non-autonomous (a) host-only model (A1) and (b) host-vector model (1) with parameter , with
and
respectively. For model (1) in forcing function (3)
and
and in asimilar way for model (A1) in the forcing function in [Citation2], Equation (2)
and
. The bullets in (b) mark the globalmaximum and minimum values for limit cycles for total infected I. Red indicates stable and blue unstable solutions.
![Figure 16. Two-strain non-autonomous (a) host-only model (A1) and (b) host-vector model (1) with parameter α=2, with β=104 and B=52 respectively. For model (1) in forcing function (3) M0=10000 and η=0.1 and in asimilar way for model (A1) in the forcing function in [Citation2], Equation (2) β0=104 and η=0.5. The bullets in (b) mark the globalmaximum and minimum values for limit cycles for total infected I. Red indicates stable and blue unstable solutions.](/cms/asset/2a1028a3-db9f-4be5-89b3-643f32eb2d7e/tjbd_a_1864038_f0016_oc.jpg)
These results from the numerical bifurcation analysis indicate the importance of seasonal fluctuations in the vector population dynamics for , leading to the emergence of a chaotic regime in correspondance with the irregular yearly spikes in dengue haemorrhagic fever incidence reported from Thailand [Citation2]. For the autonomous model, the value
is also associated with chaotic behaviour (Figures and (b)), but in real-world scenarios such values would not be realistic. The value
means that the likelihood in host-to-vector transmission of the dengue virus is equal for the primary and secondary infections, while the value
that can also lead to chaotic behaviour is linked with increased likelihood for secondary infections. A value
would imply that secondary infected hosts are more likely to transmit the virus to mosquitos, yet many such individuals may end up hospitalized or under home observation, and contribute less to the overall force of infection [Citation33].
To summarize, our numerical results demonstrate that the model presented in this work possesses a simpler dynamic structure compared to models in [Citation1,Citation22,Citation23] if hosts with a secondary infection have on average a lower contribution to the overall force of infection upon mosquitos than hosts with a primary infection and if no seasonality of the total vector population is taken into account. Deterministic chaos appears only under the assumption that if hosts with a secondary infection have on average a much higher contribution to the overall force of infection upon mosquitos than hosts with a primary infection. This at first may be surprising due to time series data on dengue haemorrhagic fever incidence reported from Thailand [Citation2] that resemble chaotic behaviour. However, if we include ecological considerations such as seasonal fluctuations into the vector population, the parameter space exhibits regions with deterministic chaos which fits into the context of irregular yearly spikes in the incidence of dengue haemorrhagic fever.
This heuristic approach leads us to the conclusion that inclusion of vector ecology in the model for a vector-borne disease is essential for understanding the complexity and should be considered in models that aim at studying various approaches for control and reduction of the disease burden. Of course, climate parameters do not follow a perfectly regular short-term trend anywhere on Earth, but the longer term climate pattern holds a trend, so the sesonality effect on the vector remains more regular and predictable to a degree, as noted in [Citation19]. Hence, in order to achieve gains in reducing the disease burden, efforts should be directed at examining the interplay of control measures efficacy and the long-term seasonal changes in the vector population. This is important especially for control based on personal protection or treatment of homes and working spaces with targeted release of repellents or insecticides.
6. Conclusion
Often the use of dimension reduction (host-only instead of host-vector) is motivated in order to obtain a model which is numerically better tractable because the time-scale difference between the host and vector yields a stiff ODE system (see [Citation27]). However, use of sufficiently robust numerical techniques allows us to analyse both host-only and host-vector with the same numerical methods.
From the viewpoint of control strategies for vector-borne diseases, we have to incorporate the vector mosquito population in the modelling and to study its impact on the dynamics. In this study, we have included epidemiological factors such as temporary cross-immunity between strains and likelihood of host-to-vector transmission of the dengue virus. In contrast to the model of [Citation15], our model (1) allows for long-term or periodic co-circulation of the two dengue virus strains, which is confirmed by data sets on dengue fever from Thailand.
We have compared the dynamic behaviour of models based on different approaches to dimension and complexity reduction of the full epidemiological model (1): one given by a QSSA for the fast variables and another on an asymptotic expansion based on the theory of geometric singular perturbation. The second-order asymptotic expansion gives a better approximation of the limit cycle orbits for some parameter values than the QSSA model, which tends to overshoot significantly in amplitude (compare the time series for the classes with primary infection from the model simulations in online Appendix Figure F1). When we compare the performance of the asymptotic expansion approximation of this two-strain dengue model to that in the single-strain SIR model [Citation27], we notice that robustness and approximation in each model depends on the order of the truncation. This approach for dimension reduction should be handled carefully in applications after studying the quality of approximation and the complexity involved in computing the higher-order coefficients.
Comparison of bifurcation structures of the full epidemiological model (1) and the host-only model (A1) reveals that the autonomous host-vector system has less complex dynamics. Only stable and unstable equilibria and stable and unstable limit cycles are involved for . Deterministic chaos occurs from a catastrophic transcritical bifurcation when
, but such large values may be difficult to motivate in real-world applications. Torus bifurcations were not found in our study in the host-vector autonomous system (1) in the region we are interested in, and neither was chaos.
For the non-autonomous, host-vector system chaotic dynamics occurs for as in the case of the non-autonomous host-only system via the Ruelle-Takens-Newhouse route to chaos. However, it originates now from a torus bifurcation were the period is fixed by the forcing namely 1 y
while in the host-only model it was free and described by the dynamics on the limit cycle where the torus bifurcation occurs.
We remark that in real-world applications the dengue epidemic dynamics is yearly periodical due to periodicity in the mosquito vector populations. Our analysis leads us to the conclusion that incorporating vector population dynamics in the model for a vector-borne disease is essential for understanding the dynamic complexity and should be an ingredient in models that aim at studying various approaches for control and reduction of the disease burden. However, validation of the models against disease data depends on the depth of the available data sets; in particular, more light needs to be shed on the role of asymptomatic dengue carriers in the epidemiological dynamics because such numbers are not easy to estimate [Citation13]. The bifurcation analysis suggests that with a repellent effect that reduces sufficiently the infection rate per vector, it is possible to stabilise the oscillatory behaviour of the model. This in turn might make the disease more manageable for the healthcare system in the long term as the numbers of infected seeking hospitalization would be less variable in time. As future work, we will use the developed epidemiological model and results from the stability analysis and numerical bifurcation to evaluate the efficacy of control measures against dengue. Such include vaccination and vector control by application of insecticides, repellents and potential release of genetically modified mosquitos [Citation8].
Supplemental Material
Download Zip (503.2 KB)Acknowledgements
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.
Disclosure statement
No potential conflict of interest was reported by the authors.
Additional information
Funding
References
- M. Aguiar, B.W. Kooi, and N. Stollenwerk, Epidemiology of dengue fever: a model with temporary cross-immunity and possible secondary infection shows bifurcations and chaotic behaviour in wide parameter regions, Math Model Nat. Phenom. 3 (2008), pp. 48–70.
- M. Aguiar, S. Ballesteros, B.W. Kooi, and N. Stollenwerk, The role of seasonality and import in a minimalistic multi-strain dengue model capturing differences between primary and secondary infections: complex dynamics and its implications for data analysis, J. Theor. Biol. 289 (2011), pp. 181–196.
- S. Bedini, G. Flamini, R. Ascrizzi, F. Venturi, G. Ferroni, A. Bader, J. Girardi, and B. Conti, Essential oils sensory quality and their bioactivity against the mosquito Aedes albopictus, Sci. Rep. 8 (2018), pp. 17857.
- L. Billings, I.B. Schwartz, L.B. Shaw, M. McCrary, D.S. Burke, and D.A.T. Cummings, Instabilities in multiserotype disease models with antibody-dependent enhancement, J. Theor. Biol. 246 (2007), pp. 18–27.
- B. Buonomo and R. DellaMarca, Optimal bed net use for a dengue disease model with mosquito seasonal pattern, Math. Meth. Appl. Sci. 41(2) (2018), pp. 573–592.
- L.B. Carrington and C.P. Simmons, Human to mosquito transmission of dengue viruses, Front. Immunol. 5 (2014), pp. 290.
- K.H. Changal, A.H. Raina, A. Raina, M. Raina, R. Bashir, M. Latief, T. Mir, and Q.H. Changal, Differentiating secondary from primary dengue using IgG to IgM ratio in early dengue: an observational hospital based clinico-serological study from north india, BMC Infect. Dis. 16 (2016), pp. 715.
- M.R.W. de Valdez, D. Nimmo, J. Betz, H.-F. Gong, A.A. James, L. Alphey, and W.C. Black IV, Genetic elimination of dengue vector mosquitoes, Proc. Natl Acad. Sci. USA 108(12) (2011), pp. 4772–4775.
- O. Diekmann and J.A.P. Heesterbeek, Mathematical Epidemiology of Infectious Diseases, Wiley series in Mathematical and Computational Biology. Wiley, Chichester, 2000.
- O. Diekmann, J.A.P. Heesterbeek, and J.A.J. Metz, On the definition and the computation of the basic reproduction ratio r0 in models for infectious diseases in heterogeneous populations, J. Math. Biol.28 (1990), pp. 365–382.
- O. Diekmann, M. de Jong, A. de Koeijer, and P. Reijnders, The force of infection in populations of varying size: a modelling problem. Technical Report AM-R9403, Centrum voor Wiskunde en Informatica, 1994.
- E.J. Doedel and B. Oldeman, AUTO 07p: Continuation and bifurcation software for ordinary differential equations, Technical report, Concordia University, Montréal, Canada, 2012.
- V. Duong, L. Lambrechts, R.E. Paul, S. Ly, R.S. Lay, K.C. Long, R. Huy, A. Tarantola, T.W. Scott, A. Sakuntabhai, and P. Buchya, Asymptomatic humans transmit dengue virus to mosquitoes, Proc. Natl Acad. Sci. USA 11 (2015), pp. 14688–14693.
- L. Esteva and C. Vargas, Analysis of a dengue disease transmission model, Math. Biosci. 150(2) (1998), pp. 131–151.
- Z. Feng and J.X. Velasco-Hernandez, Competitive exclusion in a vector-host model for the dengue fever, J. Math. Biol. 35 (1997), pp. 523–544.
- N. Fenichel, Persistence and smoothness of invariant manifolds for flows, Indiana. U Math. J. 21 (1971), pp. 193–226.
- N. Fenichel, Geometric singular perturbation theory, J. Differ. Equ. 31 (1979), pp. 53–98.
- N. Ferguson, R. Anderson, and S. Gupta, The effect of antibody-dependent enhancement on the transmission dynamics and persistence of multiple-strain pathogens, Proc. Natl Acad. Sci. USA 96(9) (1999), pp. 790–794.
- T. Götz, K.P. Wijaya, and E. Venturino, Introducing seasonality in an SIR-UV epidemic model: an application to dengue. 18th International Conference on Computational and Mathematical Methods in Science and Engineering, CMMSE 2018, Wiley Online Library, 2018
- M.G. Guzman, M. Alvarez, and S.B. Halstead, Secondary infection as a risk factor for dengue hemorrhagic fever/dengue shock syndrome: an historical perspective and role of antibody-dependent enhancement of infection, Arch. Virol. 158 (2013), pp. 1445–59.
- G. Hek, Geometric singular perturbation theory in biological practice, J. Math. Biol. 60 (2010), pp. 347–386.
- B.W. Kooi, M. Aguiar, and N. Stollenwerk, Bifurcation analysis of a family of multi-strain epidemiology models, J. Comput. Appl. Math. 252 (2013), pp. 148–158.
- B.W. Kooi, M. Aguiar, and N. Stollenwerk, Analysis of an asymmetric two-strain dengue model, Math. Biosci. 248 (2014), pp. 128–139.
- A. Korobeinikov, Lyapunov functions and global properties for SEIR and SEIS epidemic models, Math. Med. Biol. 21 (2004), pp. 75–83.
- Y.A. Kuznetsov, Elements of Applied Bifurcation Theory, 3rd ed., Vol. 112 of Applied Mathematical Sciences. Springer-Verlag, New York, 2004.
- N.M. Nguyen, D. Thi Hue Kien, T.V. Tuan, N.T.H. Quyen, C.N.B. Tran, L. Vo Thi, D.L. Thi, H.L.Nguyen, J.J. Farrar, E.C. Holmes, M.A. Rabaa, J.E. Bryant, T.T. Nguyen, H.T.C. Nguyen, L.T.H.Nguyen, M.P. Pham, H.T. Nguyen, T.T.H. Luong, B. Wills, C.V.V. Nguyen, M. Wolbers, and C.P.Simmons, Host and viral features of human dengue cases shape the population of infected and infectious Aedes aegypti mosquitoes, Proc. Natl Acad. Sci. USA 110 (2013), pp. 9072–9077.
- P. Rashkov, E. Venturino, M. Aguiar, N. Stollenwerk, and B.W. Kooi, On the role of vector modeling in a minimalistic epidemic model, Math. Biosci. Eng. 5(16) (2019), pp. 4314–4338.
- M. Recker, K.B. Blyuss, C.P. Simmons, T.T. Hien, B. Wills, J. Farrar, and S. Gupta, Immunological serotype interactions and their effect on the epidemiological pattern of dengue, Proc. R. Soc. B276(1667) (2009), pp. 2541–2548.
- F. Rocha, M. Aguiar, M. Souza, and N. Stollenwerk, Time-scale separation and centre manifold analysis describing vector-borne disease dynamics, Int. J. Comput. Math. 90(10) (2013), pp. 2105–2125.
- H.S. Rodrigues, M.T.T. Monteiro, D.F.M. Torres, A.C. Silva, C. Sousa, and C. Conceio, Dengue in Madeira island. In Dynamics, Games and Science -- International Conference and Advanced School Planet Earth, J.-P. Bourguignon, R. Jeltsch, A. A. Pinto, and M. Viana, eds., DGS II, Portugal, pp. 593–605, August 28–September 6, 2013.
- H.S. Rodrigues, M.T.T. Monteiro, and D.F.M. Torres, Seasonality effects on dengue: basic reproduction number, sensitivity analysis and optimal control, Math. Meth. Appl. Sci. 39(16) (2016), pp. 4671–4679.
- D. Ruelle and F. Takens, On the nature of turbulence, Comm. Math. Phys. 20 (1971), pp. 167–192.
- Q.A. ten Bosch, H.E. Clapham, L. Lambrechts, V. Duong, P. Buchy, B.M. Althouse, A.L. Lloyd, L.A. Waller, A.C. Morrison, U. Kitron, G.M. Vazquez-Prokopec, T.W. Scott, and T.A. Perkins, Contributions from the silent majority dominate dengue virus transmission, PLoS Pathog. 14(5) (2018), p. e1006965.
- H. Tian, Z. Sun, N.R. Faria, J. Yang, B. Cazelles, S. Huang, B. Xu, Q. Yang, O.G. Pybus, and B. Xu, Increasing airline travel may facilitate co-circulation of multiple dengue virus serotypes in Asia, PLoS Negl. Trop. Dis. 11(8) (2017), p. e0005694.
- WHO, Fact sheet: dengue and severe dengue. Available at https://www.who.int/. Accessed 4 November 2019.
- T.-T. Zheng and L.-F. Nie, Modelling the transmission dynamics of two-strain dengue in the presence awareness and vector control, J. Theor. Biol. 443 (2018), pp. 82–91.
Appendices
Appendix 1. Host-only Model
The two-strain host-only model from [Citation1,Citation2] is given by the system
(A1a)
(A1a)
(A1b)
(A1b)
(A1c)
(A1c)
(A1d)
(A1d)
(A1e)
(A1e)
(A1f)
(A1f)
(A1g)
(A1g)
(A1h)
(A1h)
(A1i)
(A1i)
(A1j)
(A1j)
Appendix 2. Solving for the exclusion equilibrium ![](//:0)
![](//:0)
Setting the right-hand side of (1) equal to 0 yields
Further, we have the following linear system in
,
which can be solved explicitly via algebraic manipulation and gives the corresponding equilibrium values for
.
Appendix 3. Jacobian matrices
The Jacobian matrices used in the local stability analysis of steady states are
(A3)
(A3)
(A4)
(A4) with
.
Appendix 4. Proofs
A.1. Proof of proposition 3.1
The Jacobian of the right-hand side of (1) (see Appendix) evaluated at
has eigenvalues
(counted with multiplicity) as well as the eigenvalues of the
-minor
which are the roots λ of the polynomial
(A5)
(A5) Polynomial (EquationA5
(A5)
(A5) ) has a positive real root if the free term is negative, a condition equivalent to
. Hence,
is locally asymptotically unstable if the above condition is met.
A.2. Proofs of proposition 3.6 and 3.7
Let . Then at the interior equilibrium
(A6a)
(A6a)
(A6b)
(A6b)
(A6c)
(A6c)
(A6d)
(A6d)
(A6e)
(A6e)
(A6f)
(A6f)
(A6g)
(A6g)
Substituting the equilibrium values for the infected classes into the vectors (Equation1j(1j)
(1j) )–(Equation1k
(1k)
(1k) ), we have to solve the nonlinear system
Dividing both sides by
, respectively, and setting
we arrive to the equivalent system
and since the right side is not identically zero, if there exist positive solutions
, it must hold that
or
when
. Hence, we have shown that the interior equilibrium
, whenever it exists, is a fixed equilibrium for
. The values of the respective classes are
etc.
If , the above system is under-determined in
, and we solve for
and the remaining classes can be expressed accordingly using (EquationA6
(A6a)
(A6a) ). This determines the curve
.
Else, for a substitution into the equation for the vector (Equation1j
(1j)
(1j) ) yields for x:
We must solve, equivalently, the quadratic
(A7)
(A7) for solutions x>0, whose coefficients are given by
We examine the number of positive roots of (EquationA7
(A7)
(A7) ).
Because the leading coefficient of (EquationA7
(A7)
(A7) ) is always negative, Vieta's rule tells us the number of positive roots depends on the signs of the other terms. In particular, Equation (EquationA7
(A7)
(A7) ) has exactly one positive root if and only if
Hence, it is sufficient to have
. We express the relations for the signs of the coefficients
as inequalities for
.
Whenever , or
, we have
, and otherwise,
if
.
Equation (EquationA7(A7)
(A7) ) has exactly two positive roots if and only if
We note that if
for all
. Otherwise when
holds.
The discriminant of (EquationA7(A7)
(A7) )
if the following conditions are satisfied
giving the boundary between the regions
and
.
To show local stability of the curve (for the case
), we compute the the Jacobian at a point
:
where the entries satisfy the following:
, such that
and
are computed from (EquationA6
(A6a)
(A6a) ).
From 's structure it is obvious that
are eigenvalues (the first two have algebraic multiplicity 2). The characteristic polynomial of the remaining
-minor of
is
where for the sake of shortness, we denote
. This polynomial has only non-negative coefficients, and does not have any positive real roots. Hence, the curve
consists of points which are local attractors.
A.3. Proof of lemma 3.9
Note
for the range
. Hence, the existing solution x must be monotonically increasing with ϕ to maintain the equality in
since the solution is unique.
Since , we can approximate the solution of the equation by replacing x by
into the numerator of the term on left-hand side of the equation. Then
giving
(A8)
(A8) This means that the interior equilibrium value for the infected vector classes
.
A.4. Normal hyperbolicity of the set of critical points
Recall the fast system
The Jacobian
of the fast system in
on the set of critical points is given by
This matrix is constant on the positive half-lines
, and the eigenvalues are
are bounded uniformly away from the imaginary axis, which concludes the proof.
Appendix 5. Asymptotic expansion
The substitution of the power series in ε into the invariance equation (Equation11(11)
(11) ) gives
Taking the expansions of f, g up to order ε produces the following system:
(A9)
(A9) Collecting terms of order
gives the system
whose solutions are the expressions of the free terms in (Equation12
(12a)
(12a) ):
(A10a)
(A10a)
(A10b)
(A10b)
To find the coefficients for terms of order in (EquationA9
(A9)
(A9) ), we use the fact that many partial derivatives vanish identically
which results in the vector
to be used in Equation (Equation13
(13)
(13) ), and where the partial derivatives are given by
The coefficients
are given by
(A11a)
(A11a)
(A11b)
(A11b)
For the second-order term, the vector has entries
where the partial derivatives are computed according to the expressions for
.
are computed from (Equation13
(13)
(13) ).