INTRODUCTION
The purpose of this paper is to examine the role of mathematical models in evaluating control strategies for diseases, illuminating their limitations and possible errors that can occur if these limitations are not carefully considered.
For diseases involving vectors, the number of variables can be large and it is sometimes unnecessary and misleading to use all the complexities of the real biological system.
As noted by Mazilu et al. [1], this is far from a trivial matter: ‘Building simple models that capture the essence of a physical phenomenon is not an easy task: there is a fine line between an accurate model that is as simple as possible, and an unrealistic model.’
This paper is arranged in three main sections.
In the following section (General considerations about modelling vectorborne infections), we show how, considering the purpose of the model, we can simplify the description of a complicated system as a vectorborne disease. We then illustrate how to limit the model variables according to our needs. In so doing, we aim to explain why there are many models of the same (or similar) systems with different numbers of variables [2–4].
The next section (Estimating R
_{0} from the initial growing phase of an outbreak: several pitfalls), examines one example found in a recent paper [5] that illustrates the dangers of basing control strategies on models without considering their limitations. We show that the method used cannot be applied as done by the authors, and that as a consequence, some of the conclusions reached are not correct. Unfortunately, the incorrect results, if taken as true by public health authorities, would greatly harm the affected populations and the reputation of the use of models in this field.
Finally, in the ‘Backward bifurcation’ section, we review an interesting paper by Garba et al. [6] that contains some biological assumptions that are inappropriate for dengue but may apply to other vectorborne infections. In addition, it contains some algebraic misprints that make it difficult to read. The paper's conclusion is that a backward bifurcation exists in the system; this is examined in our paper from a different point of view. Furthermore, the conclusions of Garba et al., which do not apply to dengue, would make the control of even those diseases to which Garba's parameters do apply extremely difficult. We argue that the values of the biological parameters used by Garba et al. are not realistic for dengue.
We will show in the ‘Estimating R
_{0} …’ section that the variables in the equations of the models are densities, i.e. vectors/humans per unit area. Therefore, if these densities are spatially homogeneous and the area investigated is small, total numbers can be obtained by multiplying the variables by the area of the region under consideration. However, when we investigate the outbreak of an epidemic, we should consider that the initial infection is not distributed homogeneously throughout the area under consideration. In the literature, this is often disregarded and epidemics are investigated by assuming that the disease invades an area homogeneously. In Appendix I, we show the effects of the initial distribution of the disease on the size and duration of the epidemic.
GENERAL CONSIDERATIONS ABOUT MODELLING VECTORBORNE INFECTIONS
We will use dengue as an example of a vectorborne disease to illustrate points that are also valid for other vectorborne diseases, such as malaria and yellow fever. We do this to simplify the notation and make the points we want to emphasize simpler to understand.
In the dengue system, there are several populations involved: a human population, an adult mosquito population, at least six aquatic forms and mosquito eggs. In addition, there is the virus, which in dengue's case will be one of four known serotypes.
Our first point is the question of how to choose the populations that will enter the model. In the seminal papers by Ross [2] and Macdonald [3] (that address malaria), only adult mosquitoes, humans and parasites were considered. This was because they wanted to investigate the existence of critical sizes of these populations, below which the disease would disappear. In the mathematical literature, a model exhibiting this feature is said to have a threshold and, in the case of diseases, this threshold is given by the wellknown basic reproduction number, R
_{0} [7]. Let us briefly recall a slightly generalized form of the classical Ross–Macdonald model [2, 3]. The variables are described in Table 1 and are usually denoted as ‘compartments’ of the model.
Table 1. Variables of the model and their biological description
The equations that govern the dynamics of the system are given below and explained immediately afterwards.
The model parameters and their biological interpretation are given in Table 2.
Table 2. Model parameters and their biological interpretation
Let us explain the meaning and limitations of the above model. First, as explained, the variables are densities, i.e. the number of humans/vectors per unit area. Therefore, to use the model as it is written above, we should consider an area where the population is approximately homogenously distributed and multiply each variable by this area. A qualitatively comprehensive attempt to discuss the challenge of calculating parameters taking into account spatial heterogeneities can be found in [8]. One particularly important point is raised by the term
Let us explain the meaning of this term. The parameter a is a composed quantity. Let a be the area explored by a mosquito via the joint movement of humans and mosquitoes. Let ξ be the number of bites a mosquito inflicts per unit time and per unit area in the human population. Then, ξ
AI
_{
M
} is the number of bites that AI
_{
M
} infected mosquitoes inflict on N
_{
H
}
A people. Hence, the fraction of bites inflicted on susceptible humans is ξAI
_{
M
}(S
_{
H
}
A/N_{H}A) = αI
_{
M
}(S
_{
H
}
/N_{H}
), where
Therefore, the number of susceptible humans that acquire the infection from infected mosquitoes per unit time is abI
_{
M
}(S
_{
H
}/N
_{
H
}), where b is the probability that a bite from an infected mosquito results in an infected (latent) human. Assuming that latent mosquitoes also transmit the infection, although with a lower probability, expressed by η
_{
M
}, we see that equation (2) represents the density of new infections per unit time due to mosquito bites.
Analogously, the term
$ac{\textstyle{{({I_H} + {\eta _H}{L_H})} \over {{N_H}}}}{S_M}$
represents the density of new infections per unit time in mosquitoes due to mosquito bites on infective humans.
The two quantities Λ_{
H
} and Λ_{
M
} are the number of humans and vectors born or otherwise introduced per unit area per unit time into the susceptible compartments. It is usual to assume that Λ_{
H
} is a logistic term of the type
${r_H}{N_H}(1  {\textstyle{{{N_H}} \over {{K_H}}}})$
, where r
_{
H
} is the Malthusian parameter and K
_{
H
} is the carrying capacity. Analogously, Λ_{
M
} is the quantity of vectors that are born or otherwise introduced per unit area per unit time.
The other terms are transition terms between the compartments as explained, for example, in [9].
From equation (1), we can deduce [4, 9, 10] that the disease cannot invade the host population if R
_{0} is <1, where R
_{0} is
The terms f
_{
H
} and f
_{
M
} are defined as
In the limit when η
_{
M
} = 0 and δ
_{
H
}→∞, we obtain the expression of R
_{0} of the original Macdonald model [3]:
The above result is very important because it indicates that we do not have to completely eradicate the mosquito population to prevent the disease from invading the host population [2]. Note, however, that the model is valid only when applied to a region where the populations involved (mosquitoes and humans) are spatially homogeneous, and in any case, border effects were neglected. We shall expand on this later in the paper.
In some regions of the world, seasonality affects the mosquito population considerably and, in some cases, the mosquito density falls so much that transmission is interrupted in the dry season. The disease, however, overwinters and reappears in the next season. A hypothesis to explain this phenomenon is transovarian transmission in the mosquitoes [4]. In this case, we are therefore forced to consider in the model the aquatic stages of the vector. Control measures that aim for the destruction of breeding places also force us to include aquatic forms in the model, even if seasonality is negligible. The model described in equation (1) has to be modified to include the immature forms. This can be done by modifying the last three equations related to the vectors in system (1) as follows:
First, we introduce two new system variables (compartments) representing noninfected and infected aquatic forms: S
_{
E
} and I
_{
E
}, respectively. For simplicity, we consider all the aquatic forms, i.e. eggs, larvae and pupae, as a single compartment (we use a subscript E to represent eggs). The new system reads
The new parameters in system (7) are described in Table 3.
Table 3. Parameters and their biological meaning in the model with aquatic forms
Let us now describe the meaning of some terms of system (7). The term pc
_{
s
}(t)S
_{
E
} represents the number of noninfected eggs per unit area per unit time that reach the adult stage. The parameter c
_{
s
}(t) was introduced to mimic seasonality (see [4]). The term pc
_{
s
}(t)I
_{
E
} represents the number of infected eggs per unit area per unit time that reach the infected adult stage. The term
$[{r_M}{S_M} + (1  g){r_M}({I_M} + {L_M})]\left( {1  {\textstyle{{{S_E} + {I_E}} \over {{K_E}}}}} \right)$
is a logisticlike term that represents the rate of mosquito reproduction in the form of noninfected eggs. Similarly,
$[g{r_M}({I_M} + {L_M})]\left( {1  {\textstyle{{{S_E} + {I_E}} \over {{K_E}}}}} \right)$
is a logisticlike term that represents the rate of mosquito reproduction in the form of infected eggs. K
_{
E
} is the carrying capacity of the environment. The remaining terms are transition rates between the compartments. This modified system also predicts a threshold for the infection to invade the host population (for details, see [4]). Note that we assume that the mortality of infected eggs is the same as the mortality of noninfected eggs. This is easily modified but, in the spirit of avoiding unnecessary experimentally unknown facts, we prefer to assume them equal. The proportional distribution of infected eggs into noninfected and infected mosquitoes was briefly discussed in [4]. In this reference, it was shown that this parameter is very important to explain semiquantitatively dengue overwintering.
ESTIMATING R
_{0} FROM THE INITIAL GROWING PHASE OF AN OUTBREAK: SEVERAL PITFALLS
From the discussion in the previous section, we can see that calculating R
_{0} for a given population is an important task. This can be done in several ways. A simple way first proposed by May & Anderson [11] for directly transmitted infection was adapted by Massad et al. [12] and further refined by Favier et al. [13] for vectorborne infections. The first pitfall we will discuss was the use made by Pinho et al. [5] of this method as described below.
In their paper, Pinho et al. [5] analysed the dengue epidemic in Salvador, Brazil, and concluded that
The value of R
_{0} is greater than 1 for the epidemic in 1995–1996 for any chosen value of the vectorcontrol parameter, indicating that other strategies would be necessary besides the adult vectorcontrol, as such as the control of the mosquito's aquatic phase, to reduce its force of infection and therefore to control the epidemic.
The aim of this section is to show that this conclusion cannot be deduced from either the model they used or from the calculations presented in their paper.
In fact, the equilibrium solutions of the model given by their equations (2.1) can be deduced analytically [10]. It is also possible, and easier, to deduce two thresholds. The first (denoted ‘basic offspring’ or demographic reproduction number) determines the permanence or extinction of the mosquito population. The second is a variant of the classical Macdonald [3] basic reproduction number (R
_{0}) and determines the establishment or extinction of the infection. From these results, and also from numerical simulations of the model, it is easy to see that by increasing the adult mortality rate, μ
_{
M
}, or the aquatic phase mortality rate, μ
_{
a
}, we can reach Macdonald's threshold before reaching the basic offspring threshold. Therefore, the disease dies out without eliminating the mosquito population (see Fig. 1 below from Pinho et al.'s model). Thus, Macdonald's threshold is useful. It would be useless if Pinho et al. [5, pp. 5692] were correct.
Fig. 1. The values of c
_{
M
} such that the disease dies out but the mosquito population is still present.
The authors reached their conclusion by using their equation (3.8)
This equation, or rather a simplified version of it, was obtained in [13–15]. This technique assumes a force of infection Λ > 0 and, therefore, no value of the control parameter, c
_{
M
} (which Pinho et al. [5] added to the adult mortality rate, μ
_{
M
}), can reduce R
_{0} below unity. To investigate the case in which Λ < 0 and, therefore, R
_{0} < 1, another method must be used, as described in [15]. Therefore, Pinho et al.'s conclusion that ‘other strategies would be necessary besides the adult vectorcontrol […] to reduce the force of infection’ cannot be deduced from their equation (3.8) and is, in fact, wrong. Again, we stress that Macdonald's threshold would be useless if Pinho et al.'s [5] conclusion were true.
When R
_{0} < 1, as explained in [15, 16], although an outbreak can occur, it will naturally fade away.
We shall return to the above method of estimating R
_{0}. Here we only mention that this method is correct if the disease is introduced in a small portion of the region, as it usually is, that is going to be invaded. This is because the exponential growth is deduced by assuming that the disease is introduced in a uniform region. However, as we shall see in Appendix I, when introduced in just a small portion of the region, the development of the disease proceeds as a wave. Since what is usually reported is the number of cases per unit time (epidemiological week, for instance) this method can only be used if the disease propagates slowly in space because (see Appendix I) the entire epidemic curve is not exponential but is so only in the very first initial phase.
Finally, to end this section, we note that the sensitivity of the mortality rates of mosquito phases was investigated in [10, 17, 18]. From these papers, we conclude that the adult mosquito mortality rate is the parameter to which R
_{0} (and also the force of infection) is, by far, most sensitive. Compared with the aquatic phase mortality rate, this is rather intuitive because, in equilibrium, killing one mosquito prevents the appearance of hundreds of eggs and other aquatic phases, whereas killing one egg prevents the appearance of just one mosquito. Orders of magnitude of the relative effect between control measures directed against mosquitoes vs. control measures directed against eggs were given in [10]. In this paper, we compared the sensitivity of R
_{0}, the force of infection and the steadystate prevalence to the control parameters. In general, the sensitivity to control measures (increase of mortality rates) directed against mosquitoes is more efficient by 2–3 orders of magnitude. However, these calculations did not take into account the logistic and financial costs of any of the abovementioned strategies.
BACKWARD BIFURCATIONS
In an influential and interesting paper, Garba et al. [6] claim that a model very similar to the one given by equation (1) exhibits socalled backward bifurcation (subcritical bifurcation). In this section, we find the stationary solutions of system (1) and prefer not to call the observed phenomena backward bifurcation. We examine this very interesting system from a different point of view and also discuss the possible origin of the phenomena exhibited by Garba et al.'s model [6].
The stationary solution for the prevalence of the infection in humans, I
_{
H
}/N
_{
H
}, is obtained by replacing the derivatives on the right side of the system (1) equations and solving the resulting nonlinear system of equations. The results (see [10, 19] for more details) are:
which is the equilibrium prevalence of the infection in humans, where
and
To obtain equation (9), we assumed a logistic birth (or immigration) rate for Λ_{
H
}
where r
_{
H
} is the per capita birth (or immigration) rate of humans and K
_{
H
} is the carrying capacity of humans.
The stationary solution for the total human population, N
_{
H
}, is
where
and
The stationary solution for the specific case of dengue prevalence in humans, I
_{
H
}/N
_{
H
}, is obtained by making δ
_{
H
} → ∞ and θ
_{
H
} = σ
_{
H
} = 0 in equation (9)
The basic reproduction number can be deduced by making the numerator of the previous equation equal to zero, resulting in equation (5).
In the specific case of dengue, the value of N
_{
H
} reduces to:
where
and
Equation (17) does not show any backward bifurcation.
A possible misprint in Garba et al.'s [6] paper results from writing Λ
_{
H
} [their equation (1)] as
The correct equation should have N
_{
H
} instead of N
_{
M
} in the denominator because, as can be noted in system (1), the total number of bites inflicted by infected mosquitoes is a(I
_{
M
}+η
_{
M
}
L
_{
M
}), a fraction of which, S
_{
H
}/N
_{
H
}, are on susceptible humans, of which a fraction b is actually infective.
Garba et al.'s [6] assumption that the total number of bites that mosquitoes inflict in people is, of course, equal to the total number of bites that people receive is correct. However, the force of infection from mosquito to human, ab/N
_{
H
}(I
_{
M
} + η
_{
M
}
L_{M}
), depends on b, whereas the force of infection from human to mosquito, ac/N
_{
H
}(I
_{
H
} + η
_{
H
}
L_{H}
), depends on c. Therefore, equation (4) in Garba et al.'s [6] paper is correct only in the particular case when b = c. In addition, Garba et al. [6] assumed that susceptible and infected mosquitoes bite with different rates, although this is not explicitly used in their model. The points mentioned in the above remarks on Garba et al.'s [6] paper have no influence on the existence of backward bifurcation. We mention these facts for completeness and to describe the Garba et al. [6] paper fairly, i.e. describing all the biological realities that they introduce in their model.
To compare our model with that in the paper by Garba et al. [6], in more detail, we take, Λ_{
H
} = const., instead of the one given by equation (13), and make Λ_{
M
} also constant. The system can be solved, and the stationary state value of I
_{
H
}/N
_{
H
} (the translation to Garba et al.'s [6] notation is shown in Appendix II) is:
where
It is easy to demonstrate that one of the roots of equation (21) is negative. Making the positive solution equal to zero, we can deduce the threshold. In fact, the threshold (and therefore R
_{0}) can be deduced by making
$  {B_g} + \sqrt {B_g^2  4{A_g}{C_g}} = 0$
, i.e. C
_{
g
} = 0. The result is
This value of R
_{0} can also be deduced by linearizing the system (1), using Garba et al.'s [6] notation, around the trivial solution (no disease). The value of N
_{
M
}/N
_{
H
} is crucial. Garba et al. [6] use Λ_{
M
}
μ_{H}
/Λ_{
H
}
μ
_{M}
(see [6, p. 14]). However, to obtain the other points in their figure 2, it is necessary to increase the initial value of N
_{
M
} with respect to N
_{
H
}. This is what we think Garba et al. [6] did: they varied the parameters and calculated R
_{0} using N
_{
M
}
/N_{H}
= Λ_{
M
}
μ
_{H}
/Λ_{
H
}
μ
_{M}
, but introduced infected mosquitoes into the initial condition. Introducing infected mosquitoes into the initial condition changes R
_{0} by increasing N
_{
M
}/N
_{
H
}. Therefore, the values given by R
_{0} in Figure 2 are actually greater than 1.
Fig. 2. The three behaviours of the system described in the text.
Consider the value of R
_{0} [equation (23)]. Because the disease affects both populations involved, we must note that, if N
_{
M
}/N
_{
H
} calculated at time t = 0 makes R
_{0} > 1, then the disease will invade the population in the sense that immediately after 0 the values of I
_{
M
} and I
_{
H
} will be greater than their values at t = 0.
If the disease invades the population, two things can happen. If the value N
_{
M
}(t)/N
_{
H
}(t) is such that, at some value of t, the value of R
_{0} falls below 1, the disease will disappear. If this does not happen, the disease will go to the steady state. This depends on the initial condition. For the values given by Garba et al. [6], who do not use densities, the transition occurs for S
_{
V
}(0) = 2717.8445, E
_{
V
}(0) = 0, I
_{
V
}(0) = 100, S
_{
H
}(0) = 511.82, E
_{
H
}(0) = 0, I
_{
H
}(0) = 1, R
_{
H
}(0) = 0. Note that without disease, the equilibrium values are N
_{
V
} = 1875 and N
_{
H
} = 512.82. The resulting figure is shown below.
A summary of the above results will now be given. First, we define a function of t, R(t), as
We can understand the behaviour of this system by using this function as explained below.

(1) If R(t = 0) (which happens to be R
_{0}) is less than 1, then the disease will not invade (see the two lowest curves in Fig. 2);

(2) If R(t = 0) is greater than 1, the disease will invade the population, but we have to distinguish between two cases:

(2.1) If, for some t, R(t) drops below 1, then the disease disappears (see curves 3–7 in Figure 2)

(2.2) If R(t) is greater than 1 for all t, then the disease goes to a steady state different from zero (see curves 8–15 in Fig. 2)

(3) There is a threshold value for N
_{
M
}(0)/N
_{
H
}(0), below which case (2.1) occurs and above which case (2.2) occurs. We found this threshold numerically but we suspect that, by using Garba et al.'s [6] method, this can be obtained analytically.
In Appendix II, we show in detail why we prefer not to call the phenomena described above and in Garba et al.'s [6] paper a backward bifurcation, as is suggested by Garba et al. [6].
CONCLUSIONS
Mathematical models in sciences must have their complexity adjusted to their goals, and, as we have seen, we have basically two classes of models. At one extreme we have models that are intended to check if our intuition about why a certain phenomenon occurs is correct. This model must be as simple as possible, otherwise it may fail. An example of such a model is a model designed to test dengue overwintering [4]. At the other extreme, we have models whose goals are to predict future outcomes. These models are necessarily very complex. An example of such a model (for dengue vaccination) is given by [20]. Of course there are models in between these classes.
In this paper, the model belonging to the second class of models (to estimate the basic reproduction number from the initial phase of the infection) fails by not including sufficient complexities that are necessary to deal with the data.
We have seen in Appendix I why the effects of the spatial distribution of a vectorborne infection should be taken into account when estimating the basic reproduction number from the initial size of the epidemics. Another crucial effect that has to be taken into account is that the size and duration of the epidemic depend on the initial distribution of the invasion of the infection, i.e. where and how it was introduced in a given region previously free of the infection.
The model belonging to the first class that is examined in this paper is one that tries to clarify the role of the socalled backward bifurcation. In the ‘Backward bifurcations’ section, we examined in detail the paper by Garba et al. [6] which analyses this possible role of the socalled backward bifurcation in dengue.
To summarize the above considerations, we list some recommendations, using vectorborne diseases models, to avoid some of our misgivings about the use of models:

(1) Complicated models of vectorborne infections should be used when attempting to make predictions about future events. They are errorprone if not all necessary biological realities are included otherwise their predictions may be unreliable. Many models of vectorborne diseases are designed to estimate parameters that are important to plan control measures. In the example we use in this paper, the estimation of the basic reproduction number from the initial phase of an epidemic is examined and examples of biological realities that should be included are given. Alternatively, we suggested another method to collect data: ideally a uniform region should be chosen with both populations uniformly distributed and only count cases from this region, from the beginning of the outbreak. Because people commute into and out this region, the above method underestimates the basic reproduction number. This error depends on the size of the region that should be carefully estimated. The region cannot be too large because of the propagation of epidemic waves as explained above. On the other hand, it cannot be too small to minimize the effect of people commuting into and out of the region. We plan to publish details of how to choose such a region in another paper. We should stress that most published calculations of the basic reproduction number are not bad approximations of R
_{0}, but should be carefully interpreted to account for error bounds.

(2) Models that have as a goal examination of the mechanism behind certain phenomenon should be as simply as possible. The model of this type examined in this paper deals with the existence of backward bifurcation. The model has, perhaps, too many variables. For instance, differential mortality in humans by dengue is known to be small and in the model studied is perhaps the cause of the socalled backward bifurcation. We conclude that this additional mortality should have not been introduced in the model for dengue. Of course the model may be useful for other diseases.
APPENDIX I. The initial condition
As explained before, most models of vectorborne infections assume (sometimes without saying so) that the variables are population densities and usually homogeneously distributed over a sufficiently large area. This approach has two problems: the first, solved in this Appendix, is that modifications of this model are needed before it can be applied to situations where inhomogeneous populations occur; the second problem is that authors usually assume that the disease is also introduced into the region in a uniformly distributed manner. In fact, the disease is introduced in a small area of the region and the disease propagates through the region as a wave [21]. Let us explain the latter statement in detail.
When the disease is introduced in a virgin area (village, city, etc.), it is introduced in a small geographical area of this region and then propagates as a wave to other places of the region. Since what is usually reported is the number of cases per unit time (epidemiological week, for instance) the method discussed in the ‘Estimating R
_{0} …’ section can only be used if the disease propagates slowly in space. This is the reason why only a few points of the epidemic curve can be used. When not used properly, the method can fail to produce a good measure of the basic reproduction number and in any case the number so derived applies only to the small region where the infection was introduced. Of course if the region is uniform this number is valid for the whole region. In this Appendix we show how to model the propagation of the infection although in this paper this is done in a very sketchy manner; however, this is important as explained above. However, the propagation of an epidemic throughout a region is of central importance in understanding the infection dynamics and in interpreting the available data. It is extremely difficult to include the mobility of humans and vector populations and we suspect that in some cases no general statements can be made. This section, however, is important in highlighting that some general consequences can already be deduced without investigating in detail some points, namely, land use, population density, overlay, etc.
In this Appendix, we compare the outcomes of two epidemics, one of which assumes that the disease is introduced uniformly in the region and the second which assumes the disease is introduced in just a small area and then propagates. Thus, by comparing the two results we can estimate the effect of the fact that epidemics usually start in small regions.
We begin estimating this effect on a onedimensional road of length D. In the case in which we assume that the disease is introduced uniformly, the system of equations (1) can be numerically integrated, and the resulting epidemic compared. To calculate the effect of introducing the disease in just small regions, the system of equations (1) has to be modified to consider the spatial dimension. Let S
_{
H
}(x,t)dx be the number of susceptible humans living between x and x + dx at time t. The other variables are defined similarly.
where
and
In equation (26), a
β
_{
H
}(x,x′) is the number of bites per unit time that I
_{
M
}(x′,t)dx′ infected mosquitoes inflict in S
_{
H
}(x,t)dx susceptible humans (note that infected mosquitoes are in position x′ and susceptible humans are in position x). Similarly, in equation (a
β
_{
M
}(x,x′) is the number of bites per unit time that S
_{
M
}(x′,t)dx′ susceptible mosquitoes inflict in I
_{
H
}(x,t)dx infected humans (susceptible mosquitoes are in position x′ and infected humans are in position x). As a first approximation, we may assume that β
_{
H
}(x,x′) and β
_{
M
}(x,x′) are the same functions. For simplicity, we also assume that β
_{
H
}(x,x′) = β
_{
H
}(x−x′) and β
_{
M
}(x,x′) = β
_{
M
}(x−x′) are functions of the distance x−x′ only.
The standard way to solve system (25) is by choosing one of its variables [say, I
_{
H
}(x,t) and writing an integral equation for it. The result in this case is an integral equation with 16 terms that will be analysed in a future publication.
We can, however, estimate the phenomena we are investigating by selecting a less complicated system [21] that describes a simple SIR model, given by equation (3) of the above reference.
As mentioned by Postnikov & Sokolov [22, pp. 209], equations (3), (4) and (10) of [21] are more general than the equations most commonly seen in the literature.
In a twodimensional region, which, for simplicity, will be considered circular, the system of equations (25) is modified by replacing x with r, the distance from an arbitrary point to the centre, and adding an angular variable θ. Thus, S
_{
H
}(r,θ,t)ds is the number of susceptible humans living in a small area ds = r
dr
dθ around the point (r, θ) at time t. This unfolds similarly for the other variables.
On the other hand, equations (26) and (27) become
and
Returning to the onedimensional case, we obtain, by solving numerically equation (10) of reference [21], the results reproduced in Figure 3a
for the case where the initial condition covers the entire road. The epidemic,
$I(t) = \int_0^L {I(x,t)}dx$
, consists of a huge peak that almost disappears and later returns after a significantly long time. This is not what is observed in the real world.
Fig. 3. Numerical solution for the onedimensional case when: (a) the initial condition covers the entire road; (b) the disease is introduced on one end of the road. Figures not to the same scale.
When the disease is introduced only at one end of the road, Figure 3b
shows the epidemic [number of cases vs. time t, i.e.
$I(t) = \int_0^L I(x,\,t)dx$
]. In this case, we can see that the epidemic rises much slower and lasts for a very long time, being interrupted only if seasonality is considered. This can be understood in the following way: the epidemic wave travels from, for example, the left side of the road to the right side. The epidemic that started because of the initial condition, at the left side, vanishes, but it is still recorded by the cases that are happening on the wave front. It may happen that the disease reappears on the left side before the wave reaches the right side. We then have the impression that the disease reaches steady state.
Returning to the case of the calculation of R
_{0} from the initial phase of the outbreak, we can see that, in nature, the disease is frequently introduced in a small part of the environment and that because of what was stated above, one has to consider only the very beginning of the outbreak. Alternatively, we can modify the way data are collected. Ideally one should choose a uniform region with both populations uniformly distributed and collect cases only from this region from the beginning of the outbreak. Because people commute into and out this region, the above method underestimates the basic reproduction number. This error depends on the size of the region that should be carefully estimated. Other challenges of introducing spatial heterogeneities in epidemic models are qualitatively analysed in [8].
APPENDIX II. Further comments on Garba et al.'s [6] findings
For a more detailed comparison with the paper by Garba et al. [6], we consider the following system of differential equations:
We intend to show that the above system presents a threshold that has a branch with negative values.
In the case of dengue, the diseaseinduced mortality rates in the populations, α
_{
M
} in the vector population and α_{
H
} in the human population, are usually assumed to be zero. In the calculations below, we keep those as different from zero, for completeness.
The stationary solution for human prevalence is given by
where
When A
_{3} = 0, we have
which is equivalent to
the threshold condition for the basic reproduction number is R
_{0} = 1. In this case (A
_{3} = 0), the coefficient A
_{2} becomes
For A
_{3} = 0, one root of I
_{
H
}/N
_{
H
} is zero and the second root of equation (31) is
The expression for X−α
_{
H
} is
and the expression for Y−α
_{
M
} is
Both X−α
_{
H
} and Y−α
_{
M
} are positive, and, as a consequence, the coefficients A
_{1} and A
_{2} are positive. Hence, the second root [equation (36)] of I
_{
H
}/N
_{
H
} is negative. Thus, we prefer not to call the interesting findings by Garba et al. [6] a backward bifurcation. Our interpretation of the phenomena observed by Garba et al. is described in the present paper in the last four paragraphs of the section ‘Backward bifurcations’.
The translation from the notation of this paper to the notation used by Garba et al. [6] is presented in Tables 4 and 5.
Table 4. Translation from the notation of this paper to the notation used by Garba et al. for the variables
Table 5. Translation from the notation of this paper to the notation used by Garba et al. for the parameters
The results obtained by Garba et al. [6] may not be entirely correct because they come from writing λ
_{
H
} [their equation (1)] as
The correct equation should have N
_{
H
} instead of N
_{
M
} in the denominator because, as can be noted in system (1), the total number of bites inflicted by infected mosquitoes is a(I
_{
M
} + η_{
M
}
L
_{
M
}), a fraction S
_{
H
}/N
_{
H
} of which are on susceptible humans, of which a fraction b is actually infective. Therefore, as mentioned previously, equation (4) in Garba et al.'s [6] paper is correct only if b = c.
ACKNOWLEDGEMENTS
This work was partially supported by LIM01 HCFMUSP, Fapesp, CNPq, Dengue Tools under the Seventh Framework Programme of the European Community, grant agreement no. 282589, and MS/FNS (grant no. 27835/2012).
DECLARATION OF INTEREST
None.
REFERENCES
1.
Mazilu, DA, Zamora, G, Mazilu, I.
From complex to simple: interdisciplinary stochastic models. European Journal of Physics
2012; 33: 793–803.
2.
Ross, R.
The Prevention of Malaria, 2nd edn.
London: John Murray, 1911, 772 pp.
3.
Macdonald, G.
The analysis of equilibrium in malaria. Tropical Disesase Bulletin
1952; 49: 813–828.
4.
Coutinho, FAB, et al.
Threshold conditions for a nonautonomous epidemic system describing the population dynamics of dengue. Bulletin of Mathematical Biology
2006; 68: 2263–2282.
5.
Pinho, STR, et al.
Modelling the dynamics of dengue real epidemics. Philosophical Transactions of the Royal Society of London, Series A
2010; 368: 5679–5693.
6.
Garba, SM, Gumel, AB, Abu Bakar, MR.
Backward bifurcations in dengue transmission dynamics. Mathematical Biosciences
2008; 215: 11–25.
7.
Anderson, RM, May, RM.
Infectious Diseases of Humans: Dynamics and Control. Oxford: Oxford University Press, 1991.
8.
Riley, S, et al.
Five challenges for spatial epidemic models. Epidemics (in press).
9.
Lopez, LF, et al.
Threshold conditions for infection persistence in complex hostvectors interactions. Comptes Rendus Biologies Académie des Sciences Paris
2002; 325: 1073–1084.
10.
Amaku, M, et al.
A comparative analysis of the relative efficacy of vectorcontrol strategies against dengue fever. Bulletin of Mathematical Biology
2014; 76: 697–717.
11.
May, RM, Anderson, RM.
The transmission dynamics of human immunodeficiency virus (HIV). Philosophical Transactions of the Royal Society of London, Series B: Biological Sciences
1988; 321: 565–607.
12.
Massad, E, et al.
The risk of yellow fever in a dengueinfested area. Transactions of the Royal Society of Tropical Medicine and Hygiene
2001; 95: 370–374.
13.
Favier, C, et al.
Early determination of the reproductive number for vectorborne diseases: the case of dengue in Brazil. Tropical Medicine and International Health
2006; 11: 332–340.
14.
Marques, CA, Forattini, OP, Massad, E.
The basic reproduction number for dengue fever in São Paulo state, Brazil: 1990–1991 epidemics. Transactions of the Royal Society of Tropical Medicine and Hygiene
1994; 88: 58–59.
15.
Massad, E, et al.
Estimation of R
_{0} from the initial phase of an outbreak of a vectorborne infection. Tropical Medicine and International Health
2010; 15: 120–126.
16.
Burattini, MN, Coutinho, FAB, Massad, E.
A hypothesis for explaining single outbreaks (like the Black Death in European cities) of vectorborne infections. Medical Hypotheses
2009; 73: 110–114.
17.
Massad, E, et al.
Modeling the risk of malaria for travelers to areas with stable malaria transmission. Malaria Journal
2009; 8: 296.
18.
Massad, E, Coutinho, FAB.
The cost of dengue control. Lancet
2011; 377: 1630–1631.
19.
Amaku, M, et al.
Maximum equilibrium prevalence of mosquitoborne microparasite infections in humans. Computational and Mathematical Methods in Medicine
2013; 2013: Article ID 659038.
20.
Coudeville, L, Garnett, GP.
Transmission dynamics of the four dengue serotypes in Southern Vietnam and the potential impact of vaccination. PLoS ONE
2012; 7: e51244.
21.
Lopez, LF, et al.
Modeling the spread of infections when the contact rate among individuals is short ranged: propagation of epidemic waves. Mathematical and Computer Modelling
1999; 29: 55–69.
22.
Postnikov, EB, Sokolov, IM.
Continuum description of a contact infection spread in a SIR model. Mathematical Biosciences
2007; 208: 205–215.