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][2][3][4][5][6][7][8][9][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 (R 0 ), 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 R 0 = 1 defines a bifurcation (or critical) point, whereby when R 0 > 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 R 0 < 1 [15][16][17][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 R 0 < 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][20][21][22][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 R 0 = 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][28][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][33][34]. In compartmental models for malaria, it has been shown that periodic changes in mosquito populations can reduce the basic reproductive ratio [35][36][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][33][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 SI − M S M I 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.

General representation of compartmental models
A typical compartmental model describes infection progression between host classes representing susceptibles, exposed, infectives and recovered individuals [11,12] 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.
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 i } has the general form: where F i represents input rate of newly infected individuals in the ith compartment and with V ± i 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 R 0 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 (M S ), exposed (M E ), infective (M I ) (SEI) model assuming that vectors do not recover from the disease. The equations for the model are as follows: Definition of R 0 in presence of seasonality is not as straightforward. The linear operator method [41,42] can be used to obtain an accurate estimate of R 0 ; however in most cases it involves numerical solutions.
Since the total population of mosquitoes is governed by dN M /dt = Π M (t) − μ M N, at equilibrium this is equivalent to a sinusoidally varying population and opens up the methods of Bacar [37] to be applied for calculating R 0 as a function of seasonality strength ε, i.e. R 0 (ε). 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 R 0 (ε) is obtained as: R 0 is the basic reproduction ratio, in the absence of seasonality, as obtained using the survival function approach [43], and Ω = 2π /T 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 R 0 (ε) < R 0 . The endemic equilibrium of the equations (2.2)-(2.4) is given by: where Q 5 = μ H2 + σ H2 , Q 6 = τ H2 + d H2 + μ H2 and For the case of small secondary host population N H2 ≈ 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 and It is reasonable to assume that a quadratic equation of the above type continues to exist even when N H2 is finite and is confirmed in numerical simulations, which we discuss below.

Complexity
Complexity of a model can be broadly divided into two categories: 1. The number of distinct host/vector species affecting the disease spread. This corresponds to distinct nodes as shown in figure 1a. 2. 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.
To  [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]. 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 R c using sensitivity analysis. The sensitivity index (γ u p ) 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]: (2.20)

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: (iv) no endemic equilibrium otherwise.
Using the condition b 0 2 − 4a 0 c 0 = 0 the critical value of R 0 (ε) for backward bifurcation, denoted by R c (ε), is then obtained as: where A and B are functions of parameters except the seasonality strength (ε), see electronic supplementary material, appendix. The range R c ≤ R ≤ 1 is the backward bifurcation region. 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 SI − M S M I approximation of the model shows that the endemic state coexists with the DFE (black horizontal line) in the backward bifurcation region R c ≤ R ≤ 1. Clearly, since R 0 < 1 is not a sufficient criterion for disease elimination, it is instructive to look at the backward bifurcation region in the Π M − μ M space, since Π M and μ M can be changed by mosquito control measures. In the Π M − μ M parameter space the area between the curves C 1 and C RC represents the backward bifurcation region as shown in figure 3a for the SI − M S M I 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 C 1 and C RC so that a larger value of Π M (for a given μ M ) 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 R c with addition of more compartments is easily seen in the R 0 − μ M space where it translates into an increase in R c 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 R c w.r.t the basic SI − M S M I model. The sensitivity indices γ R c p (see equation (2.20)) of R c w.r.t parameters (p) are tabulated in table 2. A positive (or negative) value of γ R c p indicates that an increase in parameter p also increases (or decreases) R c . According to table 2 the parameters of the model clearly form two groups based on whether they increase or decrease the value of R c , assuming an increase in these parameters. Increasing the values in the first group (b 1 , β 1 , μ H1 , ε, τ H1 ) increases while increasing the values of the second group (d H1 , μ M , σ H1 , α H1 ) decreases the threshold R c . Counterintuitively, (β H1 , Π H2 ) do not change R c , while a weakly positive but negligible dependence on σ M is observed. The parameter regimes of WNV, dengue, malaria (table 3) produce different values of R c in the R 0 − μ M space, as shown in figure 4, therefore the relative change in R c is calculated from their respective base values of R c . As an example, for WNV R c ∼ 0.609 and γ R c μ H1 ∼ 0.6415, which means that a 6.4% change in μ H1 produces a 10% change in the value of R c .
When the dead-end host is also integrated into the model, as the number of secondary hosts (given by N H2 = Π H2 /μ H2 ) is increased the backward bifurcation region decreases, as shown in figure 5   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). When seasonality is switched on (ε > 0), it is easy to see from equation (3.1) that the threshold R c increases with increasing seasonality strength. A representative bifurcation diagram shown in figure 6 illustrates this result. The bifurcation diagram was obtained by replacing Π M0 with Π M0 (1 − (ε 2 /2)G) in the nonseasonal equations. This substitution simplifies the calculation of the bifurcation diagram and is justified from the following two observations: (i) R 0 in equation (2.12) (non-seasonal case) then  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.      becomes equal to R 0 (ε) in equation (2.10) (seasonal case), and (ii) a bifurcation diagram using the linear operator method for R 0 (ε) 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 γ R c ε in table 2 confirms this result and it also shows that the WNV parameters lead to a stronger change in R c when compared to dengue/malaria for the seasonal case. 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 λ H1 . This is easily deduced from the bifurcation diagrams in figures 4-6.

Discussion
The condition R 0 < 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 SI − M S M I approximation of the model is a basic vector-host model exhibiting backward bifurcation. Then the impact of increasing complexity and seasonality on R c was studied. For instance the addition of exposed or immune compartments in hosts changes the SI − M S M I model to an SEI − M S M I or SIR − M S M I model and in both cases the threshold R c 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 SIR − M S M I model then the resulting SIRS − M S M I model has a relatively lower threshold value R c though still higher than the original SI − M S M I 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 R c w.r.t the basic SIR − M S M I 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) diseaseinduced 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 R c . 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 R 0 resulting in an increase in the value of R c ; 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 R c 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 R 0 (ε) 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.
Data accessibility. No data were used or generated during the study. Only numerical simulations were performed based on standard numerical algorithms. The parameters used in the simulation are available in the papers cited.