Abstract
A mathematical model for the dynamics of cholera transmission with permissible controls between two connected communities is developed and analysed. The dynamics of the disease in the adjacent communities are assumed to be similar, with the main differences only reflected in the transmission and disease related parameters. This assumption is based on the fact that adjacent communities often have different living conditions and movement is inclined toward the community with better living conditions. Community specific reproduction numbers are given assuming movement of those susceptible, infected, and recovered, between communities. We carry out sensitivity analysis of the model parameters using the Latin Hypercube Sampling scheme to ascertain the degree of effect the parameters and controls have on progression of the infection. Using principles from optimal control theory, a temporal relationship between the distribution of controls and severity of the infection is ascertained. Our results indicate that implementation of controls such as proper hygiene, sanitation, and vaccination across both affected communities is likely to annihilate the infection within half the time it would take through self-limitation. In addition, although an infection may still break out in the presence of controls, it may be up to 8 times less devastating when compared with the case when no controls are in place.
1. Introduction
Cholera is a gastroenteritic infection contracted after consuming the infectious dose or inoculum size of pathogenic Vibrio cholerae. The ability of V. cholerae to invade the gastroenteritis, proliferate, and damage the host characterises their virulence [1]. After invasion, the capacity of the vibrios to cause symptoms in the host characterises their pathogenicity. The transmission of vibrios from their natural habitat, the aquatic reservoirs, is through consumption of Vibrio infested water and seafood. This route forms the primary route of transmission. Cholera can also be transmitted from one person to another through consumption of unhygienic or soiled food that may be infested with pathogenic vibrios from an infected person. This route is characterised as the secondary route of transmission. For both of the mentioned routes, oral ingestion is the portal of entry of the pathogen into the human body.
Various controls for cholera have been recommended by the World Health Organization (WHO). These controls range from preventive measures to treatment protocols. Preventive measures are aimed at averting new infections by preventing the immunologically naive people from consuming or coming into contact with the bacteria. Treatment and control target the infected persons reducing the number of those infected as well as reducing the case fatality rate. Such preventive and control measures are explained in detail with regard to the metapopulation model of two connected communities as follows. First, an Oral Cholera Vaccine (OCV) recently recommended by the World Health Organization (WHO) [2] is now in use. Notably, it was recently used during the cholera outbreak that affected Haiti after the 2010-2011 earthquake [3]. Secondly, sanitation and hygiene reduce the rate of pathogen ingestion. Such controls include chemical treatment and boiling of water for cooking and drinking and proper storage and preparation of foods to prevent contamination. Thirdly, for the infected case management and treatment procedure which involves administering oral rehydration salts (ORS) to restore ion balance, intravenous administration of fluids in serious cases and use of antibiotics are recommended.
Mathematical modelling of cholera transmission dynamics dates back to the 1970s with the work done by Capasso and Paveri-Fontana [4], but more profound work has been done in the last 1 to 2 decades. Most of the work done so far focuses on specific communities and not the metapopulation framework of the infection. An extensive highlight of the metapopulation framework in epidemic models involving both cross community infections and exchange of populations between communities is given in [5] but not a single case is specific to cholera. Optimal control of the infection has only been studied recently in [6] when comparing cholera transmission between two separate Indian communities of Bogra and Calcutta. It is also important to note that optimal control of cholera in single community settings was studied in [6, 7]. To the best of our knowledge, no optimal control study has been considered in the metapopulation framework for cholera transmission to date. However, general modelling of control of epidemics in metapopulations was recently done using an SIS model [8]. In this same work, the authors highlighted the likely difficulty and mathematical intractability to be faced if a SIR or SIRS model is to be used. In this paper, we consider a SIRS model with controls since cholera transmission typically follows such dynamics. In addition a model of cholera transmission in a metapopulation setting published recently [9] highlights the detrimental effect of movement of epidemiologically naive and infected individuals across communities. In [9] however, the movement of recovered individuals was not considered. Given the length of the modelling time, we believe such movement is also possible and it is given consideration in this paper. We aim to consider the dynamics of closed and isolated communities in the presence of controls.
The paper is organised as follows. In Section 2, we present a metapopulation model with permissible controls. In the same section, vital mathematical analyses of the model are given including the community specific reproduction numbers followed by numerical results in Section 3 and in Section 5, we conclude the paper.
2. Mathematical Model
In the model, we consider two routes of transmission. The primary route is characterised by consumption of Vibrio infected water from aquatic sources. The secondary route (also referred to as person-to-person transmission) is characterised by consumption of food that may be contaminated with vibrios from faecal matter. The human population is subdivided into three compartments (, , and ) depending on their status with respect to the infection. The susceptible are those who are at risk of contracting cholera either through contaminated water or by the secondary route. Once infected, individuals move into compartment . Those who recover from the infection move into compartment at rates and for the first and second compartments, respectively. Infected persons who show symptoms receive oral rehydration salts to restore the ion balance and this increases the recovery rate by a rate . This increased recovery is assumed to be similar in both communities. The infection confers some temporary immunity which wanes at a rate . In the infection dynamics, the disease may be endemic or nonendemic. In the former case, the immunity of those infected wanes at a faster rate resulting in a SIRS type of model compared to a SIR model in the latter case. It is assumed that the time delay between consumption of Vibrio infected food or contaminated water and the commencement of infectiousness is negligible (see also [4]). The infectious dose or inoculum size must be consumed if an immunologically naive person is to be infected. Given a high concentration of the infectious dose ( [10]) and the relatively low probability of infection, the dose-response relationship is given by the Monod function for , where (the inoculum size) is concentration of the pathogen required to cause 50% chance of infection. and are the specific contact rates of individuals in the first and second communities with the aquatic reservoirs in their corresponding communities. The use of a saturated incidence function for ensures boundedness of the incidence rate and indicates that the increase in incidence rate is more gradual than linear. In the model we denote vaccination at any time by . We note that the immunologically naive individuals once vaccinated move straight to the compartment of those recovered. The likelihood of infection is reduced by fractions and related to contact with the aquatic reservoir and person-to-person contact, respectively. Here it is assumed that fraction uses clean treated water and fraction practises proper hygiene. We assume that the implementation of controls is by a government health organisation independent of the community specific health groups. Therefore, the implementation of controls may not be greatly influenced by the level of living condition in the adjacent communities but the need to contain the infection.
The flow diagram of disease progression in two communities is given in Figure 1.
The terms and are the incidence functions for the first and second communities, respectively. These two terms describe the rate of apparition of new cholera cases in the respective communities.
The parameters , , and are the respective rates of movement of susceptible, infected, and recovered individuals from the first community to the second community. Similarly , , and are the rates of movement of the susceptible, infected, and recovered individuals from the second to the first community. The movement across communities is assumed to be instantaneous and therefore there is no change in epidemiological state during travel [5]. When either the susceptible, infected, or recovered individuals move from one community to the other, they follow the dynamics of the destination community. The secondary infections are only generated from within a community and there is no cross community infection.
The populations in the first and second communities are replenished through recruitments, at rates and , respectively. Individuals in the two communities are, however, subjected to natural mortality at rates and , respectively. The infected individuals in the two communities suffer disease induced mortality at rates and , respectively. The pathogen concentration in the aquatic environment is replenished through shedding by the infected individuals at rates and , for the two respective communities.
The model system of equations for the two communities in presence of controls is given bywhere , , , and . Note that all constants in the balance equations are nonnegative. In addition and are positive, indicating that, in absence of faecal contribution from infected individuals, the bacteria cannot sustain themselves in the aquatic environment; see also [4]. The initial conditions of the model are , , , , , , , and and are all nonnegative.
2.1. Model Analysis
The solutions to model system of (1a)–(1h) with nonnegative initial conditions are all nonnegative and bounded. Interested readers can investigate positivity and boundedness of the solutions using the method outlined in ([5, Chapter 3]). We note that if both communities are free of the infection, no treatment control protocol may be implemented. However, vaccination may still be in place as will be indicated in the steady states. Although the permissible controls vary with time, to analyse the steady states, we assume that the controls are constant thereby analysing a nonautonomous system of differential equations; see also [7]. Therefore, the disease-free equilibrium is given bywhereThe term depicts the proportion of susceptible individuals who move back and forth in compartments and . Therefore, is the fraction of susceptible individuals who do not move from their home compartments. We note also that the proportions of the susceptible and fall off as quadratic terms of the vaccination (in the denominator) and increase linearly (in the numerator). Therefore, the higher the vaccination coverage, the lower the fraction of the population that remains naive to the infection.
The community specific disease threshold numbers can be obtained using the next generation matrix method outlined in [20]. When computing the disease thresholds however, we assume that the controls are constant so that the model system of equations is nonautonomous:for the first community andfor the second community. Then the model basic reproduction is given by
Theorem 1. The disease-free steady state (2) of model of (1a)–(1h) is globally asymptotically stable whenever and unstable otherwise.
The proof of Theorem 1 can be given using the approach to the proof for Lemma in [9]. However, in this case the community specific reproduction numbers to be used are (4) and (5).
2.2. Optimal Control
The general procedure of optimal control process in an epidemiological model involves the following processes: identifying permissible controls applicable to the model, setting up the objective function with controls, constructing the Hamiltonian, evaluating costate variable (adjoint functions), and identifying the threshold controls that minimise the Hamiltonian. This optimal control minimisation procedure follows Pontryagin’s Maximum/Minimum principle [21].
The objective function for our minimisation problem is given byThe coefficients , , , , and are the coefficients associated with the costs over a finite period of time . and indicate the cost associated with minimising the infection in the first and second communities, respectively. , , and are relative cost weights for the respective control measures. We assume that the cost of controls is nonlinear, hence the use of the quadratic terms. The main goal is to minimise the number of those infected in both communities at the same time minimising the cost of implementing the controls. In this respect, we seek optimal controls , , and such thatwhereTo evaluate the integral in (7), we use the Hamiltonian constructed for the model system of equations with the necessary adjoint functions. Let , , , , , , , and be the adjoint functions or costate variables associated with the evolution of the corresponding state variables. We multiply each of the adjoint functions with the right side of the equation describing the evolution of each state variable:
Theorem 2 (see [22]). Let , , and be the optimal controls for system (1a)–(1h), the state space at equilibrium, and positive semidefinite piecewise differentiable functions for all . If one supposes that, for all ,thenholds for the optimal controls , , and .
To find the differential equations with respect to the associated adjoint functions, we differentiate the Hamiltonian with respect to each of the state variables and obtainwith transversality conditionWe note that for given transversality condition (14) the following hold:The optimal controls are characterised by the following expressions:The standard control arguments on the controls [23] are such thatThe upper bound of indicates that drinkable water is least likely to be responsible for the oral faecal route of transmitting the infection, especially if water is chlorinated, disinfected, or boiled:The value is characteristic of a perfectly effective vaccine. Quite often cholera vaccines have low protective efficacy and have adverse effects associated with them. In one study by the Public Health Agency of Canada [24], the efficacy of the cholera and diarrhoeal vaccine was observed to range between 64% and 85% against Vibrio cholerae 01 El Tor. Therefore, the upper bound () of this control may not be necessarily attainable:Similarly, would signify no transmission of the pathogen through consumption of foods and that all foodstuffs are well prepared and hygienically stored and free from contamination. Unfortunately, this is not usually the case in communities where cholera is endemic.
Differentiating with respect to each of the admissible controls, we obtainThe control characterisations, , , and , of the optimal controls , , and are obtained from :It can be noted that the control characterisations are inversely proportional to the associated costs. The costs are therefore influential in establishing the levels and trajectories of the controls.
The optimal values of the controls are the ones anticipated to contain the disease. This is also in accordance with the other model parameters describing the dynamics of the disease. The incidence of the disease is also key in determining the trajectory of the state variables for the model without controls, the model with controls, and the control variables over time.
3. Numerical Simulations
To numerically solve the optimal control problem, we modify a code developed by Lenhart and Workman [22], to accommodate the model system of (1a)–(1h). We note that this optimal control problem is classified as a quadratic programming problem, since the controls are nonlinear and are of a quadratic form. The state space variables are solved forward in time and the adjoints associated with the state variables solved backward in time. The model system of equations is numerically integrated using the ODE integration routine (ode45), a fourth-order Runge-Kutta Method in MATLAB.
3.1. Parameter Estimation
Due to lack of data, the parameter values used are estimated from published literature and some intuitively selected from information related to cholera transmission dynamics. The unit of the parameters is per day except for those where otherwise is indicated. Nominal values are selected from within the specified ranges to cater for the differences in the communities under study as well as the disease transmission dynamics. In the model, we assume that the first community has better living conditions and better facilities. Therefore, the transmission parameters are assumed to be relatively lower and the infection is assumed to be less devastating in the first community. The second community, however, is assumed to be slightly less disadvantaged and the infection more devastating. Because of the difference in assumed living conditions, migration of individuals is assumed to be inclined towards the first community. The parameters related to movement are set such that , , and . It is also necessary to determine the parameters related to the cost in the objective function (7). The estimated cost of a vaccine dose per person from the study in [19] was about $5. This gives an estimated cost of 50 rands within the South African context. We note that access to clean water and proper sanitation are influenced by a number of factors ranging from the cost associated, the standard of water treatment, and the sophistication of the technology used. According to the study conducted in Johannesburg [18], if all households were to access clean in-house piped water, with proper monitoring and in-house sewage connection with partial treatment, it would cost about US $136.5 billion per year. According to Statistics South Africa [25], the 2014 mid-year population estimate in the country is 54 million. This implies that the estimated cost of access to proper sanitation per person per day is approximately US $6.925. This is equivalent to approximately R80. In the case of infection, the average cost per percentage reduction in the infection in a community is influenced by factors such as the cost of oral rehydration salts, advertisement on social media, transportation, and remuneration. Given the difference of severity of the infection in adjacent communities, we assume that the cost is lower in a community with better facilities. We assume an average cost of R200 in a more devastated community compared to R120 in a community with improved facilities. The estimated costs associated with the controls per person are given in Table 2.
4. Sensitivity Analysis
We analyse the sensitivity of the model community specific reproduction numbers ( and ) to the model parameters using the Latin Hypercube Sampling (LHS) scheme with 1000 runs per simulation. The LHS is a Monte Carlo stratified sampling method that allows us to simultaneously obtain an unbiased estimate of the model output for a given set of input parameter values; see [26]. This statistical sampling method for simultaneously carrying out sensitivity and uncertainty analysis was first applied in an epidemiological model by Blower and Dowlatabadi [27]. The approach has since been employed in a number of epidemiological models including [9, 26].
The Partial Rank Correlation Coefficients (PRCCs) for the full range of parameters are shown in the tornado plots, Figure 2(a) for the first community and Figure 2(b) for the second community. The PRCCs for the community specific reproduction numbers in Figures 2(a) and 2(b) are produced using the expressions and . The process associated with parameters with positive (negative) PRCCs has the potential of aggravating (reducing) the infection when increased. From the equations of the community specific reproduction numbers and , the controls parameters are evident. From the tornado plots in Figure 2 for community specific disease thresholds, it is evident that the controls related to increasing access to clean water and enhancing proper hand hygiene reducing the disease threshold are all associated with negative PRCCs. Therefore, with enhancement of processes related to such controls, the subsequent numbers of those infected will be less than that of their predecessors hence reducing the severity of the infection. It should be noted that, of all the considered controls, vaccination and improvement in hygiene have the greatest potential of easily containing the infection. In addition emigration (immigration) of infected individuals out of (into) a community reduces (increases) the disease burden in that community.
(a)
(b)
4.1. Numerical Results
In the model simulations, the population is not constant and the initial conditions are set as proportions as follows: ; ; ; ; ; ; ; .
The algorithm for simulating the system using the “forward-backward sweep method” is adopted from [22]. We outline the steps followed during the simulation for convenience of the reader.
Step 1. The time interval over which the system is to be integrated is subdivided into equal subintervals. The initial guesses of piecewise-constant controls , , and on the interval , for , are assumed.
Step 2. Using the assumed controls , , and and the initial conditions of the model, the system is numerically integrated forward in time according to its system of differential equations in the optimality system.
Step 3. Using the assumed controls , , and and values of , the costate system in accordance with its system of differential equations with the transversality condition is integrated backward in time
Step 4. The controls vector is updated by entering the new state vector and the costate values in the optimal characterisations.
Step 5. Convergence of the iterations is then checked by ascertaining whether successive iterations are negligibly close for some specified level of tolerance. Once this is achieved, then optimal control is said to have been achieved, the iteration is stopped, and the output values are recorded as solutions. Otherwise, the algorithm is set back to Step 2.
4.1.1. Isolated Communities in Presence of Controls
The communities are assumed to be isolated in this case if there is no movement or exchange of individuals between adjacent communities. However, we allow for recruitment from other areas and therefore the recruitment parameters are nonzero. In the simulations, we use the parameters in Table 1 and set all the parameters, describing movement to zero. The trajectories of the infected populations in both communities are given in Figure 3.
(a)
(b)
When all the controls are implemented, our results suggest that the time taken for the disease to be contained may be approximately half of the time required to contain the disease without controls, that is, by self-limitation. Although the infection is less devastating in the first community, it stays relatively longer in the community compared to the second community (cf. the solid lines in Figures 3(a) and 3(b)). It is plausible that, due to a high transmission in the second community, the disease depletes the susceptible pool at a much faster rate and in about 100 days, there are virtually no more susceptible people to infect for the chosen parameter values. In the model with controls however, the disease is contained at approximately the same time in both communities. Taking a closer look at dashed curves in Figures 3(a) and 3(b), one observes that the disease is contained slightly earlier in the first than the second community for the same level of effort. This scenario is a typical comparison of infections with high and low incidence.
The trajectories of the state space variables for the model with control, the adjoint, and the control values converge around relatively the same time. The iteration process terminates when these different variables converge sufficiently.
4.1.2. Connected Communities in Presence of Controls
When the communities are connected, the parameters describing movement are all nonzero and there is exchange of individuals between the communities.
Movement of susceptible and infected individuals is associated with not only an increase in severity but also a long duration of the infection in the two communities; see Figures 3 and 4. In presence of controls however, there is no observed significant difference in the disease severity and duration whether in isolated or connected community cases. For comparison, see dashed lines in Figures 3 and 4.
(a) Infected group in the first community
(b) Infected group in the second community
The profiles of the controls investigated in the model are indicated in Figures 5 and 6. The profiles indicate the need for more fundamental action at the beginning of the epidemic if the disease is to be contained. All the three indicated controls reduce the likelihood of susceptible individuals getting infected. The profiles for access to clean water and vaccination, Figures 5(a) and 5(b), respectively, indicate the importance of the two controls in containing the infection.
(a)
(b)
The profile for the control related to hygiene is deemed the least significant; see Figure 6. However, we reiterate that no effort is insignificant as long as it leads to reduction in the number of infections.
The susceptible groups in both communities are initially characterised by depletion followed by an increase, Figure 7. We note that the rapid depletion of the susceptible population is due to vaccination, which was assumed to be exponential in the model. As a result, the recovered group is replenished and reduces when the control is halted, Figure 8. By the end of the control period, the community is comprised of those susceptible and recovered. Once the infection is contained, it is assumed that the small group of people who could still be carrying the disease may pose a low threat of reintroducing it.
(a) Susceptible group in the first community
(b) Susceptible group in the second community
(a) Recovered group in the first community
(b) Recovered group in the second community
Although the general trajectories of the susceptible groups in the two communities are similar, there is a striking difference that ought to be highlighted. Even though there is an observed initial increase in the trajectory of the susceptible group in the first community, that of the second community only decreases. This trend is attributed to movement which is inclined towards the first community due to the better living conditions assumed, as well as the high incidence of the disease in the second community.
5. Conclusion
In this paper, a metapopulation model for cholera transmission characterised by exchange of individuals between communities is analysed. The community specific disease thresholds are given and their importance determining the existence and stability of equilibria of the model is highlighted. The conditions for annihilating the infection based on specific controls while keeping other parameters constant have been indicated. Our numerical results indicate that, in presence of controls, in case of an outbreak of the disease, the infection may be 8 times less devastating when compared to the case without controls. In addition, the duration of the infection in the community cannot be more than half the time it would in absence of controls. From the model results, we assert that implementation of preventive measures and treatment control protocols reduces the likelihood of recurrence of the infection in the communities. The recurrence of the infection is typical in the model without controls presented in [9], characterised by semiperiodic and synchronous fluctuations of the involved subpopulations.
The model presented in this paper is not short of shortfalls. We acknowledge the fact that controls across communities may not be uniform as assumed in the models. Controls efforts are often influenced by demographics of the community (including ageing which is associated with immune dysfunction) as well as the behaviour of individuals involved in the disease transmission dynamics. In addition, politics greatly affects the minority groups as well as social inequality which impacts the vulnerable groups, urbanization associated with overcrowding, inadequate infrastructure, and economic stochasticity among others. All these affect the bureaucracy associated with implementation of controls, putting up sustainable solutions resulting in heterogeneity. The model also assumes a negligible time difference between infection and becoming infectious as well as apparition of symptoms. Although such simplifying assumptions may be presumed feasible, their investigation is necessary to ascertain the effect on the epidemic size.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors acknowledge with gratitude the support from the Department of Mathematical Sciences, Faculty of Science, University of Stellenbosch, and the African Institute for Mathematical Sciences (AIMS).