Impacts of Seasonal and Spatial Variations on the Transmission of Typhoid Fever

and


Introduction
Typhoid fever is an infection caused by the bacteria Salmonella Typhi (S. Typhi), which is usually spread by ingesting contaminated food/water. According to the World Health Organization, typhoid is still endemic in several developing countries. The infection routes of Typhoid include both direct (i.e.human-to-human) transmission and indirect (i.e. environment-to-human) transmission, which is associated with the ingestion of contaminated food/water. It is worth pointing out that asymptomatically infected individuals also play an important role in the transmission of typhoid since they are not experiencing any symptoms but they are able to shed S. Typhi into the environment for many years, thereby sustaining transmission [15,16]. Those observations motivate the authors in [16] developed a system of ordinary differential equations to model the spread of typhoid, where the factors of limited treatment resources on the spread of typhoid was further included.
Our aim of this paper is to incorporate spatial and temporal effects into the model proposed in [16]. During the rainy season, a large increase of flooding may occur by the rainfall and the water resources is contaminated by the excreta with pathogenic bacteria, resulting in a higher incidence of typhoid fever [12,13]. It was evident that the supply of clean and safe drinking water plays an important role in controlling typhoid infection, and people living near water bodies (e.g. rivers) may have a higher rate of typhoid infection [4]. Thus, the effect of spatial homogeneity should be included. To develop spatially explicit models, there are two common approaches. The first one is a continuum approach using an advection-dispersion-reaction system to describe the transport and interaction of population in a (bounded) habitat. Here we adopt the second approach, namely, we can divide the habitat into several areas, and the population gradient between areas is simply described by the migration of population. To make mathematical analysis more tractable, we will focus on the case that the habitat is divided into two areas, which is also referred to as "two-patch model". On the other hand, it is confirmed that temperature changes have significant impacts on the enteric diseases [9,10], and the increase in rainfall and temperature lead to more typhoid fever cases in the study area [4]. Thus, it is natural to further explore the seasonal variations in temperature and rainfall on the transmission of typhoid. Based on those aforementioned facts, we shall propose and analyze a time-periodic system in a two-patch environment which is modified from the one in [16].
The population of human at time t in path i is divided into five classes: susceptible individuals (S i (t)), infected individuals (I i (t) ), carrier individuals (C i (t)), individuals under treatment (Q i (t) ) and recovered individuals (R i (t)). Besides, B i (t) stands for the density of bacteria in the environment in patch i at time t. Then the extended version of the model in [16] takes the following form: (1) The constant Λ i represents the recruitment of susceptible population, and µ i represents the natural death rate for the general population in the environment of patch i. Susceptible people are infected either through human-tohuman transmission at the rate Si+Ii+Ci+Qi+Ri or through the environmental bacteria from contaminated drinking water/food at the rate . The parameter β Ci (t) is the so-called typhoid transmission rate for susceptible individuals and infected/carrier individuals, and η i is used to measure the relative infectiousness of carriers C i compared to infected individuals I i . We will assume 0 < η i < 1 when the carriers C i have less infectious ability than infected individuals I i . Otherwise, we will assume that η i ≥ 1. The parameter β Bi (t) is the per capita contact rate between susceptible individuals and the contaminated environment, and K Bi is the saturation constant. Infected individuals progress to the carrier class at the rate σ i ; the naturally recovery rate for infected (resp. carrier) individuals is denoted by Ii (resp. Ci ); the mortality rate due to disease for infected (resp. carrier) individuals is denoted by δ Ii (resp. δ Ci ); the parameter α Ii (t) (resp. α Ci (t)) represents the shedding rate of bacteria by infected (resp. carrier) individuals. The recruitment into treatment class is denoted by θ i I i . The mortality rate due to illness in patients under treatment is δ Qi , and its recovery rate is γ i . The recovered individuals will only be temporarily immune to typhoid, leading to the individual being susceptible again at the rate ρ i . The generation rate of bacteria is expressed in terms of g i B i , where g i is a constant; the production rates of bacteria from infected persons and carriers are denoted by α Ii (t) and α Ci (t), respectively. We further assume that the bacteria in the environment becomes non-infectious at a rate µ Bi . Here, for w = S, I, C, R, B, m w 21 represents the immigration rate of population w from patch 2 to patch 1, while m w 12 represents the immigration rate of population w from patch 1 to patch 2. We also point out that the class Q i , i = 1, 2, represents the typhoid patients who are detected and quarantined symptomatic and chronic enteric carriers. Thus, Q i , i = 1, 2 is supposed to be on treatment, and those terms −m Q 12 Q 1 +m Q 21 Q 2 and m Q 12 Q 1 −m Q 21 Q 2 are ignored in system (1), due to the fact that the population Q i cannot move in the environment. This makes our mathematical analysis more difficult and challenging.
In the whole paper, we always assume that β Ci (t), β Bi (t), α Ii (t), and α Ci (t) are ω-periodic functions, for i = 1, 2, and which coincides with the parameters given in the table on page 664 of [16]. The organization of the rest of this paper is as follows. The well-posedness of our proposed model and the basic reproduction number, R 0 , are discussed in the next section. In Section 3, we show that the global dynamics of our proposed model can be determined in terms of the basic reproduction number, R 0 . Numerical simulations and biological interpretations are presented in Section 4 and Section 5, respectively.

The Reproduction Number
We first consider the following system Then we can verify that system (3) is a cooperative system (see, e.g., [19]) and (S * 1 , S * 2 ) is the only positive steady state. Thus, we have the following result concerning with the global stability of (S * 1 , S * 2 ) (see, e.g., [8]). Lemma 2.1. The positive steady state (S * 1 , S * 2 ) is globally attractive in R 2 + for system (3). That is, we have We further have the following result: Lemma 2.2. System (1) has a unique and bounded solution with the initial value in R 12 + , which is positively invariant. Moreover, system (1) has a connected global attractor on R 12 + in the sense that it attracts all positive orbits in R 12 + . Proof We show that system (1) admits a unique noncontinuable solution and the solutions to (1) remain nonnegative if they are non-negative initially. In view of [19,Theorem 5.2.1], we can show that for any x 0 ∈ R 12 + , system (1) has a unique local solution u(t, x 0 ) ∈ R 12 + with u(0, x 0 ) = x 0 . This shows the positive invariance of R 12 + for system (1). Next, we show that the solutions u(t, x 0 ) for system (1) are eventually bounded. For this end, we let Then we substitute N (t) = N 1 (t) + N 2 (t) into system (1), and it follows that By the positivity of solutions, we see that where µ = min{µ 1 , µ 2 }. Thus, lim sup t→∞ N (t) ≤ Λ1+Λ2 µ . This means that N (t) is ultimately bounded. Since solutions of (1) are nonnegative, it follows from (5) that S i (t), I i (t), C i (t), Q i (t), and R i (t) are ultimately bounded, for i = 1, 2. Then there exists a τ 0 > 0 and χ i > 0 such that This inequality together with the sixth and twelfth equations of (1) imply that In view of the assumption (2), we see that µ 01 := µ B1 − g 1 > 0 and µ 02 := µ B2 − g 2 > 0. , which reveals that B 1 (t) and B 2 (t) are also eventually bounded. Therefore, we have established the global existence for the solutions of system (1), and the rest of the result follows [7, Theorem 3.4.8].
The disease-free state of system (1), E 0 , takes the following form where (S * 1 , S * 2 ) is given in (4). We linearize system (1) at the disease-free state E 0 and we arrive at the following linear system Note that system (7) is cooperative (see, e.g., [19]) and irreducible (see a simple test on page 256 of [20]). From system (7), we define and where Assume Y (t, s), t ≥ s, is the evolution operator of the linear ω-periodic system and I stands for the identity matrix with size 6. Then, for each s ∈ R, Assume that C ω represents the ordered Banach space of all ω-periodic functions from R to R 6 , which is equipped with the maximum norm · . The associated positive cone C + ω is defined by Suppose that φ(s), ω-periodic in s, is the initial distribution of infectious individuals. Then the next infection operator Then the basic reproduction number is given by R 0 := r(L), the spectral radius of L.

Global Dynamics
We first consider the linear ordinary differential system where A(t) is a continuous, cooperative, irreducible, and ωperiodic k × k matrix function. Assume that Φ A(·) (t) is the monodromy matrix of (10) and r(Φ A(·) (ω)) is the spectral radius of Φ A(·) (ω). In view of [1, Lemma 2] (see also [6, Theorem 1.1]) and the Perron-Frobenius theorem [19], we see that We further have the following results: . Then there exists a positive, ω-periodic function v(t) such that e λt v(t) is a solution of (10).
For further discussions, we need the following property.
is a solution of the system (1) with initial value and Then Proof In view of the first equation in system (1), it follows that Then the matrix J(t) is cooperative (see, e.g., [19]) and irreducible (see a simple test on page 256 of [20]), where we have used the fact that S i (t) > 0, ∀ t > 0, i = 1, 2. By (11) and the irreducibility of the cooperative matrix J(t), we can use a generalized version of [19, Theorem 4.1.1] to show that In view of the equations of Q 1 , Q 2 , R 1 and R 2 in system (1), together with (12), we can further show that Q i (t) > 0, R i (t) > 0, ∀ t > 0, i = 1, 2. We complete the proof.
Let X = R 12 + . Define a family of maps {Ψ(t)} t≥0 from X to X by where u(t, x 0 ) is the unique solution of system (1) with u(0, x 0 ) = x 0 (see Lemma 2.2 and the proof therein). Suppose P : X → X is the Poincaré map associated with system (1), that is,  Proof Assume R 0 > 1. Then r(Φ F(·)−V(·) (ω)) > 1 ( see Lemma 2.4). Thus, we find a small ξ 0 > 0 such that Since solutions are continuous with respect to the initial values, we can find a ς 0 > 0 such that for any Suppose that the above claim is not true. Then we have for some x 0 ∈ X 0 with x 0 − E 0 ≤ ς 0 . Without loss of generality, we assume that It follows that Given any t ≥ 0, we rewrite t = mω + t , where t ∈ [0, ω), and m is the largest integer less than or equal to t ω . Therefore, it follows that Note that (13), for all t ≥ 0, we have From the equation of I 1 in (1), we see that From the equation of I 2 in (1), we can use the same arguments to show that Then, for t ≥ 0, we further have the following inequalities In view of the fact x 0 ∈ X 0 and Lemma 3.2, it follows that Thus, we find a fixedt 0 > 0 such that In view of Lemma 3.1, we may find a positive, ω-periodic function J(t) andJ(t) :=deλ (t−t0) J(t) is a solution of
In view of system (21), we see that lim t→∞ Q i (t) = 0, i = 1, 2. Then (R 1 (t), R 2 (t)) in (1) is asymptotic to system (19), and hence, lim t→∞ (R 1 (t), R 2 (t)) = (0, 0). Thus, (S 1 (t), S 2 (t)) in system (1) is asymptotic to system (3). By Lemma 2.1, it follows that lim t→∞ (S 1 (t), S 2 (t)) = (S * 1 , S * 2 ). Thus, we have shown that (x 0 ) = E 0 . In view of above claims and Lemma 3.3, we see that {E 0 } is an isolated invariant set in X and W s (E 0 ) ∩ X 0 = ∅, where W s (E 0 ) is the stable set of E 0 , and {E 0 } is acyclic in M ∂ . In view of [25, Theorem 1.3.1], we see that {P n } n≥0 is uniformly persistent with respect to (X 0 , ∂X 0 ) in the sense that there existsζ > 0 such that where d is the norm-induced distance in R 12 . By [14, Theorem 3.7], we know that P has a global attractor A 0 in X 0 . Since A 0 = P A 0 , we have that I 0 i > 0, C 0 i > 0, B 0 i > 0, i = 1, 2, for all Ψ(t)A 0 . Then B 0 ⊂ X 0 and lim t→∞ d(Ψ(t)φ, B 0 ) = 0, ∀x 0 ∈ X 0 . Then there exists an ζ > 0 such that for any solution (S 1 (t), Furthermore, [25,Theorem 1.3.6] implies that P has a fixed point is an ω-periodic solution of (1). We can further show that due to the similar arguments to those in Lemma 3.2. The rest of the mathematical arguments were motivated by the work [21].

Numerical Simulations
This section is devoted to the numerical investigation of seasonal and spatial influences on the dynamics of system (1). Except the migration rates, the parameters used in our numerical simulations were given in Table 1. In order to study the spatial influence on system (1), we first put m w 21 = m w 12 = 0, w = S, I, C, R, B, into system (1), and we obtain the following model in a single patch i: where either i = 1 or i = 2. Next, we define the basic reproduction number, R (i) 0 , for system (22). Let By the same ideas in Lemma 2.3 (see also [22,Theorem 2.1]), we let W (i) (t, s, λ), t ≥ s, s ∈ R be the evolution operator of the linear ω-periodic system on R 3 , where ρ(W (i) (ω, 0, λ)) is the spectral radius of W (i) (ω, 0, λ). Use the property (23) and the parameters in Table 1, we can numerically observe that R (1) 0 < 1 and R (2) 0 < 1 whenever 0 ≤ d 1 , d 2 ≤ 1 (see Figure 1 ). This reveals that the disease cannot spread in a single patch with parameters in Table 1 and 0 ≤ d 1 , d 2 ≤ 1. Table 1. Parameter values (except migration rates) used in our numerical simulations, where the values of d1 and d2 will be given.

Parameters
Mean value References θ1, θ2 0.2827 [11] δ I 1 , δ I 2 0.06 [15] δ C 1 , δ C 2 0.004 [15] δ Q 1 , δ Q 2 0.0033 [15] σ1, σ2 0.04 [2] I 1 , I 2 0.1 [15] C 1 , C 2 0.001 [16] α g1, g2 0.014 [17] Next, we numerically investigate the influence of migration on the transmission of typhoid. We will show that typhoid can become epidemic if we incorporate suitable migration rates between patch 1 and patch 2 with parameters in Table  1 into system (1) whenever d 1 and d 2 are close to 1. That is, we observe that typhoid can become epidemic with suitable migration rates between patch 1 and patch 2, although it cannot be epidemic in each isolated single patch (without migration). More precisely, the immigration rates are chosen as follows: where the migration rates from patch 2 to patch 1 are much bigger than those from patch 1 to patch 2. We first assume the parameters in system (1) take those in Table 1 and (24), then the basic reproduction number R 0 of system (1) can be numerically computed by the ideas in Lemma 2.3, and an interesting phenomenon occurs that R 0 > 1 when d 1 , d 2 are close to 1 (see Figure 2 ); nevertheless, R (i) 0 < 1, i = 1, 2 (see Figure 1 ). This indicates that migration between patches does play a central role in the transmission of typhoid. Finally, we also numerically show the extinction of the isolated single patch system (22) with i = 1, 2, respectively, but persistence of the two-patch system (1) occurs after the migration rates are included. The parameters we used are in Table 1 and (24). Putting d 1 = 0.96 (resp. d 2 = 0.95), and the extinction of system (22) with i = 1 (resp. i = 2 ) is illustrated in Figure 3 (resp. Figure 4). After we put d 1 = 0.96 and d 2 = 0.95 and the immigration rates are chosen in (24), persistence of system (1) occurs (see Figure 5).  Table 1 and we observe that R   Table 1 and (24).  Table 1 and only the dynamics of (I1(t), C1(t), B1(t)) are shown.  Table 1 and only the dynamics of (I2(t), C2(t), B2(t)) are shown.  Table 1 and (24), and only the dynamics of (I1(t), C1(t), B1(t), I2(t), C2(t), B2(t)) are shown.

Conclusion
Typhoid can be a serious problem in several developing countries, and the infection of typhoid is through both direct (i.e. human-to-human) transmission and indirect (i.e. environment-to-human) transmission, which is associated with the ingestion of contaminated food or water. Asymptomatically infected individuals are not experiencing any symptoms but they can shed S. Typhi into the environment for years, thereby sustaining transmission [15,16]. A higher incidence of typhoid fever [12,13] is usually observed during the rainy season since a large increase of flooding may occur by the rainfall and the water resources is contaminated by the excreta with pathogenic bacteria. It is also reported that people living near water bodies may have a higher rate of typhoid infection [4]. Thus, this study constructed and analyzed a system that incorporates seasonal effects into a two-patch model describing the transmission of Typhoid fever in a seasonally and spatially variable environment, in which the bacteria in the environment was included and the population of human is divided into five classes, namely, susceptible individuals, infected individuals, carrier individuals, individuals under treatment and recovered individuals.
The well-posedness of model (1) (see Lemma 2.2) is first established, and the associated reproduction number, R 0 is also provided. Since the population Q i , i = 1, 2, in system (1) is supposed to be on treatment and it cannot move in the environment, leading to that those terms −m Q 12 Q 1 + m Q 21 Q 2 and m Q 12 Q 1 − m Q 21 Q 2 must be ignored. This makes mathematical analysis more difficult and we can prove the extinction of Typhoid fever only when R 0 < 1 and the additional condition ρ i = 0, i = 1, 2, is imposed, where ρ i represents the immunity waning rate in patch i (see Theorem 3.1 (i)). When R 0 > 1, the Typhoid fever can persist in the environment if one of the initial values of infected individuals ( I 0 i ), carrier individuals (C 0 i ) and bacteria (B 0 i ) is non-zero, for some patch i (see Theorem 3.1 (ii)).
The numerical results in Section 4 further indicate that migration rates between patches may reverse the outcome of the persistence of Typhoid fever. The used parameters of system (1) takes the form in Table 1 with d 1 = 0.96 and d 2 = 0.95. When migration rates are further chosen as those in (24), it is observed that Typhoid fever persists ( see Figure  5), however, typhoid may become extinct when migration rates are changed into zeros (see Figure 3 and Figure 4). Thus, our simulation confirms that spatial effects do play an important role in the transmission of Typhoid fever.
Finally, the influences of seasonal factors on the transmission of typhoid are also performed numerically. For this purpose, the used parameters of system (1) takes the form in Table 1 with d 1 , d 2 ∈ [0, 1], then the dependence of R 0 on d 1 , d 2 ∈ [0, 1] is illustrated in Figure 2. The basic reproduction number, [R 0 ], of the time-averaged autonomous system corresponding to (1) is also calculated. For a continuous periodic function g(t) with the period ω, we define its average as [g] := 1 ω ω 0 g(t)dt. In view of  (8) and (9), respectively, in which β Ci (t), β Bi (t), α Ii (t) and α Ci (t) are replaced by  Figure 2, it follows that one may underestimate the infection risks if the seasonal effects are ignored.