Analysis of a model for waterborne diseases with Allee effect on bacteria

A limitation of current modeling studies in waterborne diseases (one of the leading causes of death worldwide) is that the intrinsic dynamics of the pathogens is poorly addressed, leading to incomplete, and often, inadequate understanding of the pathogen evolution and its impact on disease transmission and spread. To overcome these limitations, in this paper, we consider an ODEs model with bacterial growth inducing Allee effect. We adopt an adequate functional response to significantly express the shape of indirect transmission. The existence and stability of biologically meaningful equilibria is investigated through a detailed discussion of both backward and Hopf bifurcations. The sensitivity analysis of the basic reproduction number is performed. Numerical simulations confirming the obtained results in two different scenarios are shown.


Introduction
Mortality from infectious diseases is still high worldwide; even if it has declined in highincome countries, it represents a crucial issue in low-income countries. A fundamental property of infectious diseases, including diseases caused by waterborne pathogens, is that these complex interactions always result from an infectious individual or environmental source transmitting the pathogen to a susceptible individual. Examples of diseases in common waterborne infections include Cholera, Giardia, Cryptosporidium, Campylobacter, Typhoid and Paratyphoid fevers, hepatitis A and E, norovirus, rotavirus and many others [1]; these can be caused by a variety of pathogenic microbes (bacteria, protozoa, etc.) in contaminated water. Unfortunately, waterborne diseases remain a serious public health concern today, resulting in more than 3.4 million deaths a year according to WHO estimates. Severe waterborne disease outbreaks continue to occur such as Cholera in Haiti from October 2010 to January 2014; the Zimbabwe Cholera epidemic from August 2008 to July 2009. Several different factors must be considered in attempting to understand waterborne disease dynamics including different transmission pathways such as the ingestion of contaminated water or food or consumption and use of fecally contaminated water. Due to their huge impact on public health and social and economic development, waterborne diseases have been the subject of extensive studies in clinical, experimental and theoretical fields. A number of different approaches (mathematical modeling, analysis, simulation) have been used for modeling the disease transmission, and a large number of mathematical models have been published [10,12,20,21,24]. Such models provide challenges and ideas in many other fields of applied mathematics such as ecology, economics in which nonlinear mathematical models having a similar structure are considered [2][3][4][5]22,23,25]. The mathematical models allow to obtain an estimate for the behavior of epidemics and to predict the asymptotic behavior of infection in order to take suitable actions to control epidemics. On introducing a few number of infected in a population of susceptible individuals, the question that arises is if the epidemic will persist or die out. To this aim, the stability of the disease-free equilibrium (i.e. the equilibrium with no infection) and of the endemic equilibria (i.e. equilibria with nonnull population components) are analyzed. In particular: if the disease-free equilibrium is stable, then epidemic will decay; if endemic equilibria exist and are stable, then epidemic will persist. A limitation of current modeling studies in waterborne diseases, however, is that the intrinsic dynamics of the waterborne pathogens are poorly addressed, leading to incomplete, and often, inadequate, understanding of the pathogen evolution and its impact on disease transmission and spread. Most cholera models in the literature are based on the standard assumption that the pathogens cannot sustain themselves in the absence of human contribution. The rate of change for the bacterial density, in this case, is simply the sum of a positive contribution from the infected human population and a negative contribution due to natural death of the bacteria, and both contributions are assumed as linear. But there have been strong empirical evidences that the pathogens "can independently persist in the environment and, consequently, their intrinsic growth and decay may play an essential role in shaping cholera epidemics" (see [28] and further references therein). Studies from a few research groups, however [14,27], have considered that the interaction between the pathogen and the host could be more complicated than linear. Also, the bacterial growth outside of human hosts does not have to follow linear dynamics; all the cited literature assumes more general (nonlinear) functional forms for the bacterial growth. Moreover, many recent works on the impact of climate change on the pathogens growth in relation to waterborne diseases (see the two review articles [16,17]) suggest evidence for such nonlinear models. Richer models could be more adequate in representing such complex phenomena. Recently, a modified model describing the transmission dynamics of a waterborne bacterial infection, which sheds light on the importance of the type of intrinsic bacterial dynamics into the pathogen evolution equation, has been proposed http://www.journals.vu.lt/nonlinear-analysis in [29]. The authors incorporate both direct (human-to-human) and indirect (environmentto-human) transmission pathways, each represented by a bilinear incidence, and examine two types of bacterial intrinsic dynamics: a logistic growth and a cubic growth. In all studies concerning Allee effects in microbial populations [15], a minimal density (the "Allee threshold") is needed to initiate positive population development. Moreover, to better express the shape of indirect transmission, we adopt a nonlinear function (Holling type II functional response) for the incidence [6,10,19]. Such a choice implies that there must be the ingestion of a certain amount of bacteria to contract the infection. These modeling assumptions highlight more interesting and complicated dynamics. In [28], a general function (depicting also the Holling type II functional) has already been considered to express the indirect transmission, but that model simply excludes the logistic growth with a threshold that we analyze in the present paper. The plan of the paper is the following. In Section 2, the mathematical model is introduced, and the basic reproduction number is recovered. The existence and uniqueness of biologically meaningful equilibria is investigated in Section 3. The linear stability analysis of the disease-free equilibrium and of the endemic equilibria (when they exist) is performed in Section 4. In Sections 5.1 and 5.2, respectively, the possible occurrence of backward/forward and Hopf bifurcations is investigated. Section 6 deals with the sensitivity analysis, while numerical simulations on the obtained results are shown in Section 7, and concluding remarks are given in Section 8.

Mathematical model
The model governing a waterborne disease transmission with cubic growth for the bacteria, assuming a Holling type II force of infection, iṡ where the variables S(t), I(t), R(t) denote, respectively, the susceptible, infective and removed individuals, while B(t) denotes the concentration of bacteria in contaminated water, and "·" denotes the time derivative. The cubic growth for the bacteria is essential to induce an Allee effect experimentally demonstrated in many microbial populations [15]. The probability to catch the infection [6,28] is λ(B) = B/(K B + B), being K B (cells/ml) the constant indicating the half saturation rate. The constants appearing in (1)-(4) are positive, and for their biological meaning, we refer to Table 1.
We append to (1)-(4) smooth positive initial data Allee threshold when τ = 0 cells ml −1 Setting N (t) = S(t) + I(t) + R(t) the total host population size at time t, by adding (1), (2) and (4), it turns out thatṄ + µN = µN 0 , and hence i.e. the total host population size is constant. From (5) it follows that each component of the host population is bounded by N 0 and We remark that the cubic growth of bacteria induces an Allee effect on the bacteria. (3), it turns out that B(t) goes into extinction in the absence of infected individuals (i.e. bacteria population can not sustain itself in the absence of infected individuals). If then there exist two positive constants b 1 = (r(k+b) − r 2 (k−b) 2 −4krτ )/(2r), k 1 = (r(k + b) + r 2 (k−b) 2 −4krτ )/(2r) such that, in the absence of infected individuals: (i) if B(t) < b 1 , then B(t) goes into extinction; (ii) if b 1 < B(t) < k 1 , then B(t) exponentially increases up to k 1 ; (iii) if B(t) > k 1 , then B(t) exponentially decreases up to k 1 . In this sense, b 1 and k 1 represent the Allee threshold and the carrying capacity, respectively. Since we are interested in the case in which bacteria can sustain themselves even in the absence of infected individuals, in the sequel, we assume that (6) holds. Model (1)-(4) admits the (unique) disease-free equilibrium E 0 = (N 0 , 0, 0, 0). By using the next generation method [26], we determine the basic reproduction number R 0 defined as "the expected number of secondary cases produced by a typical infected individual during its entire period of infectiousness in a completely susceptible population" [11]. I and B are the only compartments directly related to the disease. The matrices denoting the generation of new infections and the transfer among infectious compartments are, respectively, given by http://www.journals.vu.lt/nonlinear-analysis then the basic reproduction number is given by

Nontrivial equilibria
In this section, we determine the equilibria of (1)-(4) having all positive components. A nontrivial equilibrium is a solution of the system βλ(B)S = (µ + δ)I, δI − µR = 0 (10) given byS beingB > 0 a positive root of the equation where Let us investigate how many equilibria are admitted by system (7)- (10). Setting θ = µN 0 / (µ + δ), from (7) and (8) it follows that From (9) it follows that and in view of the positiveness of I, when (6) holds, it provides the following necessary condition for the existence of a nontrivial equilibrium: In order to find nontrivial equilibria, from (13) and (14) one has to solve g(B) = q(B), B > 0.
The following theorem holds.
Proof. Let us start by remarking that g(B) > 0 for all B > 0, while q(B) > 0 if and only if (15) holds. Furthermore, .
Let us distinguish three cases: In this case, T (0) > 0 and T (B) = 0 admits at least one positive root B ∈ (0, b 1 ). Let us prove that this root is unique. First, T is a convex function of B and has at most two zeros.
is monotone decreasing and admits a unique solution in (0, b 1 ); • If T (B) has a unique root y ∈ (0, b 1 ), then T (B) is monotone in (0, y) and in (y, b 1 ). It follows that T (B) has at most one zero in (0, y) and one zero in has two solutions in (0, b 1 ), say y 1 and y 2 . Let us assume that Since T (B) has only one root in (0, b 1 ) and T is increasing, T (B) has only one root in (0, b 1 ).

If
, we can follow the procedure used in the case R 0 < 1 to conclude that T (B) has only one root in (0, b 1 ).
3. If R 0 > 1, it follows that T (0) < 0. By fixing all the parameters, except for the intrinsic growth rate of the bacteria, r, g(B) does not vary with r, while In conclusion, if (6) holds, there exist a unique endemic equilibrium for B > k 1 for all R 0 , while there exist 0, 1 or 2 endemic equilibria for B < b 1 according to specific values of R 0 , and the thesis is reached. The first thesis of Theorem 1 indicates the possible existence of backward bifurcation. We will investigate it in Section 5.
The following sufficient condition ensures the uniqueness of the endemic equilibrium.
Theorem 2. If (6) holds true and then there exists a unique endemic equilibrium.

Linear stability of the biologically meaningful equilibria
LetS,Ī,B,R be a generic equilibrium. Let us introduce the perturbation fields The system governing the evolution of the perturbation fields iṡ In order to perform the linear stability of the biologically meaningful equilibria, let us consider the linearized version of (18). The spectral equation of L is where I j (j = 1, 2, 3) are the principal invariants of Equation (19) admits the root λ = a 44 < 0 and the roots of the equation The null solution of linearized system is stable if and only if all the roots of (21) have negative real part. The necessary and sufficient conditions, guaranteeing that all the roots of (21) have negative real part, are the Routh-Hurwitz conditions [18] I 1 < 0, I 3 < 0, which imply necessarily that I 2 > 0.
Let us start by investigating the linear stability of the disease-free equilibrium (N 0 , 0, 0, 0). WhenS = N 0 ,Ī =B =R = 0, invariants can be written as The following theorem holds.
Proof. In view of (23), it turns out that I 1 < 0 and Going on to examine now the stability of the endemic equilibriumĒ = (S,Ī,B,R), the principal invariants become Finding necessary and sufficient conditions guaranteeing that (22) are verified is quite complicated due to the presence of a lot of parameters. Sufficient conditions ensuring linear stability of the endemic equilibrium have been found, but we prefer not to report them so as not to weigh down the paper.

Backward bifurcation
In general, it is observed that disease will eradicate (i.e. the disease free equilibrium will be stable) when R 0 < 1, while disease persists when R 0 > 1. The analysis of some epidemic models depicts that a stable endemic equilibrium exists even when R 0 < 1. This phenomenon is known as backward bifurcation. In other words, a backward bifurcation at R 0 = 1 may qualitatively be described as follows. In the neighboring of 1, for R 0 < 1, an equilibrium corresponding to a smaller number of infective individuals, which is unstable, appears, while the disease-free equilibrium and an equilibrium corresponding to a larger number of infective individuals are locally asymptotically stable. Epidemiologically, a backward bifurcation means that it is not enough to simply reduce the basic reproductive number to a value less than one to eradicate a disease. The occurrence of either a forward or a backward bifurcation can have important consequences from the point of view of the disease control and thus for the eradication of the disease. In order to investigate the possibility of backward bifurcation, we choose β as bifurcation parameter, and we obtain the value β * of β at R 0 = 1 given by The Jacobian of (1)-(4) evaluated in E 0 at β = β * is given by Since ψ 3 = 0 is a simple zero eigenvalue of J E0 β * , then the disease free equilibrium E 0 is non hyperbolic equilibrium. Hence, the model at the hand is an eligible candidate to apply the center manifold theory [8]. Let us denote by w = (w 1 , w 2 , w 3 ) T a right eigenvector corresponding to the zero eigenvalue ψ 3 = 0. It follows Furthermore, the left eigenvector θ = (θ 1 , θ 2 , θ 3 ) satisfying w · θ = 1 is given by The left eigenvector is We can thus compute the coefficients from Theorem 4.1 given in [8] as follows: http://www.journals.vu.lt/nonlinear-analysis Considering inā andb only the nonzero terms, it follows that Using the values of left and right eigenvectors, we get .
Observe that the coefficientb is always positive andā > 0 when using Theorem 4.1 in [8], we can state the following theorem.

Investigation on the onset of Hopf bifurcation
Hopf bifurcation is represented as the appearance or disappearance of a periodic orbit through a local change in the stability properties of an equilibrium point. In this section, we explore the possibility of occurrence of Hopf bifurcation and the direction of Hopf bifurcation around the interior equilibrium pointĒ with respect to the bifurcating parameter ξ. If I 1 < 0, in correspondence of an interior equilibrium, Eq. (21) gives the following eigenvalues: In view of Remark 1, let a 33 < 0 and ξ * be the critical value of the Hopf bifurcation parameter given by ξ * = (a 11 + a 22 )(a 11 + a 33 )(a 22 + a 33 ) a 23 (a 22 + a 33 − a 21 ) > 0.
From (27) we find ω, and replacing in (26), we obtain, with Differentiating v with respect to ξ and evaluating for ξ = ξ * , we get so the transversality condition for Hopf bifurcation is satisfied. Thus, it appears that a Hopf bifurcation around the interior equilibrium occurs at ξ = ξ * . The above results can be summarized in the following theorem.
Theorem 5. For model (1)-(4), if a 33 < 0 and I 3 < 0, then there exists a Hopf bifurcation emerging from its positive interior equilibrium when the contribution of each infected person to the population of bacteria, ξ, passes through the critical value ξ * .
We have obtained some conditions under which system (1)-(4) undergoes Hopf bifurcation from an interior equilibrium at ξ = ξ * . Now, we discuss the direction of Hopf bifurcation and the stability of the bifurcating periodic solution. The method used is based on the normal form theory and the center manifold reduction.
Denoting byX ≡ (X 1 , X 2 , X 3 ) T , expansion of Taylor's series up to terms of order 3 (by neglecting the higher-order terms of degree 4 and above), omitting the bar, produces the following system:Ẋ =LX +Ã(X), whereL is given by (20), and Under this transformation, system (30) becomeṡ where P −1 denotes the inverse Then system (32) can be written as There exists a center manifold for (33), which can be represented as Through approximate computation for the center manifold, we obtain 11 + ω 0 f System (33) restricted to the center manifold is given by http://www.journals.vu.lt/nonlinear-analysis The 1st Lyapunov coefficient based on the normal form (34), which determines the stability and direction of periodic solution, is given by where all the derivatives are calculated at the bifurcation point (Y 1 , Y 2 , ξ) = (0, 0, ξ * ).

Sensitivity analysis
To understand the importance of parameters, which are responsible for the transmission of a disease, we perform a sensitivity analysis. It tells which parameters deserve the most numerical attention: a highly sensitive parameter should be carefully estimated as a small variation in it may lead to large changes in the quantity of interest and qualitatively different results. An insensitive parameter does not require as much effort to be estimated. The initial disease transmission is directly related to R 0 ; for this reason, we calculate the sensitivity indices of the basic reproductive number R 0 . Let p a generic parameter of model (1)-(4). The normalized forward sensitivity index of R 0 (which is differentiable with respect to p) is defined as If this index is negative (positive), then the relationship between the parameter and R 0 is of inverse (direct) proportion. The index modulus is an indicator of the size of the effect of changes in that parameter. As we have an explicit expression of R 0 , we derive an analytical expression for the sensitivity of R 0 to each of the different parameters of model (1)-(4). Precisely, From (37) we can see that the sensitivity indices of R 0 with respect to β, ξ, K B do not depend on any parameter values, while the other indices have an obvious structure depending on some parameters. From an initial analysis, R 0 is most sensitive to the contact rate with contaminated water (β) and to the contribution of each infected person to the population of bacteria (ξ) (as shown in (37)). This result shows that any increase (decrease) by a given percentage in β or ξ will increase (decrease) by the same percentage the value of R 0 . The other parameter with highest sensitivity index is K B with Π R0 K B = −1. Increasing (decreasing) K B by a given percentage will decrease (increase) R 0 by the same percentage. The other parameters have a minor effect on R 0 . In fact, from (37) the analytical expression for the sensitivity of R 0 to any of the other parameters, evaluated in absolute value, is less than 1. For all parameters, the sign of the sensitivity indices of R 0 agrees with an intuitive expectation. In order to effectively measure, the relative change in R 0 when the remaining parameters change, one has to evaluate the sensitivity indices at the chosen baseline parameter values. However, following the procedure used by [9], we can replace τ and r (keeping b fixed) by the parameters ζ = rb + τ ,θ = r/τ , where ζ can be seen as the linear loss rate of bacteria. Measuring the sensitivity of R 0 with respect to ζ, keepingθ fixed (allowing r and τ to vary, while their ratio remains fixed), it follows that Π R0 ζ = Π R0 r + Π R0 τ = −1 so providing a good estimate of the joint impact of r and τ . This means that the aggregate parameter ζ is one of the most sensitive parameters. (Similarly, one can proceed replacing τ and b, keeping r fixed, with Γ = rb + τ andθ = b/τ , respectively, and obtain Π R0 With a similar procedure, one can estimate the joint effect of δ and µ by the sensitivity index of R 0 with respect to η = δ + µ, i.e. Π R0 The parameter µ + δ can be seen as the total loss rate of infective individuals. In Table 2, the values of all the parameters displayed in the first line are taken as the baseline values and are the same of Example 1 in Table 3. They are used to evaluate, through (36), the sensitivity indices of the remaining parameters, which are responsible for the transmission and management of a waterborne disease, in relation to the basic reproduction number R 0 . The results of such calculation are presented in the second line of Table 2. The parameters are ordered from the most sensitive to the least.

Numerical simulations and discussion
In this section, we investigate by numerical simulations some different scenarios for the proposed model and discuss the obtained results. Throughout these experiments, we also intend to highlight the impact of the assumption of a nonlinear pathogen growth on the infected dynamics: by tuning the parameter r governing this growth, we will show the higher flexibility of the present model in reproducing different epidemic outbreaks.

Disease dynamics and basic reproduction number
Specifically, we explore two situations: in the first example, we set the model parameters to obtain a basic reproduction number R 0 less than one so that the disease-free equilib-rium is stable; in the second example, we choose the parameters to obtain R 0 greater than one and to have two endemic equilibria in the interval (0, b 1 ). All the parameters values for the two examples are reported in Table 3. In the exploratory study [10], a model has been introduced to capture Cholera transmission within a community. This model has been a predecessor of several more recent ones. As discussed in [13], true values of the parameters in all these models are often difficult or impossible to estimate accurately. That work also reports the most common setting of several key parameter values across the main published models. So we choose accordingly the parameter values for the numerical experiments. Such simulations have the only purpose of showing the ability of the proposed model to represent different scenarios for the waterborne disease evolution. Of course, running such simple models against real data (as we did in a different model setting in [7]) could just reproduce the initial outbreak of an epidemy before the planning and execution of human interventions, like vaccination or water sanitation, that are not yet represented in the current formulation. More realistic models, including terms representing such interventions, could be considered in future work and undoubtedly could benefit from these preliminary studies.
Example 1. We get R 0 ≈ 0.72, b 1 ≈ 10.334, k 1 ≈ 99.667, and the polynomial T (B) introduced in Section 4 has two positive roots, one in the interval (0, b 1 ), the so-called lower branch, and one greater than k 1 (upper branch). These roots correspond to two endemic equilibria with approximate values of the variables S, I, B given by E 1 ≈ (986, 1, 2.795), E up ≈ (659, 25, 104.777), while for the disease-free equilibrium, we find E 0 = (1000, 0, 0). Now, according to the linear stability analysis, E 0 and E up are stable equilibria, but E 1 is unstable. The left panel of Fig. 1 shows trajectories for S(t), I(t) and B(t) starting from the initial point (900, 0, 5) and reaching the disease-free equilibrium (left panel). In the same parameters setting, different initial conditions with I(0) > 0 lead to the endemic equilibrium E up . This can be clearly seen in the right panel of the same Fig. 1.
When R 0 is greater than one, the polynomial T (B)/B, along with the always existent root greater than k 1 , can have either zero or two positive roots. We decided to consider this latter situation to further investigate these endemic equilibria. Then we present a second example.
Example 2. For the parameter setting reported in the second row of Table 3, we get R 0 ≈ 1.23, b 1 ≈ 302, k 1 ≈ 498, and the polynomial T (B) introduced in Section 4 has three positive roots, two in the interval (0, b 1 ) and one greater than k 1 . These roots correspond to the three endemic equilibria E 1 , E 2 and E up with approximate values of the variables S, I, B given by E 1 ≈ (5918, 302, 55), E 2 ≈ (4835, 383, 86), E up ≈ (1275, 646, 579), while for the disease-free equilibrium, we find E 0 = (10000, 0, 0). Now, according to the linear stability analysis, E 0 , E 1 and E 2 are all unstable equilibria, while E up is stable.

Bifurcation scenarios
Now we intend to illustrate through numerical experiments some of the findings of Section 5. First, let us examine in both the reported examples how the chosen parameters (and the resulting value of the reproduction number R 0 ) affect the stability of the endemic equilibria. The bifurcation diagrams in the left and right panels of Fig. 3 refer to Examples 1 and 2, respectively; in both cases, we vary the value of the contact rate β in the range [1,3] and report on the horizontal axis the resulting value of R 0 . It should be noted that in both panels, the upper branch equilibrium B up (that always exists as previously shown) is stable and its value increases along with the value of R 0 . Moreover, the left panel shows how a unique lower branch equilibrium B 1 only exists for R 0 < 1 (and it is always unstable). At R 0 = 1, this equilibrium crosses the DFE, and a backward bifurcation occurs. The right panel reports a more interesting situation: while for R 0 < 1, a single and unstable lower branch equilibrium exists, at R 0 = 1, a second equilibrium emerges (forward bifurcation), that is stable until R 0 = 1.0714 and then becomes unstable. At R 0 = 1.2508, both these equilibria disappear.    Table 3 present some numerical experiments concerning Hopf bifurcation. The analysis reported in Section 5.2 has proved how all the requirements can be met for the onset of a Hopf bifurcation on the endemic equilibrium E 1 = (S 1 .I 1 , B 1 ) with B 1 < b 1 . However, numerical experiments can hardly produce plots illustrating such a situation. This is due to the concurrent presence of a stable equilibrium E up with B up > k that eventually attracts the trajectories of the process. Then the best we can show is the extremely oscillating behaviour of the system trajectories when the parameter values are very close to the ones determining Hopf bifurcation. In the following, we will present some simulations to better clarify our findings. Let us fix all the parameter values as in Example 2 from Table 3, apart from β and ξ. Then we impose for ξ the critical value ξ * defined in Section 5 and solve for (β, B) both (11) and I 3 < 0. With these choices, we find β =   Table 3, apart from β = 1.757, ξ= 12.03. With these choice, the endemic equilibrium E 1 ≈ (8916, 80, 14) is stable, but very close to the Hopf bifurcation limit.
In the left panel, the orbit up to T = 1500, in the right panel, the complete path towards the other equilibrium value Eup ≈ (1704, 615, 587).  Table 3, apart from β = 1.757, ξ = 12.03. With these choice, the endemic equilibrium E 1 ≈ (8916, 80, 14) is stable, but very close to the Hopf bifurcation limit. In the left panel, the orbit starting from I 0 = 80 and converging (very slowly, simulation is run until T = 10000) to E 1 , in the right panel, the orbit starting from I 0 = 30 and more rapidly converging (T = 2000) to the other equilibrium value Eup ≈ (1704, 615, 587).
equilibrium E 1 ≈ (8916, 80, 14) fulfill the conditions for the onset of a Hopf bifurcation. We also evaluate the related first Lyapunov coefficient l 1 by (35) and found a negative value, so that the Hopf bifurcation near the critical point E 1 is supercritical. However, our simulations are only able to show that the system trajectories, even starting very close to E 1 , are attracted towards the stable endemic equilibrium E up ≈ (1704, 615, 587).
To further explore this situation, we then consider a slight perturbation of the parameter values (β = 1.757, ξ = 12.03), so to move the equilibrium point E 1 just within the limits of its linear stability. Even in this case, reported in the right panel of the same figure, the process trajectories, after several oscillations, are eventually attracted by the upper branch equilibrium E up . Figure 5 helps clarify this latter situation by showing the corresponding phase portrait: up to T = 1500, the orbit remains quite close to E 1 (left panel of the figure), but the enlarging oscillations lead the process to converge to E up (right panel of the same figure). Further investigations confirm that the lower branch equilibrium, even when it is stable, has a very small basin of attraction, and the system trajectories are highly sensitive to the chosen initial point: in Fig. 6, we report two phase plane portraits corresponding to the same equilibrium, in the same parameters setting, but with different initial values for the Infected population I 0 . It can be clearly seen how the corresponding orbits are dramatically different and eventually converge to different equilibrium values: E 1 (left panel) or E up (right panel).

Experiments on sensitive parameters
It is clear that the stable endemic equilibrium withB greater than k 1 , present in both scenarios, is heavily depending on the value of the carrying capacity k 1 that represents a lower bound for it. As discussed in the previous section, however, other parameters can play a significant role in the proposed model. Specifically, we further analyze the parameter setting of Example 2 and try to identify effective strategies to control the amplitude and timing of the peaks in both the infective and bacterial populations in the event of large R 0 when the infectious disease spreads out and the only stable equilibrium in the considered model is the endemic one. In Fig. 7, we show trajectories of the Infected population obtained in the parameters setting of Example 2 when only one of the parameters values (β in the left panel, δ in the right one) is modified (increased or decreased) by 10%. These parameters characterize the main contributions to the infected population dynamics. As shown in the figure and as expected by the opposite sign of their contribution in the model, their fluctuations affect the timing of the Infected peak in opposite ways and cause a shift in the peak location without sensibly modifying the equilibrium value. However, changes in β values seem to have a greater impact moving further the peak location. A similar behavior can be observed in the bacterial dynamics (not shown).
To conclude this section, let us present, in the same experiment setting, the impact of different growth parameters r on the infected population. In Fig. 8, we represented the infected-population curves corresponding to several values of r ranging from 0 to 0.25. It is clear that the linear assumption for the pathogen growth (r = 0) corresponds to an abrupt growth of the infected population, while the nonlinear assumption allows to reproduce different scenarios, where the peak in the infected population could be lowered and shifted to the right to better fit experimental data. This could be the case for an epidemic spreading at different rates in different regions of the same country as reported in several experimental studies (see, i.e. [19]).   Table 3 and different values of the r parameter ranging from 0 (dark green line) to 0.25 (blue line).

Conclusions
An ODEs model for the transmission of a waterborne disease has been formulated and investigated, both theoretically and numerically. Its main features are the cubic growth of the bacterial population, that implies an Allee effect, and a Holling type II functional response to better express the shape of the indirect transmission. For this model, the number of biologically meaning equilibria has been obtained along with their linear stability according to the basic reproduction number. The occurrence of a backward bifurcation has been investigated. Precisely, restrictions on the parameters guaranteeing the onset of a backward or a forward bifurcation, that can have consequences for the disease control, have been obtained. Hopf bifurcation has been analyzed by using methods from bifurcation theory and the center manifold theorem. Numerical simulations have confirmed the theoretical analysis and shown the rich dynamics in case of coexistence of two stable equilibria. A sensitivity analysis on R 0 has been also performed to assess the incidence of the different parameters on the initial disease transmission. The obtained results could be a useful contribution to better understand and more realistically describe the dynamics of waterborne diseases. They could also be the starting point for giving insight to future studies aimed at planning intervention strategies.