Leading indicators of mosquito-borne disease elimination.

Mosquito-borne diseases contribute significantly to the global disease burden. High-profile elimination campaigns are currently underway for many parasites, e.g., Plasmodium spp., the causal agent of malaria. Sustaining momentum near the end of elimination programs is often difficult to achieve and consequently quantitative tools that enable monitoring the effectiveness of elimination activities after the initial reduction of cases has occurred are needed. Documenting progress in vector-borne disease elimination is a potentially important application for the theory of critical transitions. Non-parametric approaches that are independent of model-fitting would advance infectious disease forecasting significantly. In this paper, we consider compartmental Ross-McDonald models that are slowly forced through a critical transition through gradually deployed control measures. We derive expressions for the behavior of candidate indicators, including the autocorrelation coefficient, variance, and coefficient of variation in the number of human cases during the approach to elimination. We conducted a simulation study to test the performance of each summary statistic as an early warning system of mosquito-borne disease elimination. Variance and coefficient of variation were highly predictive of elimination but autocorrelation performed poorly as an indicator in some control contexts. Our results suggest that tipping points (bifurcations) in mosquito-borne infectious disease systems may be foreshadowed by characteristic temporal patterns of disease prevalence.

Introduction

Vector-borne diseases constitute approximately 17 % of the estimated worldwide infectious disease burden (World Health Organization 2014). Mosquito-borne infections, in particular, contribute substantially to this burden, causing millions of deaths and hundreds of millions of cases annually (World Health Organization 2013). Mosquito-borne pathogen elimination remains a high priority on the global public health agenda, with elimination campaigns underway for malaria, onchocerciasis, and lymphatic filariasis (Hopkins 2013). The introduction of insecticides such as DDT in the mid twentieth century led to the implementation of global eradication campaigns for many of the most deleterious mosquito-borne diseases, including malaria, dengue virus, and yellow fever virus. These campaigns were partially successful in that diseases were eliminated in some countries and effectively controlled in others. For example, malaria was eliminated in 79 countries between 1945 and 2010 (Chiyaka et al. 2013) and the reduction of Aedes aegypti led to the near elimination of dengue and yellow fever in the Americas (Gratz 1999; Gubler 1998). However, the incidence of mosquito-borne infections has resurged since the 1970s due to a nexus of factors including globalization, increased global movement of people, development of vector resistance to insecticides, changes in land use, and loss of financial support and public health infrastructure for control and elimination (Institute of Medicine (US) Forum on Microbial Threats 2008; Gubler 2010; Mackenzie et al. 2004). Sustaining political will near the end of elimination campaigns is often difficult to achieve since there may be marginal return on investment (in terms of case reductions) (Cohen et al. 2012). Government and public health agencies require markers of progress to justify the billions of dollars that are spent on interventions (Cohen et al. 2010). Quantitative evaluation tools that will encourage governments and philanthropic organizations to choose the optimal level of investment in control and elimination activities after the initial reduction of cases has occurred are needed.

Mosquito-borne disease elimination, which involves halting transmission through interventions until no parasites remain (Breman et al. 2011; Klepac et al. 2013), e.g., through vector control strategies and/or treatment of infected individuals, corresponds to a critical transition in the disease transmission system. Criticality occurs at the point where the basic reproduction number, R0, the average number of secondary infected cases arising from a single infected case in an entirely susceptible population, is equal to one. Tipping points in complex systems may be described mathematically as bifurcations if the change in the external driver variable is slow relative to the characteristic speed of the internal variables. If control efforts are implemented sufficiently slowly over time such that elimination eventually occurs, the elimination strategy causes the system to cross a bifurcation point. The critical transition may be anticipated because prior to reaching the dynamical threshold, the system gradually loses stability (“critical slowing down” Strogatz 1994; Scheffer et al. 2009). In continuous-time systems, the loss of stability is measured by the dominant eigenvalue of the linearized system. Signatures of the slow return rate to the underlying equilibrium may be detectable through summary statistics. Statistical patterns of critical slowing down may be a diagnostic of the proximity to disease elimination (e.g., O’Regan and Drake 2013). To establish how close a disease is to elimination or emergence, typically, a stochastic epidemic model (e.g., a compartmental model or a branching process model) is fitted to time series data to estimate the basic reproduction number of a disease. The disadvantage of this approach is having to specify parametric models that make strong assumptions about the form of transmission. Additionally, the fit of stochastic epidemic models to data is usually evaluated using a goodness-of-fit criterion such as statistical likelihood, but many models may yield near-identical fits. Detection of critical slowing down would circumvent the problem of model structural uncertainty and improve efficiency by making use of higher order information in time series that is not used by other approaches that focus on trends in the mean. Methods that detect the presence of critical slowing down would therefore capture key structural features of transmission that are exhibited by a large family of mosquito-borne disease models. Non-parametric approaches for emergence and elimination prediction that are independent of model-fitting and evaluation would advance infectious disease forecasting significantly.

In this paper, we use stochastic differential equations to model the gradual implementation of four control activities on mosquito-borne disease elimination: the use of bed nets that hamper the per-capita human biting rate, e.g., (Nyarango et al. 2006), indoor residual insectide spraying procedures that shorten adult mosquito lifespan (Giardina et al. 2014) or late-acting insecticides (Read et al. 2009), the administration of drugs that reduce the human infectious period (Lawpoolsri et al. 2009), and the reduction of mosquito abundance due to administration of insecticides or elimination of breeding sites through larviciding or water drainage (Goodman et al. 1999; Tusting et al. 2013). We use the models to develop and validate leading indicator summary statistics for documenting disease elimination in mosquito-borne disease systems. We focus our analysis on malaria but our findings are relevant to other mosquito-borne infections such as yellow fever. Our results indicate that critical slowing down is detectable prior to mosquito-borne disease elimination.

Theory for elimination of vector-borne diseases

Mean field theory of Ross-Macdonald model

We consider the Ross-Macdonald model for vector-borne disease transmission (Keeling and Rohani 2008). The Ross-Macdonald model has a long history of use as a prototypical model that encapsulates the key properties of mosquito-borne pathogen transmission (Smith et al. 2012). Human hosts and vectors (mosquitoes) are assumed to be either susceptible to disease, or infectious. Denoting the numbers of infectious human host and infectious mosquito populations by H(t) and M(t) respectively, the model is The population size of human hosts N is assumed to be constant and consequently, the number of susceptible human hosts S is equal to (N−H), i.e., Similarly, mosquito abundance N is assumed constant and the number of susceptible mosquitoes S is (N−M), i.e., Parameters of the model and the expression for the basic reproduction number, R0 are listed in Table 1. The model has two equilibria, the disease-free equilibrium (0, 0) and the endemic equilibrium When R0 > 1, the endemic equilibrium is stable and the disease-free equilibrium is unstable. When R0<1, the stabilities of the equilibria are switched. Mathematically, a transcritical bifurcation, whereby the disease-free equilibrium and endemic equilibrium meet and exchange stability, occurs at the point where the basic reproduction number is equal to one (Keeling and Rohani 2008). Transcritical bifurcation diagrams corresponding to various control measures are shown in Fig. 1. Interventions due to the biting rate and mosquito abundance modify the endemic equilbrium differently than those that affect the per-capita mosquito mortality rate and per-capita human host recovery rate. The concave down trend in Figs. 1a, b suggests that successful pathogen extinction may only occur at the very end of elimination campaigns that focus on either of reduction of mosquito abundance or mosquito biting rates.

Table 1

Variables and parameters of the time-varying Ross-Macdonald models

Variable	Expression	Value
Number of infectious humans	H
Number of infectious mosquitoes	M
Per-capita human recovery rate	μ	0.01 day−1 (Smith and McKenzie 2004)
Per-capita mosquito mortality rate	δ	0.1 day−1 (Smith and McKenzie 2004)
Human population size	N h	1000
Per-capita mosquito biting rate	k	0.3 day−1 (Smith and McKenzie 2004)
Transmission efficiency from mosquitoes to humans	p	0.5 (Smith and McKenzie 2004)
Transmission efficiency from humans to mosquitoes	q	0.5 (Smith and McKenzie 2004)
Mosquito population size	N m	10000
Basic reproduction number	R_0=k^2pqN_mδμN_h
Time critical point is reached	t ∗
Critical biting rate	k^∗=(N_hδμ)/(pqN_m)
Critical mosquito population size	Nm∗=(N_hδμ)/(pqk^2)
Critical recovery rate	μ ∗=(k 2 p q N m)/(N h δ)
Critical mortality rate	δ ∗=(k 2 p q N m)/(N h μ)
Value of biting rate prior to application of control measures	k 0	0.3 day−1
Rate of change in biting rate	k 1	0.0001 day−1
Mosquito population size prior to application of control measures	N m0	10000
Rate of change of mosquito population size	N m1	1 day−1
Mosquito mortality rate prior to application of control measures	δ 0	0.1 day−1
Rate of change of mosquito mortality rate	δ 1	0.0025 day−1
Value of human recovery rate prior to application of control measures	μ 0	0.01 day−1
Rate of change of recovery rate	μ 1	0.001 day−1
Time-varying biting rate	k(t)=k_0−k_1t, t<t^∗k^∗t≥t^∗	t ∗=2800 days
Time-varying mosquito population size	N_m(t)=N_m0−N_m1t, t<t^∗Nm∗t≥t^∗	t ∗=9956 days
Time-varying mortality rate	δ(t)=δ_0+δ_1t, t<t^∗δ^∗t≥t^∗	t ∗=8960 days
Time-varying recovery rate	μ(t)=μ_0+μ_1t, t<t^∗μ^∗t≥t^∗	t ∗=2240 days
Environmental noise strength	σ	0.05

Fig. 1

Bifurcation diagrams for the Ross-Macdonald model. The stable equilibrium branches of the transcritical bifurcation ( mean field theory) as a function of each parameter affected by control activities (per-capita biting rate, mosquito population abundance, per-capita recovery rate and per-capita mortality rate) are shown. Bifurcation diagrams were plotted using the parameters given in Table 1

Control strategies for mosquito-borne diseases include bed net use, application of insecticides or modification of breeding habitat, and prompt administration of drug treatment. Control activities rolled out over long time scales relative to disease transmission cause disease prevalence to slowly decline, eventually pushing the transmission system over the tipping point to elimination (Fig. 2). For example, the rate of bed net use among children across four districts in Kenya steadily increased from 7 % in 2004–2005 to 67 % in 2006–2007 (Noor et al. 2007). To represent the implementation of control strategies over long time frames, we assume model parameters affected by control measures may be written as time-varying functions (Table 1). For example, if the per-capita biting rate k of mosquitoes decreases at rate k1 due to gradual increases in bed net usage, we write down the Ross-Macdonald equations as a fast-slow system (e.g., Kuehn 2011), where 𝜖 is a parameter that denotes the speed of evolution of the biting rate (assumed to be gradual, i.e., 0<𝜖<<1) and the function f describes the change in biting rate. The biting rate evolves over time to its critical value k∗ (achieved at R0=1 at time t = t∗) according to and consequently, the change in biting rate is slow relative to the time scale of transmission since . We use similar systems of equations to model the reduction of the duration of the infectious period (1/μ(t)), the reduction of mosquito population numbers (N(t)), and the increase of per-capita mosquito mortality rate (δ(t)) over time t through control applications (Table 1).

Fig. 2

Stochastic simulations of the Ross-Mcdonald system approaching elimination due to slow declines in a per-capita biting rate, b mosquito population size, and slow increases in c per-capita recovery rate d per-capita mosquito mortality rate. The dashed vertical line indicates the critical threshold for extinction of the pathogen in the deterministic system. The time to parasite extinction is longer in the fast-slow stochastic systems than in the corresponding deterministic systems

Stochastic description of Ross-Macdonald model

We assume population fluctuations are due to demographic and environmental stochasticity. All individuals of each type (humans, mosquitoes) have identical attributes, in that per-capita transmission and recovery rates for all human hosts are the same for all individuals and per-capita transmission and mortality rates are the same for all mosquito hosts. The human population size is assumed constant and equal to N (constant and equal to N for the mosquito population). Individuals may transition from being susceptible to being infectious. Table 2 shows the transition probability fluxes into and out of the infectious populations (H, M). We use the diffusion approximation outlined in Allen (2003) to derive stochastic differential equations that incorporate demographic stochasticity and time-varying control measures. In Appendix 1, we provide the derivation of the equations.

Table 2

Transition probability fluxes for Ross-Macdonald model. Numbers of infectious humans and populations are denoted by X=(H, M). The vector ΔX=(ΔX )=(H(t+Δt)−H(t), M(t+Δt)−M(t)) denotes change in state, i=1, 2, …, 5

Event i	Change in population sizes	(ΔX i)	Transition probability p i
Infection of susceptible humans	(H, M)→(H+1, M)	10	p_1=k(t)pM(N_h−H)N_hΔt
Recovery of infectious humans	(H, M)→(H−1, M)	−10	p 2 = μ(t)HΔt
Infection of susceptible mosquitoes	(H, M)→(H, M+1)	01	p_3=k(t)qH(N_m(t)−M)N_hΔt
Death of infectious mosquitoes	(H, M)→(H, M−1)	0−1	p 4 = δ(t)MΔt
No change	(H, M)→(H, M)	00	1−∑i=14p_i

Transitions are presented in their most general form by expressing parameters that may be influenced by control measures as functions of time

Environmental variations, e.g., in temperature, are important determinants of mosquito dynamics, driving fluctuations in mosquito abundance and influencing mosquito life histories and behaviors, e.g., (Rogers et al. 2002; Paaijmans et al. 2009; Paaijmans et al. 2013). Furthermore, there may be variability in individual adherence to drug treatment (e.g., Langton et al, 2015). We modify the equations obtained from the diffusion approximation to include the effects of environmental stochasticity in control activities (see Appendix 1 for details). We assume the stochastic process obtained from the diffusion approximation is perturbed by an additional environmental noise term that scales with the mean level of the rate impacted by a control activity at time t. Specifically, in a small time interval Δt, we assume each time-varying control measure g(t) is subject to environmental stochasticity as follows, where η is a normal random variate with mean zero and unit variance and σ denotes the strength of the environmental noise. Equation 5 replaces the appropriate parameter in the drift term obtained through the diffusion approximation (A.1 in Appendix 1), to represent an environmental perturbation to each process. For example, assuming per-capita biting rate is gradually reduced through control efforts, the change in the human host and mosquito populations in a small time interval [t, t+Δt) is where η are normally distributed random variables with mean zero and variance of unity. Letting Δt→0, and assuming existence and uniqueness of the stochastic integral (Allen 2003), then and the system of equations converges in the mean square sense to a system of Ito stochastic differential equations, The systems of stochastic differential equations in Table 3 corresponding to each control activity were derived in the same way (Appendix 1). Simulations of each model approaching elimination are illustrated in Fig. 2.

Table 3

Time-varying Ross-Macdonald equations with demographic and environmental stochasticity

Biting rate	dH=(k(t)pN_hM(N_h−H)−μH)dt+k(t)pN_hM(N_h−H)+μHdW_1+σMp(N_h−H)N_hdW_3
	dM=(k(t)qN_hH(N_m−M)−δM)dt+k(t)qN_hH(N_m−M)+δMdW_2+σHq(N_m−M)N_hdW_3
Mosquito population size	dH=(kpN_hM(N_h−H)−μH)dt+kpN_hM(N_h−H)+μHdW_1
	dM=(kqN_hH(N_m(t)−M)−δM)dt+kqN_hH(N_m−M)+δMdW_2+σkHqN_hdW_3
Recovery rate	dH=(kpN_hM(N_h−H)−μ(t)H)dt+kpN_hM(N_h−H)+μ(t)HdW_1+σHdW_3
	dM=(kqN_hH(N_m−M)−δM)dt+kqN_hH(N_m−M)+δMdW_2
Mortality rate	dH=(kpN_hM(N_h−H)−μH)dt+kpN_hM(N_h−H)+μHdW_1
	dM=(kqN_hH(N_m−M)−δ(t)M)dt+kqN_hH(N_m−M)+δ(t)MdW_2+σMdW_3

Five hundred simulations of each set of equations were performed. Terms under square roots represent the G 11 and G 22 entries in the diffusion matrix G(t)(A.5), whereas terms that scale with the environmental noise strength σ are the G entries

Leading indicators of elimination

To calculate leading indicators of disease elimination, we need to analyze the properties of the time-varying Ross-Macdonald system, in the neighborhood of the equilibrium point that would be present in the 𝜖→0 limit (e.g., O’Regan and Drake 2013; Kuehn 2011). In the stochastic system, the equilibrium is a quasi-stationary state. Quantifying the behavior of the deviations from the equilibrium point can be achieved by deriving a locally linear approximate description of the probability distribution that is the solution of the forward Kolmogorov equation (A.6). Since the Ross-Macdonald system does not exhibit long transients and the time-varying functions represent non-stationary processes that are changing gradually through time, linearization of the drift and diffusion terms of the equations in Table 3 gives reasonable information for the qualitative behavior of fluctuations about the quasi-stationary state.

Solutions of the linearized forward Kolmogorov equation (A.7) have the same probability distribution as the system of the stochastic differential equations (A.8). In Appendix 1, we apply the Fourier transform to these equations to derive the power spectrum of the fluctuations in human cases. Finally, the leading indicator statistics (variance, autocorrelation coefficient, and the coefficient of variation) are obtained from integration of the power spectrum. The expressions for the leading indicators are found in Table 4 and are expressed in terms of the eigenvalues of the Jacobian matrix of the system linearized at the endemic equilibrium. In Appendix 1, we show that the endemic equilibrium of system (1) is always a stable node if the parameters of the model are positive. Consequently, we present the expressions for the statistics in terms of the real and distinct dominant and subdominant eigenvalues, λ1 and λ2, respectively. Changes in these summary statistics as the bifurcation point is approached are indicators of critical slowing down.

Table 4

Analytical expressions for quasi-stationary statistics about the endemic infectious human quasi-steady state H ∗ expressed in terms of the eigenvalues

Power spectrum S I(ω)	2(α_H+D_11ω^2)(ω^2−λ_1λ_2)^2+(λ_1+λ_2)^2ω^2
Autocorrelation	1π∫0∞S_I(ω)cos(ωτ)dω
Variance σ 2	α_H+λ_1λ_2D_112λ_1λ_2(−λ_1−λ_2)
Coefficient of variation	(α_H+λ_1λ_2D_11)^122λ_1λ_2(−λ_1−λ_2)H^∗

The equilibrium of the Ross-Macdonald model, (H ∗, M ∗), is a stable node (Appendix 1) and thus the Jacobian matrix has two real, negative, distinct eigenvalues, λ 1 and λ 2. Variables for each model are described in Tables 1,  5, and 6. The expressions for the power spectrum are multiplied by 2 because they are evaluated over the frequency domain . No closed-form expression for the lag- τ autocorrelation is known and so it must be evaluated numerically

Table 5

Variables substitutions for the stochastic differential equations for the fluctuations about the endemic equilibrium (2) at time t (A.12)

Variable	Expression
H ∗	N_h(k(t)^2N_m(t)pq−N_hδ(t)μ(t))k(t)q(k(t)N_m(t)p+N_hμ(t))
M ∗	k(t)^2N_m(t)pq−N_hδ(t)μ(t)k(t)p(k(t)q+δ(t))
a 11	−k(t)p M ∗/N h−μ(t)
a 12	k(t)p(N h−H ∗)/N h
a 21	k(t)q(N m(t)−M ∗)/N h
a 22	−k(t)q H ∗/N h−δ(t)
d	a 11 a 22−a 12 a 21
T	−a 11−a 22
α H	a222D_11+a122D_22−2D_12a_12a_22

If the statistics were evaluated about the trajectory of the fast-slow system, then H ∗ and M ∗ in the table below are replaced with H(t) and M(t) respectively (e.g., Fig. 4). The expressions for the D coefficients are in Table 6. For fluctuations about an endemic equilibrium, the value of control parameters k(t), N (t), μ(t) and δ(t) are constant

Table 6

Variables substitutions for terms of the variance-covariance matrix of time-varying Ross-Macdonald models with demographic and environmental noise

Coefficient	Biting rate	Recovery rate
D 11	k(t)pN_hM^∗(N_h−H^∗)+μ(t)H^∗+(σ/Nh)^2p^2(M^∗)^2(N_h−H^∗)^2	kpN_hM^∗(N_h−H^∗)+μ(t)H^∗+(σH^∗)^2
D 12	(σ/N h)2 p q M ∗ H ∗(N h−H ∗)(N m−M ∗)	0
D 21	(σ/N h)2 p q M ∗ H ∗(N h−H ∗)(N m−M ∗)	0
D 22	k(t)qN_hH^∗(N_m−M^∗)+δM^∗+(σ/Nh)^2q^2(H^∗)^2(N_m−M^∗)^2	kqN_hH^∗(N_m−M^∗)+δM^∗
Coefficient	Mosquito population size	Mortality rate
D 11	kpN_hM^∗(N_h−H^∗)+μH^∗	kpN_hM^∗(N_h−H^∗)+μH^∗
D 12	0	0
D 21	0	0
D 22	kqN_hH^∗(N_m(t)−M^∗)+δM^∗+σ^2kqN_h(H^∗)^2	kqN_hH^∗(N_m−M^∗)+δ(t)M^∗+(σM^∗)^2

If the statistics were evaluated about the trajectory of the fast-slow system, then H ∗ and M ∗ in the table below are replaced with H(t) and M(t), respectively (e.g., Fig. 4). For fluctuations about an endemic equilibrium, the value of control parameters k(t), N (t), μ(t) and δ(t) are constant

Fig. 4

To obtain predictions for how the summary statistics behave as elimination is approached, mean leading indicators were calculated numerically using parameter values relevant for malaria (Table 1). The vertical dashed line in each figure indicates the threshold per-capita recovery rate and threshold per-capita mosquito mortality rate at R 0=1 respectively. a Lag-1 autocorrelation and coefficient of variation are predicted to increase as control measures that affect the human infectious period are applied but variance is predicted to decrease. b Lag-1 autocorrelation, variance and coefficient of variation are predicted to increase as the per-capita mosquito mortality rate increases due to control activities. Here, we compare the statistics evaluated at the equilibrium (H ∗, M ∗) and along the fast-slow trajectory (H(t), M(t)). We note that the variance is non-monotonic if evaluated along the trajectory, but there is agreement further from the threshold

Simulations

The preceding sections present an analytical theory of leading indicators of elimination for stochastic time-varying Ross-Macdonald models. To investigate the results of this theory for a particular parameter set, we calculated leading indicators of disease elimination in infectious human hosts using the models in Table 3, assuming that (a) the mean number of infectious individuals in the 𝜖→0 limit is given by the deterministic endemic equilibrium (H∗, M∗) (2) or (b) by assuming it is given by the current state (H(t), M(t) of the system approaching elimination. We selected parameters consistent with malaria (Smith and McKenzie 2004) (Table 1).

To test the robustness of the theoretical predictions in Table 4, we simulated the approach to elimination 500 times for each of the models in Table 3. For each control scenario, the infectious time series approaching elimination were sampled at weekly intervals in each simulation, assuming perfect detection (no underreporting). The transcritical bifurcation in these scenarios was approached over a long time frame (e.g., 2800 days (400 weeks) when k1=1/10, 000 day−1) in the fast-slow Ross-Macdonald model approaching elimination. The length of the time series in each simulation set depended on how long it took to reach the critical threshold R0=1 in each model (400 weeks in the changing biting rate model, 1280 weeks in the changing mosquito mortality model, 320 weeks in the changing recovery rate model, and 1424 weeks in the changing mosquito abundance model, respectively).

Analysis over a moving window

Our simulation study evaluates the performance of summary statistics (lag-1 autocorrelation, variance, and coefficient of variation) as an early warning system. Here, we describe an algorithm to process input time series data that is independent of a specific model. To investigate the robustness of the changes in the early warning theoretical predictions over a moving window, i.e., as they would be used in online analysis of surveillance data, we used Gaussian filtering to remove the influence of the slowly varying trend (Dakos et al. 2008). We fitted a Gaussian kernel smoothing function with a fixed bandwidth across the infectious human host time series up to the time that the transcritical bifurcation was predicted (t∗). We obtained the residuals by subtracting the fit from each time series. We calculated the lag-1 autocorrelation, variance, and coefficient of variation of the residuals over a moving window half the length of each time series. We calculated the lag-1 autocorrelation coefficient of each replicate using the acf function in R. The coefficient of variation was found by calculating the mean and standard deviation of each infectious replicate. The median and 95 % prediction intervals for each of the statistics were calculated over the 500 replicates of each model. The prediction intervals were calculated using the quantile function in R. To quantify the association between time and the statistic for each replicate, we used Kendall’s correlation coefficient τ, a non-parametric statistic of association between two quantities that has values between 1 (positive correlation) and -1 (negative correlation). By repeating the calculation for each realization of leading indicators, we generated distributions of the temporal correlation τ for all of the simulation sets approaching criticality.

To assess leading indicator performance, we adapted the method described in Boettiger and Hastings (2012) to distinguish between statistics obtained from systems approaching elimination from those calculated from quasi-stationary systems. Distributions of Kendall’s τ were additionally generated using the statistics obtained from realizations of quasi-stationary epidemic systems (null models) that were processed using the procedure outlined above. That is, the null models have the same environmental noise structure as the test models (the systems of stochastic differential equations in Table 3) but assume k(t) = k0, N(t) = N, δ(t) = δ0 and μ(t) = μ0 in each model respectively. Null models were initialized from the deterministic equilibrium calculated from the values of k0, N, μ0, and δ0 in Table 1 and were simulated for the same length of time as required for the transition to be approached in the test models (systems approaching elimination).

From the distributions of Kendall’s τ obtained from null and test model realizations, we calculated receiver-operating characteristic (ROC) curves. The receiver-operating curve plots diagnostic sensitivity (true positive rate, y-axis) as a function of one minus diagnostic specificity (false positive rate). The false positive rate is the integral of the distribution of the correlation coefficent under the quasi-stationary system, to the right of a threshold line, and the true positive rate is the integral of the distribution of the correlation coefficent under the system approaching elimination. These rates make up the curve. If curves lie on or near a line with a constant slope of unity, this would imply the distributions overlap completely, and they cannot be used to distinguish between quasi-stationary systems from those approaching criticality (Boettiger and Hastings 2012; O’Regan and Drake 2013). The magnitude of departure from this line is summarized by the area under the curve, which reaches its maximum at one. If the area under the curve is close to one, then there is near perfect detection of sensitivity.

Underreporting

To assess the effects of underreporting on leading indicator performance, we binomially sampled each time series in each set of stationary and supercritical intervention simulations with rates of detection ranging from 20 to 90 %. We calculated the ROC curves and reported the area under the curve for each simulation set.

Results

Predictions for mean behavior of leading indicators

Leading indicators of mosquito-borne parasite elimination exhibit systematic changes as the critical point is approached. Figure 3 shows the theoretical predictions for the lag-1 autocorrelation, variance, and coefficient of variation when the per-capita biting rate and mosquito abundance respectively gradually decline over time. In both cases, the lag-1 autocorrelation, variance, and coefficient of variation are predicted to increase as the system nears criticality. However, the increase in variance is not sustained as the biting rate of mosquitoes is lowered. Figure 4 shows the predictions as the per-capita recovery and mosquito mortality rates respectively are increased. The lag-1 autocorrelation and the coefficient of variation are both predicted to increase as control measures are applied but the variance is predicted to decline as the infectious period shortens (Fig. 4a). Variance is predicted to increase as the per-capita mosquito mortality rate increases but the variance is a non-monotonic function of mortality rate (Fig. 4b). In all cases in Figs. 3 and 4, the theoretical predictions were calculated using the eigenvalues arising from linearizing each system about its equilibrium because the equilibrium was sufficiently close to the trajectory of the fast-slow systems in Table 1. Only for the mosquito mortality rate was there a difference between the theoretical equilibrium and the observed trajectory (Fig. 4b).

Fig. 3

To obtain predictions for how the summary statistics behave as elimination is approached, mean leading indicators were calculated numerically using parameter values relevant for malaria (Table 1). The vertical dashed line in each figure indicates the threshold mosquito population abundance and threshold per-capita biting rate at R 0=1 respectively. a Lag-1 autocorrelation and coefficient of variation are predicted to increase as control measures impacting per-capita biting rate are applied but variance becomes non-monotonic close to the critical point. b Lag-1 autocorrelation, variance, and coefficient of variation are predicted to increase as control actitivies affecting mosquito population abundance are applied

Simulation study results

The predictions for the statistics obtained from the simulations are broadly consistent with the theoretical predictions for the mean statistics. In Figs. 5 and 6, the trends in the median and 95 % prediction intervals for the summary statistics (autocorrelation, variance, coefficient of variation) as per-capita biting rate and mosquito abundance respectively agree with the mean theoretical predictions (Fig. 3). Additionally, the areas under the ROC curves are close to 1, indicating that it is possible to distinguish between quasi-stationary systems and systems approaching elimination due to control activities that impact mosquito behavior and abundance. However, predictions from simulations over a moving window obtained from reducing the infectious period of human hosts and the lifespan of mosquitoes respectively are not as robust (Figs. 7 and 8). As the infectious period is shortened, theoretical predictions for variance and coefficient of variation are robust over a moving window, but the AUC value for lag-1 autocorrelation indicates that autocorrelation is a less accurate indicator of elimination in this scenario. The less accurate performance of autocorrelation as an indicator is also seen as mosquito lifespan is shortened due to control measures. Variance is predicted to decline as elimination is approached, matching the theoretical prediction along the trajectory (Fig. 4b).

Fig. 5

Simulation study results arising from arising from reduction in per-capita biting rate. Note that the value of the biting rate is continuously changing over the 200-week window, at a rate of . Panels a, c, and e show the median statistics (thick lines) and 95 % prediction intervals (shaded regions). The dashed vertical line marks the time of the transcritical bifurcation. The trends in the median statistics agree with the mean theoretical predictions (Fig. 3a). Panels b, d, and f show the results of the ROC analysis. The AUCs are high, indicating it is possible to distinguish between the stationary system and one slowly approaching elimination. A bandwidth of 80 weeks was selected for Gaussian filtering

Fig. 6

Simulation study results arising from reduction in mosquito population size. Note that the value of mosquito abundance is continuously changing over the 712-week window, at a rate of . Panels a, c, and e show the median statistics (thick lines) and 95 % prediction intervals (shaded regions). The dashed vertical line marks the time of the transcritical bifurcation. Theoretical predictions for the trends in each summary statistic are robust over a moving window. The AUCs are high, indicating it is possible to distinguish between the stationary system and one slowly approaching elimination. A bandwidth of 80 weeks was selected for Gaussian filtering

Fig. 7

Simulation study results obtained from reduction in human infectious period. Note that the value of the recovery rate is continuously changing over the 160-week window, at a rate of . Panels a, c, and e show the median statistics (thick lines) and 95 % prediction intervals (shaded regions). The dashed vertical line marks the time of the transcritical bifurcation. Theoretical predictions for trends in variance and coefficient of variation are robust over a moving window. Panels b, c, and f show the performance of the statistics, assessed through ROC analysis. The AUCs are high, indicating it is possible to distinguish between the stationary system and one slowly approaching elimination but the AUC value for the autocorrelation indicates it is less accurate in distinguishing between the stable system and one approaching criticality. A bandwidth of 80 weeks was selected for Gaussian filtering

Fig. 8

Simulation study results obtained from reduction in per-capita mosquito mortality rate. Note that the value of the mortality rate is continuously changing over the 640-week window, at a rate of . Panels a, c, and e show the median statistics (thick lines) and 95 % prediction intervals (shaded regions). The dashed vertical line marks the time of the transcritical bifurcation. Theoretical predictions evaluated about the trajectory (red lines in Fig. 4b) are robust over a moving window. Panels b, c, and f show the performance of the statistics, assessed through ROC analysis. The AUCs are high, indicating it is possible to distinguish between the stationary system and one slowly approaching elimination. A bandwidth of 80 weeks was selected for Gaussian filtering

Imperfect detection results

As expected, the accuracy of the leading indicators, measured by AUC, generally increases with detection rate (Fig. 9). The coefficient of variation appears to be the most robust indicator under the effects of observation error while the performance of variance and autocorrelation generally appear to be sensitive to the effects of underreporting.

Fig. 9

Effects of imperfect detection on leading indicator performance for each intervention: a Changing mosquito abundance, b changing biting rate, c changing recovery rate, and d changing mosquito mortality rate. Each time series was binomially sampled with detection rates ranging from 20 to 90 %. The area under the curve (AUC) is graphed as a function of detection rate. The coefficient of variation (CV) appears to be the most robust indicator to underreporting

Discussion

Detecting the onset of critical transitions in infectious disease systems such as the emergence of novel pathogens and the elimination of disease would be of tremendous value for public health. The goal of our study was to develop a theory of leading indicators for vector-borne infectious disease transmission systems that follow the assumptions of the Ross-Macdonald model and that are forced through a critical transition through gradual increases in control efforts. Our main results include analytical expressions for summary statistics for a family of time-varying Ross-Macdonald models approaching elimination (Table 4). Numerical calculations assuming time-varying control activities indicate that trends in the observable statistics are discernible. Testing the robustness of these predictions via a simulation study suggests that critical slowing down is detectable in human-host fluctuations of Ross-Macdonald type mosquito-borne disease systems. These results are in broad agreement with the findings that variance, autocorrelation, and coefficient of variation were predictive of disease elimination, in supercritical SIR and SIS compartmental models with time-varying vaccination uptake and transmission rates respectively (O’Regan and Drake 2013).

The simulation study suggests that variance and coefficient of variation perform well as leading indicators of elimination in all control contexts, but autocorrelation performs poorly if control strategies affect the removal rates of either human host or mosquito populations, manifested as reduction in human infectious period or mosquito lifespan respectively. In contrast, if control activities influence the transmission rate, either by reducing mosquito abundance or per-capita biting rate, all of the statistics, including autocorrelation, perform strongly. The poor performance of autocorrelation is somewhat surprising, as it has been suggested that autocorrelation may be the most robust indicator of critical transitions, at least for fold bifurcations (Dakos et al. 2012). One hypothesis for this finding is how the eigenvalues of the system change as criticality is approached. For the numerical parameters considered here, time-varying control measures that result in reduction of transmission result in eigenvalues that decline with control (albeit at different rates), whereas measures that reduce the removal rate can lead to non-monotonic eigenvalues, e.g., the subdominant eigenvalue exhibiting increases in magnitude, rather than decreasing in magnitude, as might be expected (Fig. 10). More research is needed to investigate how dynamical changes in eigenvalues impact the summary statistic predictions generally. Interestingly, the autocorrelation appears to respond earlier than the coefficient of variation to interventions (Figs. 5,  6,  7, and 8). The speed of the response of leading indicators to gradual parameter changes, and determining the factors most important in driving statistical patterns of critical slowing down, are subjects of ongoing research.

Fig. 10

Eigenvalues obtained from linearization about the stable node equilibrium corresponding to each control activity. The red dashed vertical line corresponds to the critical value of each parameter where R 0=1. Parameter values relevant for malaria were used (Table 1)

The ROC procedure we describe here is not appropriate for measuring leading indicator statistics in time series of case reports. What is observed in practice is a single realization (case report time series) that may lie far from the mean model prediction due to stochastic fluctuations. Rather, we use the theory to calculate the theoretical values of the mean leading indicator summary statistics and compare changes in them as the critical point is approached with changes in median statistics from 500 stochastic simulations over a moving window (“Simulations” and “Analysis over a moving window” sections; and “Results” section). Changes in summary statistics as the critical point is approached are signatures of critical slowing down. Our simulation study tests the performance of each summary statistic as an early warning system, which we define as an information processing system that returns a binary signal at each time in the simulation.

Our online processing algorithm assumes that the input data is the number of infectious individuals at time t, but in practice, data would consist of incidence time series (the number of new cases at time t). Incidence data will closely agree with the number of currently infected individuals if the incidence time series has the same resolution as the duration of the infectious period (Ferrari et al. 2005). If the system is quasi-stationary or moving towards elimination sufficiently slowly, then the number of infectious individuals at a time equal to the infectious period should equal the number of new cases. However, this issue becomes complicated for diseases such as malaria where case reports in elimination settings may only include very sick individuals that seek treatment and not those who are infected but are not sick enough to seek treatment. Moreover, our results in Fig. 9 suggest that the effects of imperfect detection are not uniform across interventions. Development of theory for detection of critical slowing down accounting for underreporting and the presence of unobserved processes is needed.

The Ross-MacDonald model has long been used as a strategic model for understanding mosquito-borne infectious disease (Smith et al. 2012), but has not been studied under the realistic scenario of a slow environmental forcing. Actual elimination requires the growth of institutions (like hospitals, central government, land conversion like ditching, and mosquito control networks) that are slow to develop and need to be maintained until elimination is reached. While many control campaigns could be considered “pulse” interventions, in that the intervention is deployed as quickly as possible, the models considered here encompass “press” strategies (interventions that are maintained until elimination is reached). Our approach ignores integrated control measures that impact multiple components of vectorial capacity simultaneously, e.g., insecticides targeted on adult mosquitoes that affect mosquito abundance and adult longevity. However, it would be straightforward to derive the theoretical predictions for any such model, simply by appropriate modification of the formulas for the eigenvalues. A limitation of our time-varying models is that they do not include the extrinsic incubation period (EIP), to which transmission can be sensitive to variation. Including a delay in the form of an exposed mosquito class could be important for anticipating elimination if evolution-proof control activities that affect mosquito lifespan are implemented since sustained transmission between hosts and vectors requires the EIP to be less than adult lifespan. Additionally, our simple model of mosquito-borne disease elimination does not include many realistic aspects of mosquito-borne dynamics including waiting time distributions (Smith and McKenzie 2004), temperature dependence in the gonotrophic cycle and mosquito development time (Paaijmans et al. 2013; Read et al. 2009), spatial location of hosts and vectors, and spatial heterogeneity in adult and larval mosquito habitats (Hollingsworth et al. 2014; Perkins et al. 2013; Reiner et al. 2013). Future tactical models examining elimination dynamics should include these complexities but since most models of vector-borne diseases are of the Ross-Macdonald type (Reiner et al. 2013), it is reasonable to develop theoretical predictions for summary statistics using the existing Ross-Macdonald theory, combined with simple assumptions of control activities applied to malaria.

In conclusion, to our knowledge, this study constitutes the first theory for non-parametrically anticipating elimination of mosquito-borne diseases. Our analysis of Ross-Macdonald models parameterized for malaria shows that critical slowing down is detectable in fluctuations in human cases for mosquito-borne parasites approaching elimination, suggesting that disease elimination may be anticipated even in the absence of a detailed understanding of underlying mechanisms. Our work suggests that online algorithms for detecting changes in leading indicators may be achievable, possibly aiding sustainment of the gains made by elimination programs.