SEASONAL FORCING AND EXPONENTIAL THRESHOLD INCIDENCE IN CHOLERA DYNAMICS

We propose a seasonal forcing iSIR (indirectly transmitted SIR) model with a modified incidence function, due to the fact that the seasonal fluctuations can be the main culprit for cholera outbreaks. For this nonautonomous system, we provide a sufficient condition for the persistence and the existence of a periodic solution. Furthermore, we provide a sufficient condition for the global stability of the periodic solution. Finally, we present some simulation examples for both autonomous and nonautonomous systems. Simulation results exhibit dynamical complexities, including the bistability of the autonomous system, an unexpected outbreak of cholera for the nonautonomous system, and possible outcomes induced by sudden weather events. Comparatively the nonautonomous system is more realistic in describing the indirect transmission of cholera. Our study reveals that the relative difference between the value of immunological threshold and the peak value of bacterial biomass is critical in determining the dynamical behaviors of the system.

1. Introduction.Cholera is a severe intestinal disease caused by ingesting water contaminated with the bacterium Vibrio cholerae.The main symptom of cholera infection is profuse watery diarrhea that rapidly causes dehydration and can lead to death within days if not promptly treated.Although listed as one of the oldest known diseases, cholera remains a serious public health burden in developing countries where poverty and poor sanitation and hygiene are prevalent.Major cholera outbreaks in recent years include those in Zimbabwe in 2008, Vietnam in 2009, Nigeria in 2010, Ghana in 2011, as well as the one in Haiti during 2010-2012 which is regarded as one of the largest cholera epidemics in modern history, with more than 530,000 reported cases and over 7,000 deaths [24].
The remainder of this paper is organized as follows.In Section 2, we introduce the new exponential saturation incidence and incorporate periodicity into the iSIR model.In Section 3, we provide the forward invariance, the uniform persistence, the existence and stability (local and global) of periodic solutions, for the nonautonomous system.In Section 4, we conduct numerical simulations to validate the analytical results and to test the impact of sudden weather events on cholera outbreaks.Finally, we summarize and discuss our conclusions in Section 5.
2. Model derivation.The autonomous iSIR model was introduced in [10]: where the saturation function α(B) takes the form as follows: , where n is a positive integer.
( Although as n becomes larger, the curve will become smoother, these curves are not sufficiently smooth, posing challenges in mathematical analysis.we thus introduce the following novel exponential threshold incidence function which is infinitely smooth: , where n is a positive number.(3) Obviously this new function satisfies the half saturation condition (when B −c = H, we have α(B) = a/2).
In addition, we define the value of the i-th derivative of α(•) at c as 0 for i = 1, 2, • • • .It is easy to show that for every positive number n, the function α(•) is infinitely smooth.The advantage of such defined functions is reflected not only in theoretical analysis but also in numerical simulations.We will highlight this point throughout the paper.
To our knowledge, such an exponential threshold incidence from is first of its kind.It is clear from Figures 1 and 2 that these exponential functions match corresponding Holling type functions very well, yet with better mathematical properties that enable us to deepen our understanding of the rich dynamics of the cholera model.In addition, these exponential saturation functions can be potentially applied to the modeling of many other infectious diseases where the minimal infection doses (MID) are relevant.For these reasons, we carefully analyze the exponential incidence incorporated in system (1), under a relatively simple, autonomous setting; the details are provided in Appendices.The global stability analysis repeats almost the same process in Kong et al. [11], thus we do not present global results here.The results of the autonomous system build a solid base for our mathematical analysis of the more complex, nonautonomous system in Section 3. The parameters in system (1) are described in Table 1, taken from Jensen et al. [9].Seasonal forcing is believed to be an important culprit for cholera outbreaks due to the indirect transmission from the reservoir that heavily depends on weather events.Here we explicitly incorporate seasonality into the bacterial intrinsic growth rate and carrying capacity in our iSIR model, and obtain where α(B) takes the exponential threshold form of (3), and S, I, R satisfy S + I + R = N .We assume r(t) and K(t) are all positive continuous periodic functions, i.e., there exist ω, r m , r M , K m , and K M such that r(t + ω) = r(t), K(t + ω) = K(t), and For simplicity, we set s = S/N, i = I/N , and drop equation(4c).We rewrite system (4) in an equivalent form as 3. Mathematical analysis.

3.2.
Persistence.We proceed to show the uniform persistence of the system (5).
To that end, we first establish the following result.
Since the equation possesses a globally stable solution u = µ a+µ , comparing the two equations ( 11) and ( 12), we have lim Set ε 1 = µ a+µ , and we obtain (7).Picking µ 1 such that 0 < µ 1 < ε 1 , then there exists a T 1 > 0 such that By (5b), we have di dt It follows Integrating both sides over [T 1 , T ], we get Hence we have Let T = T 1 + nL + ∆T , where 0 ≤ ∆T < L and n is a non-negative integer.Considering (10), we have Taking due to ( 16) and ( 17), we have lim Now by (8), there exists a T 2 > 0 such that thus by (4c) we have It has a solution v = N δε2 2µ which is globally stable.Based on (19) and (20), by the comparison theorem [23], we have which implies (9).
Remark 3.3.Obviously, ( 7),( 8) and ( 9) imply the persistence of (4).The hypothesis of this theorem has a reasonable background in reality.For example, if for every year there is a period, even if very short, within which the bacterial density is higher than the threshold value c, then the infection will spread and the disease will persist.
The hypothesis can also be given in another form.Consider the equation where r(t), K(t) are positive continuous ω-periodic functions.It can be shown that this equation possesses a unique periodic solution B, which is globally asymptotically stable.Specifically, we have the following result.
According to (5c) and (21), by the comparison theorem [23] we know that for any solution of system (5), (s(t), i(t), B(t)) with positive initial values, there must exist a T > 0 such that B(t) > B(t), for t > T.
3.3.Stability.Generally speaking, stability analysis of non-autonomous dynamical systems is challenging.In what follows, we attempt to study the trivial and non-trivial periodic solutions of our model (5) under some stronger conditions.
We first present a known result below that characterizes the stability relationship between the trivial solution of a nonlinear system and that of the corresponding linear system [7] [21]. where (a) If X = 0 is uniformly asymptotically stable for then X = 0 is uniformly asymptotically stable for system (23).
We apply this lemma to system (5), and discuss two cases separately.
Case 2. Consider i = 0, B = 0.It is stated for equation ( 21) that there exists an unique solution B(t).Assume max t∈[0,ω] B(t) ≤ c, then we have an unique solution for system (5)

In this case we have
Hence the stability depends on the sign of the third Floquet exponent Now we will show that λ 3 < 0. By the definition of B(t) we simply have Obviously B(t) = 0 for t ∈ [0, ω].Dividing both sides of the above equation by B(t), we have Integrating both sides of the above equation from 0 to ω, then we have due to the fact that B(t) is ω-periodic.Hence
Now we consider the stability of a general positive ω-periodic solution of system (5).We prove the following theorem first.Theorem 3.10.Suppose that c < max t∈[0,ω] B(t).Then there exist an ε > 0 and a constant T > 0 such that, for every solution (s(t), i(t), B(t)) of system ( 5) with positive initial condition, we have Proof.By the assumption c < max t∈[0,ω] B(t) and theorem (3.5), there exist an Consider equation (5c) and let this inequality hold provided 2 . Hence It follows Taking we see that while Ecologically, this theorem implies that, as some individuals infected, the bacterial density could surpass the natural carrying capacity, and it may result in more susceptible population get infected even their threshold value c is greater than max t∈[0,ω] B(t).It could be a vicious circle if the patients could not be cured in time.And it may cause environmental deterioration in some poor area where medical facilities are insufficient.Theorem 3.11.Suppose that c < max t∈[0,ω] B(t).Then the ω-periodic solution (s * (t), i * (t), B * (t)) guaranteed by theorem (3.6) is globally asymptotically stable if Proof.Let (s(t), i(t), B(t)) be any solution of system (5) with nonnegative initial condition.Define the following Lyapunov function We calculate the right-hand derivative of L(t) along with system (5) by using the fact that |x| ′ = sign(x)x ′ : In any time t, the sign of i(t) and s(t) should be included in the following three cases, and we shall consider three cases respectively.
Case 1.By theorem (3.10), there exist an ε > 0 and a T > 0 such that Case 2. We have four subcases to discuss.
(33) By (5c) we have Hence, by lemma (3.10), for ∆B > 0 we have and for ∆B < 0 we have where T is given by lemma (3.10 we have lim t→∞ ∆B = 0. Thus in this case, (3.11) still holds.Now, after the discussion of the above three cases, we have Corollary 3.12.Suppose that c < max t∈[0,ω] B(t), then system (5) admits an unique and globally asymptotically stable ω− periodic solution if Remark 3.13.When α(•) takes the form of Holling type described in equation ( 2), it reaches its maximum at B = c the corresponding condition for the global stability of the periodic solution is Remark 3.14.From the proof of the theorem, it is clear to see that the condition (26) is equivalent to max α ′ (B) < rmδ 2N ξ .Obviously, α ′ (B) describes the rate that the incidence grows when B increases.If this growth is not so fast; i.e., bounded above by rmδ 2N ξ , then a unique and globally asymptotically stable periodic solution exists.

4.1.
Simulations for the autonomous system (40).As space is limited, we omit the case when c is much less than or is much greater than the carrying capacity of bacteria, in which the corresponding system will be stable or be extinct respectively.We only present the reader a special case, in which the system possesses a bistability, when the threshold value c is close to the carrying capacity of bacteria.
The reader may wonder whether X 2 is essentially (1, 0, K).Here we emphasize that X 2 is not (1, 0, K).Firstly, by simulation we see that the final state X 2 is locally independent of initial conditions, which means that X 2 is locally stable.This example is a special case shown in figure 15, where X 1 and X 2 are corresponding to E 3 and E 1 , both of which are locally stable.Secondly, the result is in accordance with the computational results given by (41).Thirdly, theoretically system should be permanent according to remark 3.15.
This example shows that when c is slightly greater than K, with a high density of bacteria in initial condition, solutions may approach a final state X = (1.3185× 10 −3 , 9.9768 × 10 −4 , 1.2805 × 10 6 ) which possesses a bacterial density higher than K.This is also an example to show that (1, 0, K) is not globally stable in case of c > K.

Simulations for the nonautonomous system.
In what follows, we rewrite the time-dependent bacterial intrinsic growth rate by r(t) and the carrying capacity by K(t) to avoid confusion of notations.To represent seasonal fluctuation, we consider r(t) and K(t) in a simple form as   with a period ω = 365 days.Here the constants K and r (i.e., the base values of the intrinsic growth rate and carrying capacity) take values listed in Table 1.When the threshold of immunity is significantly higher then the maximum of carrying capacity of bacteria, solutions of system (5) might tend to a disease free periodic solution (1, 0, B).See figure 6, where we take n = 1, a = 0.1, δ = 0.1, r = 10, ξ = 20, µ = 10 −4 , c = 1.5 × 10 6 , N = 10 6 , K = 10 6 , H = 10 7 .
Example 4.5.When c is close to maximum of K(t), slightly less than maximum of B(t), H is relatively low, ξ is relatively high, such as K = 1 × 10 6 , N = 1 × 10 6 , r = 0.25, δ = 0.1, ξ = 90, µ = 1 × 10 −4 , H = 5 × 10 5 , a = 0.1, n = 2, and seasonal fluctuation is mild, say r(t) = (1 + 1/10 * sin(t))r; K(t) = (1 + 1/10 * sin(t))K , system (5) may have different locally stable state, depending on initial conditions.It's easy to verify that condition (26) is not satisfied.Presently we have detected three kind of different final states.Taking the initial values x 0 (0) = (0.2 × 10 6 , 0.00001 × 10 6 , 1.2 × 10 6 )( Here B(0) is higher than the peak value of K(t)), system will approach to a stable state X = (s 0 (t), i 0 (t), B 0 (t)), which is a periodic solution fluctuating around (9.99 × 10 −4 , 9.98 × 10 −4 , 1.27 × 10 6 ), (See figure 7, right column).Taking x 1 (0) = (0.2 × 10 6 , 0.00001 × 10 6 , 0.8 × 10 6 ) (Here B(0) is less than the peak value of K(t)) system will tend to the second kind of final state X 1 (t), which is mainly very close to the disease free solution (0, 0, B(t)) , and occasionally accompanied by short time eruptions of infection(See figure7, left . Two final states of system (5) depending on different initial values column).Taking x 2 (0) = (0.2×10 6 , 0.01×10 6 , 0.8×10 6 ) system will approach to the third kind of final state X 2 (t).Noticing that while initial conditions x 1 (0) and x 2 (0) has the same fraction of susceptible population and bacterial density, X 2 (0) has a higher fraction of infective population.After having experienced a large scale of outbreak of infection, X 2 (0) leads the system to a different final state X 2 (t), which has a higher bacterial density and a sharp declining of the fraction of susceptible population than X 1 (t).Obviously X 2 (t) marks a worsen environment(due to the higher bacterial density).In the simulation from the enlarged figure of infective population, we can detect a lot of short time eruptions with different scales.These could be very small, say lower than 10 −10 (See the left column of figure 7), and it could hardly be detected.In reality this might not represent an infected individual, but still it could be treated as a very slight infection of individuals.Practically when a person is infected in such a slight degree, he may be uninfected.
Example 4.6.In this example we consider two different situations when xi is zero or nonzero.Taking c = 1.085 × 10 5 , which is slightly less than the peak value of B(t), the bacterial density in disease free system, and values of other parameters are as in Example 4.5.Taking initial value x(0) = (7.5 × 10 6 , 0.05 × 10 6 , 1 × 10 6 ).We see that when ξ = 0, there is an outbreak of epidemic every year which is synchronous with the peak value of B(t), and after this short period, the fraction of infective population is very low(nearly zero), and the average fraction of susceptible population is between 1 percent and 2 percent, namely there are about 1 to 2 percent of people remain uninfected.When ξ = 90, the fractions of susceptible and infective population are both nearly 0.001, and this shows that there are only 0.1 percent of the total population are not infected.Strictly speaking, they are totally two different models when ξ = 0 and ξ > 0. These two models are fit for different environment.In those places, when population density is low, and public sanitation is relatively good, and people have less contact, the model for ξ = 0 is suitable.In places when population density is high, public sanitation is bad, and people may have to contact frequently, the model for ξ > 0 is more suitable.In short, when ξ > 0, infection between people are considered.See figure 9.
Example 4.7.In this example we assume the biological system encounters a sudden event, say an abnormal reproduction of bacteria due to the major natural disasters.Take K = 1 × 10 6 , r = 0.25, δ = 0.1, ξ = 90, µ = 1 × 10 −4 , H = 8 × 10 5 , a = 0.1, n = 2, x(0) = (0.5 × 10 6 , 0.001 × 10 6 , 0.8 × 10 6 ).For c = 1.2 × 10 6 , which is significantly higher than the peak value of K(t), normally system will tend to a disease free state.However, when system encounters a sudden event, it may have a different behavior depending on system parameters.These simulations are in connection with total populations.In the first case, taking N = 1 × 10 6 , system experienced a severe infection, caused a sharp declining of the fraction of susceptible population.After that, the state of the system gradually recovers to the primary state (See the left column of the figure 10).In the second case, taking N = 1 × 10 7 , we find that the sudden event brings a long term influence on the system behavior.After a severe infection, system tend to a new stable state.The fraction of susceptible population has a vast slash, and the fraction of infection population has a substantial increase, and the bacterial density rises to a higher level (See the right column of figure 10).This new state marks a worsen environment caused by the sudden event.
Remark 4.8.From these examples, we observe that the threshold value c is critical in determining the dynamical behaviors.In section 3, we have obtained a sufficient condition for the stability of system (5) when c is less than the peak value of B. However when c is significantly less than the peak value of B, the condition is not important for the stability.This can be explained as follow.From the process of the above proof, we see that when c < max t∈[0,ω] B(t)the condition for the stability of the periodic solution suffices to ν > 0. When c is much less then max t∈[0,ω] B(t) or the carrying capacity of bacteria, α(B) could be much greater than 0, which implies a greater A in (10), and a greater ε 2 in (18), which in turn implies a greater i by (8).Thus by the proof of lemma (3.10), it implies a greater B. When B is much greater then B, α ′ (B) could be much less than α ′ ( B), and thus condition ν > 0 could be much easier to be satisfied, which means the stability of the periodic solution.When c is close to max t∈[0,ω] B(t) or the carrying capacity of bacteria, condition (26) is crucial to determine whether the periodic solution is stable or not.
Comparatively, the nonautonomous model is more realistic than the autonomous one.As it is shown in Example 4.1, 4.3-4.7,nonautonomous model exhibits more complexities especially when c is close to the peak value of K(t).Example 4.6 presents an interesting phenomena that an eruption of epidemic may occur unexpectedly.This is congruent to what we encounter in reality.

5.
Discussion.Cholera remains epidemic and endemic regionally in the world, especially in the locations lacking adequate sanitation and water infrastructure.With the continuing outbreaks, mathematical modeling plays an important role in deciphering its dynamics and providing suggestions for governments and health to take effective actions.Cholera is an indirectly transmitted infectious disease, and recent modeling efforts have been made in this direction, such as the addition of an environmental pool and the incorporation of an immunological threshold for infection [4,10,11].Subject to the annual variation of temperature, nutrients, rainfall, etc. in the reservoir, seasonality is believed to be pivotal in determining cholera epidemics and endemism [6,17].Seasonal drivers have been recently incorporated in cholera transmission models [5,15,16].In this paper, we incorporate the seasonal factor explicitly in bacterial growth term and discuss its impact on cholera dynamics.In addition, we originally introduce a new threshold incidence function which is infinitely smooth.Our proposed incidence function biologically mimics the standard Holling type functions, and the smoothness of the new incidence form allows us to perform rigorous and deeper mathematical analysis for the complicated dynamics of an indirectly transmitted infectious disease.
It is mathematically challenging to analyze a high-dimensional nonautonomous system.In the paper, we provide the forward invariance, the uniform persistence, the existence of periodic solutions, and their stability (local and global), for our nonautonomous model.Forward invariance and persistence illustrate the rationality of the proposed seasonal forcing model.The immunological threshold is the key parameter in determining persistence, periodic solutions and their stability, that is, the mean infection threshold of susceptible individuals in an epidemic region is critical in determining cholera outbreaks and their severity.Hence in order to prevent an outbreak of cholera, we have either to raise the immunity of susceptible individuals or to improve the local environments, bring down the carrying capacity of bacteria which determines the periodic solution B in our model.To prove the persistence and the stability of periodic solutions for a nonautonomous highdimensional dynamical system is of considerable challenge.Our proofs are new and may be applicable to other models of this type.Our numerical simulations validate our mathematical results.
Further steps to take with the seasonal forcing model would be to refine the conditions on the global stability results of periodic solutions.The imposed conditions are sufficient but not necessary for global stability.To verify and calibrate the nonautonomous model with seasonality, we should fit our theoretical outputs to empirical data from regional cholera outbreaks such as the reported cases in Bangladesh and Haiti.The seasonality is clearly required to fit a multi-year data set.In the modeling perspective, it is important to mechanistically incorporate temperature, nutrients, and sudden events such as hurricane and eddy flow, for the environmental reservoir.
Appendix A. Equilibria of the autonomous system.We start our analysis by examining the equilibria of the autonomous model.Since S + I + R = N where the total population N is a constant, we will drop R from the system.Meanwhile, introducing s = S/N, i = I/N , we obtain an equivalent form to the system (1): where the saturation function α(B) takes the exponential form of (3).An equilibrium (s * , i * , B * ) of system (40) satisfies System (40) has possibly two, three or four equilibria.The exact number of equilibrium points not only depends on the threshold value c, but also on the values of parameters of K, r, N, ξ, µ, δ, etc.However, in any case E 0 (1, 0, 0) is an equilibrium of system (40).Combining (41b) and (41c) we form an equation where ν = N Kξµ r(µ + δ) . Denote Thus equation ( 42) is equivalent to . Now we will work out the number of equilibrium points of equation ( 42) through the functions f (B) and g(B).Though f (B) is quite simple, we still need some analysis to get a few characteristics of g(B).
Denote h n (B) = nH n ln 2 (B−c) n .By (3), we have The third property of g(B) is a direct result of the following proposition.From the discussion above we can see that the shape of the curve of function g(B) is quite similar to that of α(B).What is important is that the shape is directly related to the number of equilibria of the system.
Let g c (B) = g(B + c).Given a number M x that is large enough, when c increases from 0 to M x , the curve g c (B) has a parallel translation from passing the origin, to crossing the curve f (B), and then to the positive direction of B axis.The influence of the threshold value c on the equilibria of system (40) can be seen in this process of movement.
Denote the point corresponding to B in the curve of function g(B) as M .If the position of M corresponding to c = 0 is on the left of the curve of function g(B), then as c increases from 0 to large enough, M has a parallel movement from left to right along the positive B axis and it will pass through g(B).There must be a point c such that M is on the curve of function f (B), and we denote the corresponding abscissa as BM .At BM we possibly have Case i).g ′ ( BM ) ≤ f ′ ( BM ), and Case ii).g ′ ( BM ) > f ′ ( BM ).
Both cases are possible depending on the parameter values we choose.
Proof.we only need to prove that apart from E 0 , Equation (42) has only one equilibrium E 1 .
Appendix B. Stability of the autonomous system.We now analyze the stability of the equilibria for the autonomous model.The Jacobian of system ( 40) is In particular, the Jacobian corresponding to E 0 (1, 0, 0) is Since J 0 has three eigenvalues, −µ, −µ − δ and r, two of which are negative and the third one is positive, we have Theorem B.1.The disease-free equilibrium E 0 (1, 0, 0) for system (40) is an unstable saddle-node.
In case of c > K, system (40) possesses another equilibrium E 1 (1, 0, K), and the corresponding Jacobian is The characteristic polynomial of J 1 is = (λ + µ)(λ + µ + δ)(λ + r) Since all eigenvalues of J 1 are negative, we have Theorem B.2.In case of c > K, system (40) possesses an equilibrium E 1 (1, 0, K), which is locally asymptotically stable.To proceed, we apply the following standard result from the Routh-Hurwitz stability criterion: Lemma B.3.The necessary and sufficient condition for the polynomial (52) to be stable is a 0 > 0, a 1 > 0, a 2 > 0, a 3 > 0 and a 1 a 2 > a 0 a 3 .
The condition a 0 > 0 is trivial.Since B * > K, we obtain a 1 > 0.
According to (42), we have B * (B * − K) < ν. Hence, and we have If a 3 ≤ 0, according to Lemma (B.3), the corresponding equilibrium is unstable.The condition a 3 ≤ 0 is equivalent to F ′ (B * ) ≥ 0, and also equivalent to f ′ (B * ) ≤ g ′ (B * ), which means the slope of the curve of g(B * ) is equal or steeper than that of f (B * ).This result provides an intuitive way for the judgment of the instability of the equilibrium E(s * , i * , B * ), which we refer to as a geometrical method.We state the above results as follow.Geometrically, this theorem suggests that if the slope of the tangent line of f (B) at B * is much larger than that of g(B), then the equilibrium E(s * , i * , B * ) is locally asymptotically stable.

Figure 4 .
Figure 4.An example when system (40) does not approach (1, 0, K) in the case that c is slightly greater than K

Figure 5 .
Figure 5. Populations of system (5) approach a periodic solution

Figure 11 .
Figure 11.Curves of u(B) and v(B) have an unique intersection B

Figure 12 .a µ 2 ,1+log a µ 2 ,
Figure 12.Curves of function f and g with changing threshold values c.

Table 1 .
[9]ameter values from Jensen et al.[9] |s−s * |. (t).It is similar to case 2.2.Now we assume condition (26) holds, then ν > 0. Taking ζ = min{µ, r m , ν}, and by ( * (t) and B(t) > B * 7, each of which is in the range of corresponding parameter values listed in table 1. c is significantly less than the peak value of K(t) (1.1 × 10 6 ).It's easy to verify that condition (26) is satisfied.Simulation results are shown in figure5.We clearly observe a stable periodic solution.In particular, the peak values of the susceptible and infective individuals occur periodically at times when the bacterial concentration is rapidly ascending.See figure5.
Table 1 and it means B − < B. Hence we have B ∈ (B − , B 0 ).H