Assessing the suitability for Aedes albopictus and dengue transmission risk in China with a delay differential equation model

Dengue is considered non-endemic to mainland China. However, travellers frequently import the virus from overseas and local mosquito species can then spread the disease in the population. As a consequence, mainland China still experiences large dengue outbreaks. Temperature plays a key role in these outbreaks: it affects the development and survival of the vector and the replication rate of the virus. To better understand its implication in the transmission risk of dengue, we developed a delay differential equation model that explicitly simulates temperature-dependent development periods and tested it with collected field data for the Asian tiger mosquito, Aedes albopictus. The model predicts mosquito occurrence locations with a high accuracy (Cohen’s κ of 0.78) and realistically replicates mosquito population dynamics. Analysing the infection dynamics during the 2014 dengue outbreak that occurred in Guangzhou showed that the outbreak could have lasted for another four weeks if mosquito control interventions had not been undertaken. Finally, we analyse the dengue transmission risk in mainland China. We find that southern China, including Guangzhou, can have more than seven months of dengue transmission per year while even Beijing, in the temperate north, can have dengue transmission during hot summer months. The results demonstrate the importance of using detailed vector and infection ecology, especially when vector-borne disease transmission risk is modelled over a broad range of climatic zones.


S1.1 Deriving DDEs from ODEs
The terms contributing to the rate of change for each class can be classified into three categories: recruitment, maturation, and mortality. The transformation is demonstrated with the help of class E, eggs.
We assume there are E 0 eggs at t = t 0 and, for simplicity, that no eggs are laid after t = t 0 . Eggs experience a temperature-dependent mortality of µ(T (t)) = µ(t): This equation can easily be integrated to After an egg development period of τ E , all eggs mature to the juvenile stage, thus at t = τ E : with δ τ E = 1/∆t at τ E and 0 otherwise. Finally, we include egg recruitment. Class E is recruiting from mature females A with a rate of β (1 − ω) (egg laying rate times a factor for diapausing eggs). and Combining recruitment, maturation, and mortality gives The same applies to immature female development with a development period of τ I , and transition of exposed to infectious mosquitoes after the EIP.
Due to larval density dependence and the additional mortality term, juvenile mortality does not follow an exponential decay. Integration of is achieved by substitution of x(t) := J(t) exp( t 0 µ(σ) dσ). Then Integration yields so that gives

S1.2 Time delay calculation
Before we run any simulations, we first calculate development times for immature life stages (τ Et , τ Jt , τ It ) and the extrinsic incubation period (EIP t ) for the temperatures given at all days of the year t = 0 . . . 365.
Starting from any given time point t 0 , we sum over all the following inverse development times 1 denote the fraction of daily development, until this sum reaches 1: For example, if for days t = 1, 2, 3 temperatures would drop slightly and the juvenile development times would be τ J = 10, 12, 13 days. We would calculate that after these three days, 1 10 + 1 12 + 1 13 = 26% of the development would be completed. We would continue to add inverse development times until we reach 100% development and take the final number of days as development period.
Note that this approach is independent of the step size chosen to solve the numeric equations. The approach by (author?) [1], in which the time delay itself is modelled by a differential equation, is dependent on the step size. If the step size is not chosen small enough, this could for example lead to unexpected differences in an incubation period before and after a spike in temperature.

S1.3 ODE vs. DDE dynamics
The extrinsic incubation period (EIP) is the period a mosquito spends infected but not yet infectious. During this time, the pathogen (the dengue virus in our case) has to replicate and spread in the mosquito body before its saliva gets infectious and the pathogen can be transmitted to a susceptible person during a bite.
There are different ways to model the transition from infected to infectious mosquitoes, e.g. with ordinary differential equations (ODEs) or delayed differential equations (DDEs). Using ODEs or DDEs imposes different assumptions on the infection dynamics.
For this comparison, we will assume that we only have the host and no vector, i.e. infected hosts can directly infect susceptible hosts like influenza for example. The same principles are applicable to vector-host models though.

ODE model
SEI denotes susceptible, exposed, and infectious mosquitoes for the infection dynamics. In ODEs, the rate of change of the classes u = (S, E, I) depends on a reaction function f that describes a linear or non-linear relationship of these classes; d dt u(t) = f (u(t)). The (simplified) equations are given by: with β the infection probability, µ the mortality rate that is the same for all three classes, and γ the incubation rate, i.e. the inverse of the EIP, γ = 1 EIP .

DDE model
Another approach to model this incubation period are delayed differential equations (DDEs). Here, the rate of change at a certain time not only depends on the current state but also on previous states at time t − τ ; In this case, the time delay τ is the EIP. The transition of mosquitoes into class I does not start before the EIP is over. But each time step, a certain proportion of mosquitoes die. After the EIP has passed, all the mosquitoes that have not died in the meantime move to class I. The factor e − t t−τ µdσ denotes for this mortality experienced over the EIP. Note that the incidence rate βSI is written without the delay terms in the susceptible's equation. This is necessary as otherwise infected but not yet infectious individuals remained in the susceptible class for the period τ and could then get infected again. The total number of individuals are then given by N := S + I + E.
The DDE model could only be calculated with classes S and I but in our model, we need to calculate class E as the exposed females also lay eggs.

Comparison
The difference between the two model assumptions becomes clear when comparing single infection cases directly, see Figure A. After a single infection case is introduced from outside at day 0, e.g. through a mosquito bite, the first infectious individual in the ODE model appears at day 1 while it appears only after EIP = 10 days with the DDE model.
Using ODEs assumes that after the end of the EIP, 50% of exposed mosquitoes have become infectious and changed into class I. However: before and after the end of the EIP, there are constantly mosquitoes changing from class E to class I. The distribution of transitions is not bell shaped (the transitions are not normally distributed) but rather resembles an exponential decay which is problematic: Most mosquitoes change from class E to I after day 1, regardless of the EIP, and a tiny fraction will still wait to become infectious in 100 years. In contrast, using DDEs assumes that after the end of the EIP, the mosquitoes that do not die in the exposed state, become infectious all at once. If τ is small, the models give similar results which is in accordance with (author?) [2]. The more the EIP or τ increases, the more the output of the ODE and DDE models will differ. Figure B shows the disease dynamics without mortalities for τ = 4 and τ = 18 days.
It is thus important to use DDE when analysing disease dynamics in which the EIP can take more than a week due to its temperature-dependence. Otherwise, the length and amplitude of disease outbreaks could be overestimated.

S1.4 Sensitivity Analysis
To investigate the influence of each parameter on the cumulative number of dengue cases in 2014, we perform the elementary effects test, EET [3], see Figure C. The EET measures the influence of single input parameters on model outputs, as well as their degree of interaction with other parameters. Latin hypercube sampling is used to vary parameters in the range of ±10% of the standard setting [4]. The model is then run with Guangzhou 2014 climate data until convergence is reached. The total number of cases at the end of the year is then retained as reference. Octave scripts for this method are derived from the SAFE toolbox [5]. The total number of dengue cases will be strongly associated with the basic reproduction number for vectorborne diseases: This relationship is reflected by our sensitivity analysis: The biting rate, α, and the mosquito's adult mortality, µ A , that appear twice in the formula, have the biggest effect on the simulated number of cumulative case, followed by the extrinsic incubation period, EIP, and transmission probabilities from host to vector, b HV , and vector to host, b V H . Factors that have a larger impact on the vector-to-host ratio m, including the environmental carrying capacity K, and development periods for juvenile and immature females, τ J and τ I , also have larger impacts. The other remaining factors have smaller impacts on m and the cumulative number of dengue cases. The distributions for mean and standard deviation of EEs indicate that those parameters that interact with other parameters also have a bigger effect on the model output.

S1.6 Equilibria
In order to determine model equilibria, we first write down our full equation system: We now assume that temperature and precipitation are constant and thus that all parameters are stable At equilibria DDEs behave like ODEs, as in the limit x(t − τ ) = x(t). Thus we get: Assuming E, J, I, A, E d all positive, we get a trivial, x ⋆ 1 = 0, and a non-trivial solution, and I ⋆ and J ⋆ being longer terms we do not write down.
It is now of interested for what parameter range either equilibrium is stable. If the trivial equilibrium is stable, the mosquito population will not be able to successfully reproduce in the long run.

Stability
We now analyse the stability of the trivial equilibrium. Let x ⋆ be an equilibrium of equation and let δx(t) be the displacement from this equilibrium, assumed small. We assume f is smooth around x ⋆ and thus has derivatives of all orders. Accordingly, Linearisation, using a Taylor series, yields with J 0 being the Jacobian with respect to x, and J τ i being the Jacobians with respect to x(t−τ i ), evaluated at x(t) = x(t − τ 1 ) = . . . x(t − τ n ) = x ⋆ . Following (author?) [7], we now assume that equation (2) has an exponential solutions, so that we can write δx(t) = Be λt , and thus we get Calculating the Jacobian matrices at the zero equilibrium x ⋆ = 0 gives for the normal Jacobian with respect to x(t) and Jacobians J τ E , J τ J , J τ E +τ J ,. . . with respect to x(t − τ E ), As all our DDEs in the fully written equation system only have delays for A and E d , the remaining Jacobians For our characteristic equation, we get |J 0 + e −λ(τ E +τ J +τ I ) J τ E +τ J +τ I + e −λ(τ J +τ I ) J τ J +τ I − λI| = 0, and calculating the determinant yields: Substitution of λ := u + iv and using Euler's formula gives two equations that show that λ can be real or complex, with either positive or negative real parts. We thus try to reduce the complexity and only look at the active mosquito season with ω = 0 and φ = 0 to calculate λ for our standard parameter set. We get: and thus with W (x) denoting the Lambert-W function. Calculating λ with actual values for τ E (T ), µ A (T ) etc., shows if the trivial equilibrium is stable or unstable. In this case, the population reaches a non-trivial equilibrium for the current model setting when the constant temperature is higher than 16.5°C and lower than 32.5°C.
Thus, the temperature of the active mosquito season (roughly April until end of September) should be between these two values for a successful establishment of Ae. albopictus.

S1.9 DDE vs. ODE transmission risk
We repeated the simulations to calculate the potential length of the dengue transmission season using the equivalent ODE model (compare section S1.1) with identical parameters, see Figure H.