Dengue transmission: mathematical model with discrete time delays and estimation of the reproduction number

ABSTRACT In this paper, we establish a mathematical model with two delays to reflect the intrinsic and extrinsic incubation periods of virus in dengue transmission. The basic reproduction number of the model is defined. It is proved that the disease-free equilibrium is stable when and the positive equilibrium is stable when . Next, we derive an estimation formula for the reproduction number when the human population is partially susceptible to dengue. As an application, the values of dengue transmission in Singapore in the years 2013–2015 are estimated. Our estimation method can be applied to estimating of other infectious diseases, especially when the human population is not completely susceptible to the disease.


Introduction
The incidence of dengue has grown dramatically around the world in recent years. The disease is now endemic in more than 100 countries especially in tropical and subtropical regions. It is estimated that over 40% of world's population are at risk from dengue and the dengue infections worldwide is 50-100 millions every year [56]. Dengue is a vectorborne disease transmitted by the vector Aedes aegypti, a type of mosquito inhabited in urban and suburban area. The virus is transmitted to human via the bites of infected female mosquitoes. After virus incubation period of 1-16 days (extrinsic incubation time [23]), an infected mosquito can transmit the virus for the rest of its life. Upon bitten by an infected mosquito, an infected person will show symptoms after an incubation period of 1-11 days (intrinsic incubation time [23]) and the symptoms usually last for 2-7 days [26,41,46]. One can see that both the intrinsic incubation period and the extrinsic incubation period of the virus are considerably long and hence significant.
Mathematical models have been widely used to quantitatively analyse the transmissions of various infectious diseases [1,[3][4][5][6][7]10,14,17,20,29,32,48,50]. These models are mainly compartmental models where a population is divided into different compartments corresponding to the different courses of a disease. Among models for dengue transmission, the human population is usually divided into three compartments: the susceptible (S H ), the infectious (I H ) and the recovered (R H ); while the vector population is divided into two compartments: the susceptible (S V ) and the infectious (I V ) [2,9,13,18,21,25,44].
We know that considerable incubation time of dengue virus, both intrinsic and extrinsic, is present in dengue transmission. Indeed, to model the incubation time of pathogens of various infectious diseases, such as malaria, HIV/AIDS, tick-borne diseases, etc., time delays have been widely used [32,33,38,42,48,51,52,57]. However, in the literature of dengue transmission, there are few models with time delay(s) established and systematically analysed [31]. Motivated by this, the first aim of this paper is to establish a compartmental model with two discrete delays to reflect the intrinsic and extrinsic incubation times present in dengue transmission.
We note that a term, such as S(t − τ )I(t − τ ), is usually adopted to incorporate the incubation time τ [23,32,47,51]. Meanwhile, the following delayed SIR model with nonlinear incidence rate was proposed and studied in [38] In model (1), the term βSI(t − τ )/[1 + αI(t − τ )] indicates that the infection at time t is caused by the infectious group I, infected τ units of time earlier. Hence, τ is the average incubation time of virus. In this paper, we shall deal with the intrinsic and extrinsic incubation times of dengue virus as in model (1). We shall also investigate the threshold dynamics of the delayed model established in terms of the basic reproduction number R 0 , which is defined as the average number of secondary cases induced by a typical infectious individual in a completely susceptible population [1]. The basic reproduction number is a key parameter both in the theoretical studies and the practice to control the transmission of infectious diseases [15,16,22], since the value of R 0 is associated with the transmission potential of a disease, the possible maximum number of the infected individuals and the strength of the control measures for a disease [27]. Usually, the disease will die out gradually in a population if R 0 < 1, while the disease will be endemic in the population if R 0 > 1. Indeed, we shall prove that this general understanding also holds true in the delayed model established.
In view of the importance of R 0 in threshold dynamics, it is clear that to control the prevalence of a disease, measures should be taken to reduce the value of R 0 to less than 1. Moreover, it is a well-known result that a disease will cease its prevalence if the vaccination rate for the susceptible is larger than (1 − 1/R 0 ) [16,22,49]. In the case of dengue, vaccine is on its way and vaccination will be possible in the near future. It is therefore very important to estimate the value of R 0 since vaccination is the most commonly used method to control disease outbreak. This leads to the second aim of this paper which is to propose a method to estimate R 0 of dengue transmission.
It is difficult to use the theoretical expression of R 0 derived from mathematical models to estimate its value directly. This is because, though some general methods [50,55] yield theoretical expression of R 0 in terms of the parameters of the models, in most cases some of these parameters are very difficult to estimate in practice. For example, the total number of mosquitoes is required -it is a big challenge and perhaps even impossible to get the exact number of certain types of mosquitoes in a region.
On the other hand, there is readily available information about a disease outbreak such as the number of newly infected cases announced by disease surveillance organizations daily or weekly. It is therefore more appropriate and fruitful to look for methods that employ these data to estimate R 0 of the disease. In fact, there are formulas to estimate R 0 from the initial growth data of a disease [1,23,34,35,54]. These formulas are developed under the scenario of an early outbreak of a disease, i.e. it is assumed that the whole population is susceptible to the disease. Therefore, for diseases, such as dengue, that have been prevalent in a population for decades and only a portion, not the whole, of the population is susceptible to these diseases [19], the application of formulas developed under a different scenario will inevitably cause errors.
Hereafter, we shall call 'R 0 ' in a partially susceptible population as 'reproduction number' (as in [30]) and denote it by R t . It is noted that in [30], the reproduction number is estimated via Richards model and the parameters of the model are obtained by curve fitting. In this paper, by regarding the number of infected cases periodically notified by surveillance organizations as the solution of an approximate mathematical model, we shall derive a formula to estimate the reproduction number R t of dengue transmission in a partially susceptible population. Our technique is inspired by the work of [23] and is different from the methods in [30,45]. Further, to illustrate the usefulness of the new estimation formula, we shall apply it to the dengue outbreak scenarios of Singapore in the years 2013-2015. Note that dengue has been prevalent in Singapore for decades and so the population is partially susceptible to the disease [19]. We remark that our estimation method can be generalized to estimating R t of other infectious diseases in a partially susceptible population.
The paper is organized as follows. In Section 2, a dengue mathematical model with two delays is established, then the basic reproduction number R 0 is defined and the threshold dynamics of the model is investigated. In Section 3, an approximate model is given based on the delayed model established earlier, then an estimation formula for the reproduction number R t in a partially susceptible population is derived. As an application, the R t values of dengue transmission in Singapore in the years 2013-2015 are obtained. Finally, we present the conclusion and discussion in Section 4.

Model with two delays
We first propose the following model with two delays Model (2) is based on the compartmental model (without delays) for dengue transmission established in [13]. Comparing with the model in [13], (i) we have introduced time delays τ 1 and τ 2 , which are the average extrinsic incubation time and intrinsic incubation time, respectively; (ii) the natural immunity from dengue and the disease induced death of human are neglected; and (iii) the average biting rates of susceptible and infected vectors are assumed to be the same.
In model (2), N H and N V are, respectively, the total number of human and mosquito population in a region, H is the recruitment rate of human population, μ H is the natural death rate of human, γ is the recovery rate of human from infection, and μ V is the death/birth rate of mosquito. Further, c(b) is the average transmission probability of an infectious human (mosquito) to a susceptible vector (human) and a is the average biting rate of vectors. All the parameters in model (2) are positive. We note that in model (2) the general incidence rate is adopted and two time delays τ 1 (average extrinsic incubation time) and τ 2 (average intrinsic incubation time) are introduced. The term S H (t) is the survival rate of mosquito (human) since not all those infected will survive after the time τ 1 (τ 2 ).
The human population N H is a variable since the fourth equation of (2), which is the sum of the first three equations of (2), gives On the other hand, since the birth rate and the death rate of the mosquito population are equal, N V = S V + I V is a constant in (2). It follows that and the sixth equation of model (2) can be rewritten as In view of (3) and (4), model (2) can be reduced to the following 4-dimensional system The dynamics of (6) is equivalent to that of (2). The initial values of (6) are given as where (t) and (t) are continuous in [−τ , 0]. Our first theorem shows that the biological model (6) with initial values (7) is wellposed. Theorem 2.1: All the solutions of (6) with (7) are positive and bounded.
Proof: Solving the third equation of (6), we obtain Hence, N H (t) > 0 for t > 0 and lim t→+∞ N H (t) = H /μ H . So N H (t) is positive and is bounded. Next, we shall prove that S H (t) > 0 for t > 0. Suppose the contrary and there exists t 1 such that S H (t 1 ) = 0 and dS H (t 1 )/dt ≤ 0. However, from the first equation of (6) we get dS H (t 1 )/dt = H > 0. This contradiction shows that S H (t) > 0 for t > 0.
We shall proceed to prove that I H (t) > 0 for t > 0. Suppose the contrary and there exists the smallest t 2 such that I H (t 2 ) = 0 and dI H (t 2 )/dt ≤ 0. It follows from the second equation of (6) that Therefore, we should have I V (t 2 − τ 1 ) ≤ 0. Since the solutions of (6) with (7) are continuous, there exists the smallest Then, all the solutions of (6) with (7) will enter and remain in D as t → +∞. For convenience, we now adopt the following change of variables Model (6) is then transformed to the following system and the initial values (7) are transformed to where t ∈ [−τ , 0], τ = max{τ 1 , τ 2 } and s H0 > 0, n H0 > 0, φ(t) > 0 and 0 < ψ(t) < 1.

R 0 and equilibria
By the survival method [28,51], the basic reproduction number R 0 of model (9) is given as where In (12), R 01 reflects the average number of mosquitoes infected by one typical infected person during his/her infectious period and R 02 describes the average number of human infected by an infected mosquito. Thus, R 0 is the average number of persons in the secondary class who are newly infected by one typical infected person in the first class. We shall investigate the role of R 0 in the equilibrium of (9), which is denoted byÊ = (ŝ H ,î H ,n H ,î V ). To proceed, we shall obtainÊ by solving the following algebraic system for It is immediate from the third equation of (13) that Next, summing the first two equations of (13) provideŝ while using (14) in the fourth equation of (13) giveŝ Now, substituting (14), (15) and (16) into the second equation of (13), we obtain We have the following two possible cases for the solutions ofî H . Case 1.î H = 0. It follows from (15) and (16) From the expression (17), it is clear thatî H > 0 if and only if R 0 > 1. Further, from (15) and (16), we haveŝ where R 01 and R 02 are given in (12). The above discussion results in the following theorem.

Stability of equilibria
In this section, we shall investigate the threshold dynamics of (9) with (10). In fact, we shall show that the disease-free equilibrium E 0 = (1, 0, 1, 0) is locally asymptotically stable when R 0 < 1 and the positive equilibrium Our method is to analyse the eigenvalues of the Jacobian matrix of the linearized system of (9) at the corresponding equilibriumÊ = (ŝ H ,î H ,n H ,î V ). Denote this Jacobian matrix by J(Ê), the expression of which is given in A.1 of the Appendix. The Jacobian matrix method has also been employed in [47].

Theorem 2.5:
The disease-free equilibrium E 0 of model (9) with (10) is locally asymptotically stable when R 0 < 1 and it is unstable when R 0 > 1.
Proof: Let J(E 0 ) be the Jacobian matrix of the linearized system of model (9) at the disease-free equilibrium E 0 , and let I be the identity matrix with suitable order. By direct computation, we have Clearly, J(E 0 ) has a real negative eigenvalue −μ H (multiplicity 2). The other two eigenvalues of J(E 0 ), denoted by λ 1 and λ 2 , are the roots of the equation Suppose R 0 < 1. When τ 1 = τ 2 = 0, it is clear that the real parts of the roots λ 1 and λ 2 of (19) are both negative since λ 1 + λ 2 < 0 and λ 1 λ 2 > 0. We show further in the following that no pure imaginary root of (19) exists if τ 1 > 0 or τ 2 > 0 or both. Suppose the contrary and (19) has a pure imaginary root λ = iω, ω > 0. Substituting λ into (19) and equating the real and imaginary parts to zero, respectively, yield It follows that ω > 0 is the root of However, it is obvious that (20) has no such positive root ω when R 0 < 1. This contradiction implies that (19) has no pure imaginary root if τ 1 > 0 or τ 2 > 0 or both. Note that f (λ, τ 1 , τ 2 ) is analytic, together with the fact that the real parts of the roots of (19) are all negative when τ 1 = τ 2 = 0, it implies that all the roots of (19) remain in the left half plane for all τ 1 > 0 or τ 2 > 0 or both. Hence, the real parts of all the roots of (19) are negative for any τ 1 ≥ 0 and τ 2 ≥ 0. This shows the local asymptotical stability of E 0 when R 0 < 1.
This implies that (19) always has a positive real root for any τ 1 ≥ 0 and τ 2 ≥ 0. This shows that the disease-free equilibrium E 0 is unstable when R 0 > 1.
Having proved that the disease-free equilibrium E 0 is locally asymptotically stable when R 0 < 1, we now proceed to analyse the stability of the positive equilibrium E * .

Lemma 2.6:
Proof: Under the special case τ 1 = τ 2 = 0, the model (9) reduces to the model studied in [13] when p = 0 and α H = 0. It is proved in [13] (Theorem 2(ii)) that all the eigenvalues of the Jacobian matrix of the linearized system of (9) at E * , J(E * ), have negative real parts when R 0 > 1. Hence, E * is locally asymptotically stable when R 0 > 1.
For convenience, we list the expressions of the constants A ij and B ik , i = 1,2, j = 1, 2, 3, k = 1, 2, which appear in subsequent lemmas, in A.2 and A.3 of the Appendix. Lemma 2.7: Let τ 2 = 0. Further, suppose that and Then, E * exists and it is locally asymptotically stable.
Proof: Note that (21) is equivalent to R 0 > 1 in the case of τ 2 = 0, hence E * exists by Theorem 2.4. Since the case τ 1 = τ 2 = 0 has been considered in Lemma 2.6, it remains to tackle the case when τ 2 = 0 and 0 < τ 1 < τ * 1 . To proceed, we shall analyse the local asymptotical stability of E * via the eigenvalues of the Jacobian matrix J(E * ) of the linearized system of (9) at E * . We have where A 1j and B 1k , j = 1,2,3, k = 1,2, are constants independent of λ (refer to Appendix A.2 for details). It follows that J(E * ) has a negative eigenvalue −μ H and the other eigenvalues are determined by the equation Next, we shall prove that (23) has no pure imaginary root when τ 2 = 0 and 0 < τ 1 < τ * 1 . Suppose the contrary and (23) has a pure imaginary root iω, ω > 0. Using a similar technique as in the proof of Theorem 2.5, we get Denote z = ω 2 , and rewrite (24) as It follows from (25) that g 1 (0) = A 2 13 − B 2 12 ≥ 0 (from the first inequality of (22)), and in view of the second inequality of (22), we find for z > 0. Noting that there is at most one value of z satisfying dg 1 /dz = 0, we have g 1 (z) > 0 for z > 0. Hence, (25) has no positive real root. This contradiction implies that (23) has no pure imaginary root when τ 2 = 0 and 0 < τ 1 < τ * 1 . Moreover, since the function involved in (23) is analytic and all the eigenvalues of J(E * ) have negative real parts when τ 1 = τ 2 = 0 (from the proof of Lemma 2.6), all the roots of (23) remain in the left half plane when τ 2 = 0 and 0 < τ 1 < τ * 1 . This shows that E * is locally asymptotically stable.
Similar to the proof of Lemma 2.7, we establish the next result. Lemma 2.9: Let τ 1 = 0. Further, suppose that and Then, E * exists and it is locally asymptotically stable. Theorem 2.10: Let (21), (22) and (27) be satisfied. Further, suppose that Then, E * exists and it is locally asymptotically stable.

Estimation formula for R t
Recall that the basic reproduction number 'R 0 ' in a partially susceptible population is termed reproduction number and is denoted by R t . In this section, we shall derive an estimation formula for the reproduction number R t of dengue transmission in a partially susceptible population.
To this aim, we assume in model (2) that the total number of the human population N H = S H + I H + R H is constant for all time t. Therefore, the linear recruitment rate of human, λ H , is equal to the natural death rate μ H . The model (2) is simplified to the following: The biological implications and the values of the parameters in model (31) are listed in Table 1. From (11), the basic reproduction number R 0 of model (31) is now given as Clearly, if one employs (32) to calculate R 0 , one needs to estimate the value of N V , a, b and c. As mentioned in the Introduction, it is a big challenge to estimate accurately the total number of Aedes aegypti in the region (N V ) and ubiquitous difficulties also occur in the Intrinsic incubation period of dengue virus 1-11 days [23] estimations of a, b and c. Therefore, we shall develop a new estimation formula that does not involve parameters which are difficult to determine. Indeed, we shall give a formula to estimate the reproduction number R t (instead of R 0 ) when the human population is partially susceptible to the disease. In standard practice, surveillance organizations will report newly infected cases daily or weekly to the public. The newly infected cases are clearly related to the solutions of model (31) which in turn would involve parameters such as N V , a, b and c. Therefore, instead of estimating these parameters directly, we can consider the solutions of model (31) with a view to using the information of newly infected cases.
As the newly infected cases are related to the dynamics of I H , we rewrite the second and the last equations of (31) as It is well known that partial immunity is gained when the infected individual recovers from dengue infection. Therefore, not all the human population is susceptible to certain serotype of dengue virus if the disease is prevalent in a region for decades. As an example, only less than 60% of the population in Singapore is susceptible to dengue [19]. Hence, we assume that S H = N H − I H − R H ≈ ηN H where η is the proportion of the susceptible in the whole human population at time t. Also, we set S V = N V − I V ≈ N V due to the relatively short lifespan of the mosquitoes. Consequently, from (33) we get an approximate model of the disease evolution where m = N V /N H . Note that the parameter η in model (34) indicates that the human population is partially susceptible.
The following method to derive the estimation formula for R t is inspired by the work of [23] where the human population is assumed to be completely susceptible -different from our assumption. To proceed, we assume that the cumulative number K of infected cases is exponential [1,23], i.e. K = C exp(λt), where λ is the force of infection in a partially susceptible population and C is a positive constant. Since the cumulative cases correspond to the integration of I H , it follows that I H is also exponential with the same rate λ. As such, we can assume that I H and I V take the form where C 0 and C 1 are positive constants.
Substituting (35) into (34), we have The first and the second equations of (36), respectively, give and Combining (32), (37) and (38), we get the estimation formula of R t in a partially susceptible population as

Remark 3.1:
The new formula (39) will be used to estimate R t . We observe that the presence of the incubation times τ 1 and τ 2 increases (decreases) the value of R t when λ > 0 (λ < 0), and the presence of the proportion of susceptible η increases the value of R t . As (39) does not involve the parameters N V , a, b and c, so the difficult problem of estimating these parameters is solved. Moreover, the force of infection λ can be obtained from the infected cases announced by the surveillance organizations through appropriate curving fitting. The value of η can be estimated by statistical method such as that proposed in [19]. Other parameters can also be determined from known data available in the literature.

Estimating R t : the singapore case
In this section, we shall illustrate the estimation process of R t via (39) by employing the factual data of dengue cases in Singapore in the years 2013, 2014 and 2015. The newly infected cases of dengue in Singapore are reported by the Ministry of Health on a weekly basis [40]. The data of 2013-2015 are depicted in Figure 1. In Figure 1, the unit of time (horizontal) axis is 'e-week', which is the abbreviation of 'epidemiological week', that consists of 7 days from Sunday to Saturday. The first e-week of year 2013 is from 30 December 2012 (Sunday) to 5 January 2013 (Saturday). For more details, one can refer to [40].
Using the dengue data of 2013, we shall present the detailed estimation process and thereby obtain the R t value for the year 2013. By a similar process, we shall also obtain the R t values for the years 2014 and 2015.

Year 2013
We first calculate the cumulative weekly dengue cases in 2013. Next, the curve fitting toolbox of Matlab is applied to fit the cumulative cases with the function c 1 exp(c 2 t) and the time unit is rescaled to days. Let the interval i, j denote the period from the ith e-week to the jth e-week. In practice, it is more important to estimate the R t value during a disease outbreak. Hence, noting that the 2013 outbreak occurred from the first e-week, we perform the curve fitting in the period 1, n (for different values of n) and record the associated Rsquare value. R-square, a value between 0 and 1, is a descriptive statistic of Matlab curve fitting toolbox that evaluates the goodness of fit, the larger R-square is, the more accurate is the fit. Hence, we shall choose the fit with the largest R-square to obtain the estimated value of λ (= c 2 ).
For illustration, a sample of R-square values of curve fitting in 1, n of 2013 is listed in Table 2. In fact, the largest R-square value of 0.9883 (Matlab R2011a) is achieved when n = 26, and the curve fitting is depicted in Figure 2. Thus, we obtain λ = 0.01429 (±0.00085, 95% confidence interval).
We set μ H = 3.32 × 10 −5 days −1 (life expectancy statistic of Singapore in 2013 [39]), γ = 1/6 days −1 [23] and μ V = 0.064 days −1 [8]. Knowing that the extrinsic incubation  period varies from 1 to 16 days and the intrinsic incubation period varies from 1 to 11 days [23], we take the mean and set the extrinsic incubation period τ 1 = 8.5 days and the intrinsic incubation period τ 2 = 6 days. It is observed in [19] that less than 60% of the population in Singapore was susceptible to dengue. Therefore, we set the proportion of the susceptible in the human population η = 60%. Now all the parameters in (39) are estimated, and it follows immediately from (39) that the reproduction number for the year 2013 is

Years 2014 and 2015
We shall apply the above method to estimate R t of dengue transmission in Singapore in the years 2014 and 2015. For the year 2014, the curve fitting is carried out in the period 9, n (for different values of n) since the dengue outbreak began from the 9th e-week in 2014 (refer to Figure 1). It is found that curve fitting in the period 9, 31 attains the largest R-square value of 0.9956, and we obtain λ = 0.009767 (±0.000324, 95% confidence interval). Figure 3(a) illustrates the curve fitting in 9, 31 . With γ , μ H , μ V , τ 1 , τ 2 and η taking the same values as before, formula (39) gives the reproduction number for the year 2014 as Similarly, for the year 2015, the curve fitting is carried out in the period 11, n (for different values of n) since the dengue outbreak began from the 11th e-week in 2015 (refer to Figure 1). As a result, λ = 0.006919 (±0.000113, 95% confidence interval) is obtained by curve fitting in the period 11, 24 which has the largest R-square value of 0.9994. Figure 3(b) illustrates the curve fitting in 11, 24 . Again with γ , μ H , μ V , τ 1 , τ 2 and η taking the same values as before, formula (39) gives the reproduction number for the year 2015 as

Dynamical estimation of R t
We note that the above estimation process of R t in a particular year is carried out after the whole year's data are available and the outbreak starting point is observed from the data. Therefore, the estimated R t value is kind of 'on hindsight' and may not be timely enough for decision making; moreover, the outbreak starting point is sometimes difficult to determine from the data. This leads to the necessity of estimating R t dynamically, i.e. to estimate R t over a moving short period of time, so as to provide more timely information for disease control. It should be realized that formula (39) can be applied to obtain the dynamical estimation of R t in any arbitrarily given period. As an example, suppose we fix the period as 6 e-weeks, then the force of infection λ in each time interval i, i + 5 (i ≥ 1) can be obtained via curve fitting on the infected cases of the corresponding period, subsequently the value of R t over i, i + 5 is estimated via (39). Indeed, with period fixed as 6 e-weeks, using the Singapore

Conclusion and discussion
In view of the considerable amount of incubation time presents in both human and mosquito during dengue transmission, in Section 2, we have established a mathematical model with two delays that takes into account the extrinsic incubation time (τ 1 ) and the intrinsic incubation time (τ 2 ). We have shown that the solutions of the delayed model are positive, and the model is dissipative with a nonnegative invariant set which is bounded. The basic reproduction number R 0 is defined and the equilibria of the delayed model are calculated. We have investigated the threshold dynamics of the delayed model and it has been established that the disease-free equilibrium E 0 is locally asymptotically stable when R 0 < 1, and the positive equilibrium E * exists and is locally asymptotically stable when Recalling that the definition of basic reproduction number R 0 involves a completely susceptible population, we term 'R 0 ' in a partially susceptible population as ' reproduction number' and denote it by R t . In Section 3, we have derived an estimation formula for the reproduction number R t of dengue transmission in a population (such as Singapore) that is not completely susceptible to the disease -only a fraction (η) of the population is susceptible. Our methodology in Section 3 involves introducing a simplified delayed model of that established earlier in Section 2 and recognizing the solutions of this simplified model as the factual data of dengue cases announced by surveillance organizations at regular time intervals. As a result, a new R t expression (39) of the simplified delayed model is obtained. We remark that our method can be generalized to derive estimate for the 'reproduction number' of other infectious diseases in a partially susceptible population. To illustrate the application of the new estimation formula (39), we have employed the dengue data of Singapore in the years 2013, 2014 and 2015 and performed appropriate curve fitting to obtain R t values of 2.72, 2.34 and 2.13 for the years 2013, 2014 and 2015, respectively. The estimation process is demonstrated to be feasible, convenient to apply and efficient.
From the estimation formula (39) for R t , as noted in Remark 3.1, we know that when the extrinsic or intrinsic incubation period is considered (τ 1 , τ 2 > 0), it will increase (decrease) the estimated value of R t when the force of infection λ > 0 (λ < 0); also a partially susceptible population (η < 1) will increase the estimated R t . Taking the year 2013 in Singapore as an example, a smaller R t = 2.21 (versus 2.72) is obtained when the extrinsic and intrinsic incubation periods are not considered (τ 1 = τ 2 = 0); also a smaller R t = 1.63 (versus 2.72) is obtained if we assume that the population is completely susceptible to dengue (η = 1).
In the estimation process, the parameters involved in (39) are, respectively, set by their mean values. In fact, an interval estimate of R t (rather than a point estimate) will be better since the values of the parameters in (39) are in intervals (see Table 1) and the estimation of λ by Matlab is also an interval. As an example, an interval estimate of [2.05, 3.61] is obtained for R t in the year 2013 by applying the values of parameters given in Table 1 and the 95% confidence interval [0.01344, 0.01514] of λ.
From Section 3, we know that the estimated values of R t can be different in different time intervals i, j , i.e. R t varies with respect to time. Hence, the value of R t can be estimated over a moving period of time, such as i, i + 5 (i ≥ 1) -we call this dynamical estimation. The dynamical estimation can reflect the variation of R t with respect to time.
In the literature, the reproduction number of dengue transmission in Singapore is estimated to be 2.23 in the year 2005 [30], and 1.9 in the year 2002 [36]. Note that in [36], the population is not mentioned to be partially susceptible and the so-called 'equilibrium method' is used in getting R 0 = 1/s * , where s * = s H (∞)/N H (∞); while in [30], the population is considered to be partially susceptible and the cumulative dengue cases of year 2005 in Singapore during the period 17, 52 are regarded as the solution of suitable Richards model, then the parameters of the model are obtained by curve fitting. Further, the estimated basic reproduction number for different cities of Brazil varies from 2.26 to 11.39 in [12], and from 1.3 to 11.6 in [43]. Thus, comparing with these results, the estimated R t values of 2.72, 2.34 and 2.13 obtained in this paper for the dengue transmission in Singapore in the years 2013, 2014 and 2015 are within the range.
In the process of applying (39) to estimate the reproduction number of dengue outbreak, we observe that the natural death rate (mortality rate) μ V of the vector, Aedes aegypti, varies greatly in the literature. Different values were applied by different researchers in numerical simulation. For example, μ V = 0.025 days −1 in [37], μ V = 0.9 days −1 in [17], μ V = 0.156 weeks −1 in [12,35], and μ V = 1/10.5 days −1 in [11]. In Section 3, we have set μ V = 0.064 days −1 which is the mean of laboratory experimental values that vary from 0.05 to 0.082 days −1 [8]. In [24], a DENSim model, based on experiments, was proposed to estimate the mortality rate (μ V ) of dengue vector and the extrinsic incubation time (τ 1 ). Many factors, such as temperature, vapour pressure and relative humidity, were involved in the model. This model was applied in [23] to obtain the corresponding mortality rate of dengue vector and extrinsic incubation period of dengue virus in different cities of Brazil. We remark that, to improve the estimation accuracy of R t in certain regions, not only laboratory experiments but also field experiments should be carried out to estimate μ V and τ 1 , since it is pointed out in [8] that the mortality rate of Aedes aegypti in field experiments may be quite different from that in laboratory experiments.
In the simplified model (34) for dengue transmission, we introduce the partial immunity of human population due to previous infection that is described by the parameter η. The formula (39) for the reproduction number R t in a partially susceptible population is thus deduced. Vaccination will be a new control method in the near future since dengue vaccine is on its way to be developed and this method will affect the value of η. The formula (39) makes it possible to compare the effectiveness of vaccination with the traditional control method of vector control [2,9]. The sensitivity of R t with respect to the parameters η and μ V is respectively From (43) we know that the value of R t will decrease 1% if the vaccinated human population increases 1%, and the value of R t will decrease λ/(λ + μ V )% if the death rate of the vector increases 1%. Recalling from Section 3, the estimated value of λ in Singapore for years 2013-2015 is about 0.01 while the value of μ V is about 0.06, hence λ/(λ + μ V ) ≈ 1/7. This means vaccination (if possible) is more sensitive to decrease the value of R t than vector control in Singapore for years 2013-2015 (about 7 times). In fact, as seen from (43) and (44), R t is more sensitive to vaccination since λ/(λ + μ V ) < 1 during the prevalent period of dengue. Supposeη is the smallest vaccination rate to reduce the value of R t at present to a value less than 1, from (43) we have R t (1 −η) < 1, which yields η > 1 − 1/R t . This value ofη implie theoretically the smallest vaccination rate to eradicate the dengue transmission in a region. Lastly, we note that it is more appropriate to incorporate partial immunity in the analysis of model (2) and its corresponding non-dimensional model (9) for two reasons [13]: (i) partial immunity can be gained due to previous infection, and (ii) the human population may get immunity via dengue vaccine in the near future. Also, models with variable human and vector population will be more suitable for diseases, such as dengue, that have a long epidemic cycle. Global asymptotical stability [53,58] of the disease-free equilibrium E 0 and the epidemic equilibrium E * of model (9) can be further analysed. It is desirable in practice to be able to forecast the future reproduction number based on the data obtained earlier, but our estimation method for R t cannot be used to forecast the future variation of reproduction number of dengue transmission. We remain all the above as future work.

Disclosure statement
No potential conflict of interest was reported by the authors.