Effects of complexity and seasonality on backward bifurcation in vector-host models.

We study implications of complexity and seasonality in vector-host epidemiological models exhibiting backward bifurcation. Vector-host diseases represent complex infection systems that can vary in the transmission processes and population stages involved in disease progression. Seasonal fluctuations in external forcing factors can also interact in a complex way with internal host factors to govern the transmission dynamics. In backward bifurcation, the insufficiency of R0 < 1 for predicting the stability of the disease-free equilibrium (DFE) state arises due to existence of bistability (coexisting DFE and endemic equilibria) for a range of R0 values below one. Here we report that this region of bistability decreases with increasing complexity of vector-borne disease transmission as well as with increasing seasonality strength. The decreases in the bistability region are accompanied by a reduced force of infection acting on primary hosts. As a consequence, we show counterintuitively that a more complex vector-borne disease may be easier to control in settings of high seasonality.

Introduction

Predicting the elimination/extinction of an infectious disease from a population is an important question that confronts policy makers in the event of a disease outbreak. Microparasitic vector borne diseases such as West Nile virus (WNV), chikungunya, dengue, malaria and Zika differ in complexity, are prevalent or are emergent in different parts of the world [1-10]. The transmission cycles of these diseases are maintained in various hosts such as birds, humans and other animals, with the degree of complexity differing between diseases in terms of the stages involved in the progression of disease in each host population, and specifics of the structure and processes governing the transmission process between hosts [11,12].

The essential features of the transmission dynamics of microparasitic diseases can be studied via compartmental models, whereby the populations of hosts and vectors are divided into distinct classes or compartments, each defining the state of an individual with respect to the stage of disease progression, e.g. susceptible, exposed, infectious, recovered hosts [11]. Increasing the number of compartments and transmission structures in such a system increases realism, but also behavioural complexity in these models [13,14]. For these models, the basic reproduction ratio (R0), defined as the number of secondary cases produced by a single infected case, can serve as an indicator of disease endemicity in the population. In the simplest case (e.g. the susceptible-infected-recovered (SIR) model), the condition R0 = 1 defines a bifurcation (or critical) point, whereby when R0 > 1 a stable infected equilibrium emerges and establishes, whereas when it is less than unity the disease-free equilibrium (DFE) emerges and is generally globally stable [11].

By contrast, recently it has been shown that compartmental models may describe a scenario in which a turning point of the infected equilibrium may exist in a region where all equilibrium states are positive and R0 < 1 [15-18]. This induces multiple stable equilibria and disruption of the global stability of the disease-free equilibrium, with the result that instead of converging globally to the disease free state when R0 < 1, the solution may approach an infected state that coexists with the disease-free state. Several disease models have been investigated for the existence of this phenomenon in the literature [1,5,15,19-23], and compartmental models developed for vector-borne diseases, including WNV [1,19,24], dengue [21], malaria [5], and other arboviral diseases [22] have also been shown to exhibit backward bifurcation and thus violation of the R0 = 1 threshold condition separating diseases endemicity and extinction [23].

Several mechanisms leading to backward bifurcation have been proposed, such as partially effective vaccination programmes [25,26], social differences in susceptibles [17], nonlinear incidences [15], the structure of interactions among multi-groups [27-29], and multiple stages of infection [30]. However two features appear important: a nonlinear incidence rate due to saturation effects [13,31]; and asymmetry in the susceptibility of hosts to infection [19], including in the death rates of susceptible and infected hosts. The latter disease induced increase in death rates among infected hosts along with a mass action force of infection has been shown to be a particularly important driving factor for the occurrence of backward bifurcation in WNV, dengue, malaria and other arboviral models [5,19,20].

Seasonal variation in reported cases of infectious diseases, ranging from childhood diseases (e.g. measles, diphtheria and chickenpox) to infections affecting all ages like influenza, cholera and vector borne diseases (e.g. WNV, malaria, dengue), is common across temperate and tropical geographies [32-34]. In compartmental models for malaria, it has been shown that periodic changes in mosquito populations can reduce the basic reproductive ratio [35-37]. Similarly, the recurrence of WNV in many previously affected areas and incursion into new geographical areas is also, in most cases, related to seasonal effects on vector population dynamics [24]. There is a growing body of literature [11,12,32-34,38] investigating the role of external drivers (e.g. rainfall, temperature) responsible for seasonal fluctuations in disease incidence, but few studies exist that have investigated the impact of these important external drivers in diseases suspected to show backward bifurcation [39].

In this paper, we investigate the dynamics of a general host susceptible-exposed-infective-temporarily recovered (SEIRS)–vector SEI model, and tune it to study the effect that complexity (in terms of structure, topological structure of interactions and parameters) may have on the phenomenon of backward bifurcation. We also study the influence of seasonality on the phenomenon using the model. To ground the analysis, we investigated the impact of both factors using parameter regimes derived from the literature for WNV, dengue and malaria.

The subsequent sections are organized as follows. In the upcoming section we first describe the SEIRS–SEI–SEIRS model and derive a quadratic equation for endemic equilibrium where a dead-end host is absent. In the results section, we evaluate and present the conditions for backward bifurcation to exist in the presence of seasonality in our disease systems, and compare them with the non-seasonal case. We then investigate the effect of disease complexity by adding compartments beginning from the stage of the model (see below), and follow this by quantifying the effect of including a dead-end host. We end by discussing the likely occurrence, extent, and the factors that underlie backward bifurcation in our disease systems, and the implications they have for disease control.

Material and methods

General representation of compartmental models

A typical compartmental model describes infection progression between host classes representing susceptibles, exposed, infectives and recovered individuals [11,12]. During an emerging outbreak, all individuals are susceptible; once infected they can remain dormant (exposed) for some time before infecting others or they immediately become infectives and start infecting others. Infected individuals may (i) recover, (ii) die due to infection, or (iii) become susceptible again.

In general, in this system, a heterogeneous population can thus be grouped into N homogeneous compartments. We can rearrange the compartments such that the first m of them correspond to the infected (including symptomatic and asymptomatic) while the rest correspond to uninfected classes. The time evolution of the state of the population {x} has the general form: where represents input rate of newly infected individuals in the ith compartment and with being the transfer in/out of the compartment i by all other means. Writing the compartmental model in the form of equation (2.1) facilitates the calculation of the basic reproductive ratio R0 for complex models involving multiple hosts or vectors [40] as shown in the electronic supplementary material, appendix. This framework also allows development of a hierarchical set of specific models of increasing or decreasing complexity. Below we begin by describing a maximal (SEIRS–SEI–SEIRS) model for these systems, reflecting the dynamics of transmission represented by WNV. We then show how complexity may be reduced to capture the dynamics of dengue and malaria transmission.

The seasonal SEIRS–SEI–SEIRS model

A one-vector–two-host system where the second host is a dead end is characterized by a network of three nodes where the vector node infects the two host nodes but gets infected by only one host node, as shown in figure 1. The mutually infecting vector–host nodes form the primary cycle of the disease. Here infection progression in each host node is described by a susceptible (S), exposed (E), infective (I), temporarily recovered (R) (SEIRS) model while the vector node is described by a susceptible (MS), exposed (ME), infective (MI) (SEI) model assuming that vectors do not recover from the disease. The equations for the model are as follows: where the state variables correspond to the primary host-1, correspond to secondary host-2, correspond to vectors, are force of infection on host 1, host 2, and vectors respectively. The total number of mosquito bites b are divided between hosts 1 and host 2 and are given by b1 and b2. The definition of other parameters are given in table 1. Seasonal variations are introduced by a sinusoidal mosquito birth rate as discussed in the electronic supplementary material, appendix. In all our later calculations we will assume . The model equations (2.2)–(2.4) have a disease free equilibrium given by: Definition of R0 in presence of seasonality is not as straightforward. The linear operator method [41,42] can be used to obtain an accurate estimate of R0; however in most cases it involves numerical solutions. Since the total population of mosquitoes is governed by , at equilibrium this is equivalent to a sinusoidally varying population and opens up the methods of Bacar [37] to be applied for calculating R0 as a function of seasonality strength ε, i.e. R0(ε). The method of Bacar [37] is based on using a combination of the survival function approach and Fourier analysis. It is approximate but provides a deeper analytical understanding of the effects of seasonality. We give a summary of both these procedures in electronic supplementary material, appendix, but use the later method unless specified otherwise. The expression for R0(ε) is obtained as: where R0 is the basic reproduction ratio, in the absence of seasonality, as obtained using the survival function approach [43], and is the oscillation frequency, T  =  365 is the period of oscillation in days and the function G is obtained using the methods of Bacar [37]. For the parameters of diseases considered here, equation (2.10) leads to the condition .

Figure 1.

A one-vector–two-host system where the second host is a dead end. This is equivalent to a network of three nodes where the vector node infects the two host nodes but gets infected by only one host node. (a) The host--vector interaction network. (b) Flow chart for the host--vector interaction network.

Table 1.

Parameters and their interpretation in the full SEIRS–SEI model. A blank (–) unit means the parameter is dimensionless.

parameter	interpretation	units
b	average mosquito biting rate	—
b1	mosquito biting rate on host 1	day−1
b2	mosquito biting rate on host 2	day−1
β_1	transmission probability host 1 to mosquito	—
β_H1	transmission probability mosquito to host 1	—
β_H2	transmission probability mosquito to host 2	—
Π_M	mosquito birth rate	day−1
μ_M	mosquito death rate	day−1
Π_H1	host 1 birth rate	day−1
Π_H2	host 2 birth rate	day−1
μ_H1	host 1 death rate	day−1
μ_H2	host 2 death rate	day−1
dM	disease induced death rate mosquito	day−1
d_H1	disease induced death rate in hosts 1	day−1
d_H2	disease induced death rate in hosts 2	day−1
1σ_M	exposed period mosquito	day
1σ_H1	exposed period in hosts 1	day
1σ_H2	exposed period in hosts 2	day
τ_H1	host 1 recovery rate	day−1
τ_H2	host 2 recovery rate	day−1
α_H1	host 1 loss of immunity rate	day−1
α_H2	host 2 loss of immunity rate	day−1
η_M	modification parameter mosquito	—
η_H1	modification parameter host 1	—

The endemic equilibrium of the equations (2.2)–(2.4) is given by: where and For the case of small secondary host population NH2 ≈ 0, then inserting equation (2.13) and equation (2.15) into the forces of infection equations (2.5) and (2.7), it is easy to see that λH1 satisfies the quadratic equation It is reasonable to assume that a quadratic equation of the above type continues to exist even when is finite and is confirmed in numerical simulations, which we discuss below.

Complexity

Complexity of a model can be broadly divided into two categories: To illustrate the second form of complexity we take the example of just the primary host–vector cycle (): setting implies that there are no exposed and recovered classes and hence the model is obtained (WNV [1,18]), setting the model (malaria [5]) is obtained, and reduces equation (2.2) to model (dengue one serotype [20,21]). The first type of complexity is introduced by additional host/vector species and in our case we add a dead-end host. The presence of dead-end hosts such as humans and other mammals in the vicinity of birds has relevance for WNV [24].

The number of distinct host/vector species affecting the disease spread. This corresponds to distinct nodes as shown in figure 1a.

Each node is represented by a compartmental model: the complexity can thus also increase as the number of compartments and/or parameters governing the transfer from one compartment to another is increased.

Although transmission complexity can be introduced in a variety of ways, e.g. nonlinear transmission rate arising due to saturation effects, in this work we consider a linear transmission rate to minimize system nonlinearity, allowing us to investigate complexity in terms of increasing compartments and hosts.

Sensitivity index

We investigated the effect of parameters on the backward bifurcation critical point Rc using sensitivity analysis. The sensitivity index () of a variable u with respect to the parameter p is formally defined as the ratio of relative change in a variable u to the relative change in the parameter p [44]:

Results

In light of equations (2.16)–(2.19), we can state the equivalent of Theorem 1 in (3.1) that the primary host–vector (SEIRS–SEI) model has: Using the condition the critical value of for backward bifurcation, denoted by , is then obtained as: where A and B are functions of parameters except the seasonality strength (ε), see electronic supplementary material, appendix. The range is the backward bifurcation region.

a unique endemic equilibrium if (as c0 < 0)

a unique endemic equilibrium if b0 < 0 and c0 = 0 or

two endemic equilibria if c0 > 0, b0 < 0 and

no endemic equilibrium otherwise.

Bifurcation diagrams in this paper are obtained using the Newton–Raphson root finding procedure for stable and unstable fixed points of a dynamical system [45] (unless mentioned otherwise). In the absence of seasonality, i.e. ε = 0, a bifurcation diagram in figure 2 for the approximation of the model shows that the endemic state coexists with the DFE (black horizontal line) in the backward bifurcation region . Clearly, since is not a sufficient criterion for disease elimination, it is instructive to look at the backward bifurcation region in the space, since and can be changed by mosquito control measures. In the parameter space the area between the curves C1 and CRC represents the backward bifurcation region as shown in figure 3a for the approximation of the model. As the complexity of the model is increased by bringing into play more compartments/parameters, the backward bifurcation region is altered in significant ways: for example, a simple inclusion of temporary immunity shifts the curves C1 and CRC so that a larger value of (for a given ) is required to be in the backward bifurcation region. Activating an exposed/latent state increases this threshold further, as shown in figure 3b,c. The shift in Rc with addition of more compartments is easily seen in the space where it translates into an increase in Rc as shown in figure 3d. In general, a compartment representing the delay of infectious state or a temporary (permanent) immunity in hosts increases the threshold Rc w.r.t the basic model. The sensitivity indices (see equation (2.20)) of Rc w.r.t parameters (p) are tabulated in table 2. A positive (or negative) value of indicates that an increase in parameter p also increases (or decreases) Rc. According to table 2 the parameters of the model clearly form two groups based on whether they increase or decrease the value of Rc, assuming an increase in these parameters. Increasing the values in the first group increases while increasing the values of the second group decreases the threshold Rc. Counterintuitively, do not change Rc, while a weakly positive but negligible dependence on is observed. The parameter regimes of WNV, dengue, malaria (table 3) produce different values of Rc in the space, as shown in figure 4, therefore the relative change in Rc is calculated from their respective base values of Rc. As an example, for WNV and , which means that a 6.4% change in produces a 10% change in the value of Rc.

Figure 2.

Bifurcation diagram as a function of R0 in the absence of seasonality. Here the top black and red curves are the force of infection corresponding to the stable and unstable endemic states of host 1 when R0 ≤ 1. The black horizontal line is the disease-free equilibrium (DFE). At R0 ≥ 1 the DFE becomes unstable and only one stable endemic state exists. The parameters for WNV from table 3 were used in the SEIRS–SEI–SEIRS model, with no secondary hosts, to generate this figure.

Figure 3.

The model (without seasonality) and effects of complexity. The backward bifurcation region in the space is sandwiched between and curves which correspond to and respectively. The plots in (a–c) show that as Rc increases due to increasing complexity it also shifts the C1 and curves thus allowing for better control by varying the mosquito birth/death rates. In (d) the shift of Rc curve is shown in the space as more compartments are added to the model. Here the WNV model parameters (table 3) were used in the approximation model and were used to change the complexity.

Table 2.

Sensitivity index w.r.t model parameters (p) for WNV, dengue and malaria. We took while calculating and for the rest .

parameter (p)	WNV	dengue	malaria
b1	0.08394292488	0.0002338511190	0.0009617256570
β_H1	0	0	0
β_1	0.08394292485	0.0002338511190	0.0009617256566
μ_M	−0.08394292425	−0.0002338534960	−0.0009617250206
Π_H1	0	0	0
μ_H1	0.6414893380	0.0009251279255	0.1790822307
d_H1	−0.6415417466	−0.5655065196	−0.8441519805
σ_M	0	0	0
σ_H1	−0.3345086117 × 10−5	−0.0006105349835	−0.006315525812
τ_H1	0	0.5651919267	0.8427761814
α_H1	0	0	−0.1713909052
ε	0.002159482349	0.0004596099794	0.00007150530087

Table 3.

Model parameters for WNV, dengue and malaria. Model parameters for WNV, dengue and malaria are based on the following references: [1,19], ([21], E Michael 2015, unpublished work) and [44]. The units of parameters are given in table 1.

parameter	WNV	dengue	malaria
b	0.1	0.315	0.42
β_1	0.16	0.9	0.64
β_H1	0.88	0.8	0.08
β_H2	0.88	0	0
Π_M	variable	variable	variable
μ_M	0.0625	0.0625	0.0625
Π_H1	1000	1000	1000
Π_H2	variable	0	0
μ_H1	0.001	170×365	170×365
μ_H2	170×365	—	—
d_M	0	0	0
d_H1	0.005	0.99	0.99
d_H2	10−5	—	—
σ_M	1000	0.0904	0.0909
σ_H1	1000	0.1671	0.0714
σ_H2	0.0714	—	—
τ_H1	0	0.274	0.1
τ_H2	0.0714	—	—
α_H1	0	0	0.00001
α_H2	0	—	—
η_M	0	0	0
η_H1	0	0	0

Figure 4.

(a) The backward bifurcation boundary Rc in the space for the model parameters corresponding to WNV (black), dengue (red) and malaria (green) from table 3. Although Rc varies with in all three cases the range of variation is largest for WNV parameters. The backward bifurcation region is largest for malaria followed by WNV and dengue when (assuming a mosquito lifespan of 16 days). (b,c) Bifurcation diagrams showing stable and unstable branches of (b) force of infection as a function of R0 for malaria and dengue and (c) force of infection for WNV model with along with an model with a small recovery rate (cf. Figure 3d). We plotted dengue and malaria in the same graph because of similarity in model dimension and parameters in contrast to WNV. Even though the WNV shows large force of infection as compared to malaria and dengue, this may be because no latent period was considered for WNV along with incomparable natural and disease-induced host death rates. However, the general principle that a larger force of infection increases the backward bifurcation region in both (b) and (c) can be seen to hold.

When the dead-end host is also integrated into the model, as the number of secondary hosts (given by ) is increased the backward bifurcation region decreases, as shown in figure 5 for the WNV model parameters listed in table 3. This is because involving secondary dead-end hosts will lead to a reduction in the overall force of infection (figure 5).

Figure 5.

WNV model with humans as dead-end hosts. A bifurcation diagram showing stable (black and green) and unstable (red and blue) branches of the force of infection as a function of dead-end host population (governed by the birth rate ). It shows that the force of infection reduces as the number of bites on the dead-end hosts increases (see equations (2.5)--(2.8)).

When seasonality is switched on (), it is easy to see from equation (3.1) that the threshold Rc increases with increasing seasonality strength. A representative bifurcation diagram shown in figure 6 illustrates this result. The bifurcation diagram was obtained by replacing with in the nonseasonal equations. This substitution simplifies the calculation of the bifurcation diagram and is justified from the following two observations: (i) R0 in equation (2.12) (non-seasonal case) then becomes equal to R0(ε) in equation (2.10) (seasonal case), and (ii) a bifurcation diagram using the linear operator method for R0(ε) with the original seasonal term also shows trends similar to figure 6 with increasing seasonality (see electronic supplementary material, figure S1). It also follows, as a consequence of equation (3.1), that the threshold mosquito birth rate is higher in the seasonal case as compared to the non-seasonal case. The sensitivity index in table 2 confirms this result and it also shows that the WNV parameters lead to a stronger change in Rc when compared to dengue/malaria for the seasonal case.

Figure 6.

A bifurcation diagram showing stable (colours as labelled) and unstable (violet) branches of the force of infection as a function of ε. It shows that the force of infection reduces, resulting in increase of Rc, as ε is increased. To obtain this figure we use the non-seasonal equations and replace by for a given R0 (equation (2.12)) in the model equation (2.2) and then evaluate the stable and unstable branches using the Newton–Raphson method (see text for justification of this procedure).

We also find that the reduced/increased backward bifurcation region for all the cases above is always accompanied by a reduced/increased force of infection on primary hosts . This is easily deduced from the bifurcation diagrams in figures 4–6.

Discussion

The condition R0 < 1 for stability of DFE breaks down due to the backward bifurcation phenomenon. It has been shown to exist in models of many vector-borne diseases [1,21,24]. However, previous studies ignored the impact of increasing complexity, including seasonality.

In this work we argued that the approximation of the model is a basic vector–host model exhibiting backward bifurcation. Then the impact of increasing complexity and seasonality on Rc was studied. For instance the addition of exposed or immune compartments in hosts changes the model to an or model and in both cases the threshold Rc is found to increase, with the addition of recovereds found to have the stronger effect. On the other hand, if immunity is temporary in the model then the resulting model has a relatively lower threshold value Rc though still higher than the original model. Therefore, in general, increasing model complexity (by adding delay to the infectious state or a temporary (permanent) immunity in hosts) leads to an increase in the threshold Rc w.r.t the basic model. Similarly, the addition of a dead-end host in the model also reduces the backward bifurcation region for the simple reason that if mosquitoes bite dead-end hosts more, such bites will result in reducing the overall disease transmission rate.

For the parameter regimes of WNV [1,19], malaria [5,44] and dengue ([21], E Michael 2015, unpublished work) the sensitivity analyses carried out in this study show that the host (bird) disease-induced death rate strongly affects the backward bifurcation region for all diseases with increases in death rates of infecteds under conditions of variable birth rates increasing the backward bifurcation region. However, for the WNV case, these results also show that increasing the natural death of hosts (birds) more markedly increases the threshold Rc. This indicates that culling of hosts (birds) could constitute the major strategy for managing WNV spread in areas where backward bifurcation is possible for this disease. By contrast, because increased immunity will decrease the backward bifurcation regions more for dengue and malaria (table 2), improving the rate of host recovery by vaccination may be a more effective control strategy for these diseases.

Seasonality always reduces the value of R0 resulting in an increase in the value of Rc; this implies that it is easier to control a vector-borne disease in highly seasonal areas relative to weakly seasonal ones. Based on our results we observe that although the backward bifurcation region decreases in the order of malaria, WNV and dengue, the results from the sensitivity analysis predict that seasonality has stronger effect on WNV followed by dengue and malaria.

In general, however, we note that given the observation that an increase in Rc values is a consequence of a reduced force of infection acting on primary hosts, this effect could be used to develop a guiding principle for designing complexity-based control programmes for vector-borne diseases. For example, this could involve, given the specificity of the disease, culling of primary hosts such as birds in the case of WNV, adding of secondary hosts in malaria (zooprophylaxis) or vaccination of humans (as proposed for dengue and malaria).

The simple sinusoidal variation in mosquito birth rates studied in this work may not capture the full complexity of seasonality; incorporating the exact form of this factor has been pointed out as a significant element to be considered in disease modelling (see the arguments put forth in [38,46]), and growing empirical evidence indicates that these external forcings (i.e. via rainfall, temperature fluctuations acting on vector populations as well as pathogen development rates) may indeed vary year-to-year [47,48]. However, based on the fact that any oscillatory signal could be expanded in a Fourier series with multiple frequencies, it is reasonable to assume that calculations of R0(ε) reduce to the results of Bacar [37] when all but one frequency has a dominant role.

Incorporating nonlinear transmission rates (depicting saturation effects among others), and variable rather than fixed parameter range also increase the complexity of the models. An investigation of nonlinear forms of the transmission rate, multiple host–vector models, seasonality reflected via multiple frequencies with different amplitudes, and variable parameter ranges to gain a fuller understanding of process complexity in these systems clearly form a task for future work.