Abstract
The management of the Buruli ulcer (BU) in Africa is often accompanied by limited resources, delays in treatment, and macilent capacity in medical facilities. These challenges limit the number of infected individuals that access medical facilities. While most of the mathematical models with treatment assume a treatment function proportional to the number of infected individuals, in settings with such limitations, this assumption may not be valid. To capture these challenges, a mathematical model of the Buruli ulcer with a saturated treatment function is developed and studied. The model is a coupled system of two submodels for the human population and the environment. We examine the stability of the submodels and carry out numerical simulations. The model analysis is carried out in terms of the reproduction number of the submodel of environmental dynamics. The dynamics of the human population submodel, are found to occur at the steady states of the submodel of environmental dynamics. Sensitivity analysis is carried out on the model parameters and it is observed that the BU epidemic is driven by the dynamics of the environment. The model suggests that more effort should be focused on environmental management. The paper is concluded by discussing the public implications of the results.
1. Introduction
The Buruli ulcer disease (BU) is a rapidly emerging, neglected tropical disease caused by Mycobacterium ulcerans (M. ulcerans) [1, 2]. It is a poorly understood disease that is associated with rapid environmental changes to the landscapes, such as deforestation, construction, and mining [3–6]. It is a serious necrotizing cutaneous infection which can result in contracture deformities and amputations of the affected limb [3, 7]. Very little is known about the ecology of the M. ulcerans in the environment and their distribution patterns [3]. The survival of vectors or pathogens in the environment can be directly or indirectly influenced by landscape features such as land use and cover types. These features influence the vector or pathogen's ability to survive in the environment or to be transmitted. In most cases the dynamics of the reservoirs and vector depend on the management of the environment. Research has shown that BU is highly prevalent in arsenic-enriched drainages and farmlands [8, 9].
The lack of understanding of the dynamics of the interactions of humans, the vectors, and the BU transmission processes severely hinders prevention and control programs. However, mathematical models have been used immensely as tools for understanding the epidemiology of diseases and evaluating interventions. They now play an important role in policy making, health-economic aspects, emergency planning and risk assessment, control-programs evaluation, and optimizing various detection methods [10]. The majority of mathematical models developed to date for disease epidemics are compartmental. Many of them assume that the transfer rates between compartments are proportional to the individuals in a compartment. In an environment where resources are limited and services lean, this assumption is unrealistic. In particular, the uptake rate of infected individuals into treatment programs is often influenced by the capacity of health care systems, costs, socioeconomic factors, and the efficiency of health care services. For BU, the number of people admitted for treatment is limited by the capacity of health care services, the cost of treatment, distance to hospitals, and health care facilities that are often few [11, 12]. BU treatment is by surgery and skin grafting or antibiotics. It is documented that antibiotics kill M. ulcerans bacilli, arrest the disease, and promote healing without relapse or reduce the extent of surgical excision [13]. Improved treatment options can alleviate the plight of sufferers. These challenges all stem from the fact that many of the developing countries have limited resources.
The demand for health care services often exceeds the capacity of health care provision in cases where the infected visit modern medical facilities. It will be thus plausible to use a saturated treatment function to model limited capacity in the treatment of the BU; see also [10, 14, 15]. The transmission of BU is driven by two processes: firstly, it occurs through direct contact with M. ulcerans in the environment [1, 16, 17] and, secondly, it occurs through biting by water bugs [18, 19]. In this paper we capture these two modes of transmission and also incorporate saturated treatment. The aim is to model theoretically the possible impact of the challenges associated with the treatment and management of the BU such as delays in accessing treatment, limited resources, and few medical facilities to deal with the highly complex treatment of the ulcer. We also endeavour to holistically include the main forms of transmission of the disease in humans. This makes the model richer than the few attempts made by some authors; see, for instance, [18].
This paper is arranged as follows. In Section 2, we formulate and establish the basic properties of the model. The model is analysed for stability in Section 3. Numerical simulations are given in Section 4. In fact, parameter estimation, sensitivity analysis, and some numerical results on the behavior of the model are presented in this section. The paper is concluded in Section 5.
2. Model Formulation
2.1. Description
The transmission dynamics of the BU involve three populations: that of humans, water bugs, and the M. ulcerans. Our model is thus a coupled system of two submodels. The submodel of the human population is an type model, with denoting the susceptible humans, those infected with the BU, those in treatment, and the recovered. The total human population is given by The submodel of the water bugs and M. ulcerans has three compartments. The population of water bugs is comprised of susceptible water bugs and the infected water bugs . The total water bugs population is given by The third compartment, , is that of M. ulcerans in the environment whose carrying capacity is . The possible interrelations between humans, the water bugs, and environment are represented in Figure 1. As in [14, 15], we also assume a saturation treatment function of the form where is the maximum treatment rate. A different function can, however, be chosen depending on the modelling assumptions. The function that models the interaction between humans and M. ulcerans has been used to model cholera epidemics [20] and the references cited therein. We note that if BU cases are few, then , which is a linear function assumed in many compartmental models incorporating treatment; see, for instance, [21, 22]. On the other hand, if BU cases are many, then a constant. So for very large values of of the uptake of BU patients into treatment become constant, thus reaching a saturation level. The parameters and are the effective contact rates of susceptible humans with the water bugs and the environment, respectively. Here is the product of the biting frequency of the water bugs on humans, density of water bugs per human host, and the probability that a bite will result in an infection. Also, is the product of density of M. ulcerans per human host and the probability that a contact will result in an infection. The parameter gives the concentration of M. ulcerans in the environment that yield chance of infection with BU.
The dynamics of the susceptible population for which new susceptible populations enter at a rate of are given by (4). Some BU sufferers do not recover with permanent immunity; they lose immunity at a rate and become susceptible again. The third term models the rate of infection of susceptible populations and the last term describes the natural mortality of the susceptible populations. In this model, the human population is assumed to be constant over the modelling time with the birth and death rate being the same: where and is a function that models saturation in the treatment of BU.
For the population infected with the BU, we have Equation (5) depicts changes in the infected BU cases. The first term represents individuals who enter from the susceptible pool driven by the force of infection . The second term represents the treatment of BU cases modelled by the treatment function . The last term represents the natural mortality of infected humans.
Equation (6), models the human BU cases under treatment. In this regard, the first term represents the movement of BU cases into treatment and the second term, with rates and , respectively, represents natural mortality and recovery.
For individuals who would have recovered from the infection after treatment, their dynamics are represented by the following equation: The first term denotes those who recover at a per capita rate and the second term, with rates and , respectively, represents the natural mortality and loss of immunity.
The equations for the submodel of water bugs are Equation (8) tracks susceptible water bugs. The first term is the recruitment of water bugs at a rate of . The second and third term model the infection rate of water bugs by M. ulcerans at the rate of and the natural mortality of the water bugs at a rate , respectively. Equation (9) deals with the infectious class of the water bug population. The first term simply models the infection of water bugs and the second term models the clearance rate of infected water bugs , from the environment.
The dynamics of M. ulcerans in the environment are modelled by
The first term models the shedding of M. ulcerans by infected water bugs into the environment and the second term represents the removal of M. ulcerans from the environment at the rate .
System (4)–(10) is subject to the following initial conditions:
It is easier to analyse the models (4)–(10) in dimensionless form. Using the following substitutions: and given that , and , system (4)–(10) when decomposed into its subsystems becomes where , , and . Given that the total number of bites made by the water bugs must equal the number of bites received by the humans, is a constant; see [23].
3. Model Analysis
Our model has two subsystems that are only coupled through infection term. Our analysis will thus focus on the dynamics of the environment first and then we consider how these dynamics subsequently affect the human population. We first consider the properties of the overall system before we look at the decoupled system.
3.1. Basic Properties
Since the model monitors changes in the populations of humans and water bugs and the density of M. ulcerans in the environment, the model parameters and variables are nonnegative. The biologically feasible region for the systems (13)-(14) is in and is represented by the set where the basic properties of local existence, uniqueness, and continuity of solutions are valid for the Lipschitzian systems (13)-(14). The populations described in this model are assumed to be constant over the modelling time.
We can easily establish the positive invariance of . Given that , we have . The solutions of systems (13)-(14) starting in remain in for all . The -limit sets of systems (13)-(14) are contained in . It thus suffices to consider the dynamics of our system in , where the model is epidemiologically and mathematically well posed.
3.2. Positivity of Solutions
For any nonnegative initial conditions of systems (13)-(14), the solutions remain nonnegative for all . Here, we prove that all the stated variables remain nonnegative and the solutions of the systems (13)-(14) with nonnegative initial conditions will remain positive for all . We have the following proposition.
Proposition 1. For positive initial conditions of systems (13)-(14), the solutions , , , , and are nonnegative for all .
Proof. Assume that
Thus , and it follows directly from the first equation of the subsystem (13) that
This is a first order differential equation that can easily be solved using an integrating factor. For a nonconstant force of infection , we have
Since the right-hand side of (18) is always positive, the solution will always be positive. If is constant, this result still holds.
From the second equation of subsystem (13),
The third equation of subsystem (13) yields
Similarly, we can show that and for all and this completes the proof.
3.3. Environmental Dynamics
The subsystem (14) represents the dynamics of water bugs and M. ulcerans in the environment. From the second equation, we have where In this case .
The case yields the infection free equilibrium point of the environmental dynamics submodel given by The submodel also has an endemic equilibrium given by
Remark 2. It is important to note that the is our model reproduction number for the BU epidemic in the presence of treatment driven by the dynamics of the water bug and M. ulcerans in the environment. A reproduction number, usually defined as the average of the number of secondary cases generated by an index case in a naive population, is a key threshold parameter that determines whether the BU disease persists or vanishes in the population. In this case, it represents the number of secondary cases of infected water bugs generated by the shedded M. ulcerans in the environment. determines the infection in the environment and subsequently in the human population. We can alternatively use the next generation operator method [24, 25] to derive the reproduction number. A similar value was obtained under a square root sign in this case. The reproduction number is independent of the parameters of the human population even when the two submodels are combined. It depends on the life spans of the water bugs and M. ulcerans in the environment, the shedding, and infection rates of the water bugs. So, the infection is driven by the water bug population and the density of the bacterium in the environment. The model reproduction number increases linearly with the shedding rate of the M. ulcerans into the environment and the effective contact rate between the water bugs and M. ulcerans. This implies that the control and management of the ulcer largely depend on environmental management.
3.3.1. Stability of
Theorem 3. The infection free equilibrium is globally stable when and unstable otherwise.
Proof. We propose a Lyapunov function of the form
The time derivative of (25) is
When , is negative and semidefinite, with equality at the infection free equilibrium and/or at . So the largest compact invariant set in such that when is the singleton . Therefore, by the LaSalle Invariance Principle [26], the infection free equilibrium point is globally asymptotically stable if and unstable otherwise.
3.3.2. Stability of
Theorem 4. The endemic steady state of the subsystem (14) is locally asymptotically stable if .
Proof. The Jacobian matrix of system (14) at the equilibrium point is given by Given that the trace of is negative and the determinant is negative if , we can thus conclude that the unique endemic equilibrium is locally asymptotically stable whenever .
Theorem 5. If , then the unique endemic equilibrium is globally stable in the interior of .
Proof. We now prove the global stability of endemic steady state whenever it exists, using the Dulac criterion and the Poincaré-Bendixson theorem. The proof entails the fact that we begin by ruling out the existence of periodic orbits in using the Dulac criteria [27]. Defining the right-hand side of (14) by , we can construct a Dulac function We will thus have Thus, subsystem (14) does not have a limit cycle in . From Theorem 4, if , then is locally asymptotically stable. A simple application of the classical Poincaré-Bendixson theorem and the fact that is positively invariant suffice to show that the unique endemic steady state is globally asymptotically stable in .
3.4. Dynamics of BU in the Human Population
Our ultimate interest is to determine how the dynamics of water bugs and M. ulcerans impact the human population. The overall goal is to mitigate the influence of the M. ulcerans on the human population. We can actually evaluate the force of infection so that This means that the analysis of submodel (13) is subject to . Our force of infection is thus now a function of the reproduction number of submodel (14) and is constant for any given value of the reproduction number. Figure 2 is a plot of versus .
Using the second equation of system (13), we can evaluate so that From the third equation of (13), we have
Substituting and in the first equation of (13) at the steady state yields a quadratic equation in given by where Clearly our model has two possible steady states given by where are roots of the quadratic equation (33). We note that if , we have and . By Descartes’ rule of signs, irrespective of the sign of , the quadratic equation (33) has one positive root; the endemic equilibrium . We thus have the following result.
Theorem 6. System (13) has a unique endemic equilibrium whenever .
Remark 7. It is important to note that when subsystem (14) is at its infection free steady state then the human population will also be free of the BU. We can easily establish the BU free equilibrium in humans as . The existence of is thus subject to the water bugs and the environment being free of M. ulcerans.
3.4.1. Local Stability of
Theorem 8. The disease free equilibrium whenever it exists is locally asymptotically stable if and unstable otherwise.
Proof. When , then either there are no infections in the water bugs or they are simply carriers. So exists. The Jacobian matrix of system (13) at the disease free equilibrium point is given by The eigenvalues of are , , and . We can thus conclude that the disease free equilibrium is locally asymptotically stable whenever .
3.4.2. Local Stability of
Theorem 9. The unique endemic equilibrium point is locally asymptotically stable for .
Proof. The Jacobian matrix at the endemic steady state is given by
If we let , then the eigenvalues of are given by the solutions of the characteristic polynomial
where
Using the Routh-Hurwitz criterion, we note that , and . The evaluation of yields
This establishes the necessary and sufficient conditions for all roots of the characteristic polynomial to lie on the left half of the complex plane. So the endemic equilibrium is locally asymptotically stable.
In the next section we establish the global stability of the endemic equilibrium using the approach according to Li and Muldowney [28] based on monotone dynamical systems and outlined in Appendix A of [29, 30].
3.4.3. Global Stability of the Endemic Equilibrium
We begin by stating the following theorem.
Theorem 10. If , system (13) is uniformly persistent in , the interior of .
The existence of , only if , guarantees uniform persistence [31]. System (13) is said to be uniformly persistent if there exists a positive constant such that any solution with initial conditions satisfies
The proof of uniform persistence can be done using uniform persistence results in [31, 32].
Theorem 11. If , the endemic equilibrium point of system (13) is globally asymptotically stable when .
Proof. Using the arguments in [28], system (13) satisfies assumptions and in . Let and be the vector field of system (13). The Jacobian matrix corresponding to system (13) is
The second additive compound matrix is given byWe let the matrix function take the form
We thus have
where is the diagonal element matrix derivative of with respect to time andwhere ′ represents the derivative with respect to time.
The matrix can be written as a block matrix so that
where
Let denote the vectors in and let the norm in be defined by
Also let denote the Lozinski measure with respect to this norm. Following [33] we have
where and are the matrix norms with respect to the vector norm and is the Lozinski measure with respect to the norm.
In fact
We now have
The third equation of (13) gives
Substituting (53) into (52) yields
where is the constant of uniform persistence. The inequalities follow Theorem 10.
If we impose the condition , then
where with
Hence
The imposed condition implies that the infection rate is greater than the recovery rate. The result follows based on the Bendixson criterion proved in [28].
4. Numerical Simulations
In this section we endeavour to give some simulation results for the combined subsystems (13) and (14). The simulations are performed using MALAB, and we set our time in years. We carry out sensitivity analysis to determine the effects of a chosen parameter on the state variables. Specifically, we chose to focus on the parameters that make up the model reproduction number because we are interested in parameters that aid the reduction of the BU epidemic. We now give a brief exposition on parameter estimation.
4.1. Parameter Estimation
The estimation of parameters in the model validation process is a challenging process. We make some hypothetical assumptions for the purpose of illustrating the usefulness of our model in tracking the dynamics of the BU. Demographic parameters are the easiest to estimate. For the mortality rate , we assume that the life expectancy of the human population is 61 years. This value has been the approximation of the life expectancy in Ghana [34] and is indeed applicable to sub-Saharan Africa. This translates into per year or equivalently per day. The Buruli ulcer is a vector borne disease and some of the parameters we have can be estimated from literature on vector borne diseases. Recovery rates of vector borne diseases range from to per day [35]. The rate of loss of immunity for vector borne diseases ranges between and per day [35]. Although the mortality rate of the water bugs is not known, it is assumed to be per day [18]. We assume that we have more water bugs than humans so that . The remaining parameters were reasonably estimated based on literature on vector borne diseases and the intuitive understanding of the BU disease by the first two authors.
4.2. Sensitivity Analysis
Many of the parameters used in this paper are not experimentally obtained. It is thus important to test how these parameters affect the output of the variables. This is achieved by employing sensitivity and uncertainty analysis techniques. In this subsection, we explore the sensitivity analysis of the model parameters to find out the degree to which the parameters influence the outputs of the model. We determine the partial correlation coefficients (PRCCs) of the parameters. The parameters with negative PRCCs reduce the severity of the BU epidemic while those with positive PRCCs aggravate it. Using Latin hypercube sampling (LHS) scheme with simulations for each run, we investigate only four of the most significant parameters. These parameters influence only submodel (14). The scatter plots are shown in Figure 3.
(a)
(b)
(c)
(d)
Figures 3(a) and 3(b) depict parameters with a positive correlation with the reproduction number. They show a monotonic increase of as and increase. This means that, to curtail the epidemic, the reduction in the shedding rate and infection of water bugs by M. ulcerans is of paramount importance. On the other hand, Figures 3(c) and 3(d) show a negative correlation with the reproduction number. This means that the clearance of the water bug and the M. ulcerans in the environment will reduce the spread of BU epidemic.
A more informative comparison of how the parameters influence the model is given in Figure 4. The tornado plot shows that the parameter affects the reproduction more than any of the other parameters considered. So interventions targeted towards the reduction in the shedding rate of M. ulcerans into the environment will significantly slow the epidemic.
4.3. Simulation Results
To validate our mathematical analysis results, we plot phase diagrams for less than 1 and greater than 1 for the environmental dynamics. The global properties of the steady states are confirmed in Figures 5(a) and 5(b). The black dots show the location of the steady states.
(a)
(b)
Figure 6 shows a three-dimensional phase diagram for the human population dynamics. The existence of the endemic equilibrium, when , is numerically shown here. The plot shows the trajectories of parametric solutions of (13) for randomly chosen initial conditions. The position of the endemic equilibrium point is indicated on the diagram.
To determine how the infection of the water bugs translates into the transmission of BU in humans, we plot the fraction of BU in humans over time while varying , the clearance rate of bacteria from the environment. Figure 7(a) shows how the infections of BU in humans change with variations in the value of . The infections are evaluated as the number of infected humans. The figure shows that, as the clearance of the bacteria from the environment increases, this translates into a reduction in the number of infected human cases. Similar results can be obtained if the parameter is considered, as shown in Figure 7(b). So interventions to reduce the impact of the epidemic on humans can also be instituted through the reduction of bacteria and water bugs in the environment. It is important to note that the practicality of such an intervention is a mirage. This has worked for other vector borne diseases such as malaria.
(a)
(b)
We also explore the role played by direct M. ulcerans infection on the transmission dynamics of BU. Research has shown that antecedent trauma has often been related to the lesions that characterize BU [36]. Figure 8 shows how the proportion of infected humans changes with increasing transmission rate through direct contact with the environment. The figures show that, as the transmission rate of direct contact of humans with the environment increases, the proportion of the infected also increases. This has direct implications on how humans interact with the environment.
The shedding rate of M. ulcerans into the environment and the treatment rate of humans are important to consider. In Figure 9(a) we observe that increasing the amount of M. ulcerans shed in the environment increases the infections of the BU disease in the human population. While this may sound very obvious, the quantification of the effects thereof is of particular significance. We also note that there are benefits in increasing the maximum threshold with regard to treatment. An increase in the value of will lead to a decrease in the number of infections in the human population.
(a)
(b)
5. Conclusion
We present a deterministic model whose main aim was to capture the two potential routes of transmission and treatment uptake in a resource limited population. This uptake is not linear, and hence we propose a response function of the Michaelis-Menten type. Because of the nature of the infection process, our model is divided into two submodels that are only coupled through the infection term. The model is analysed by determining the steady states. The analysis is done through the submodels. The model in this paper presented a unique challenge in which the infection in one submodel takes place at the steady state of the other submodel. The model analysis is carried out in terms of the model reproduction number . Numerical simulations are carried out. The model parameters were estimated from literature and sensitivity analysis was done because not much of the disease is understood and parameter estimation was difficult. Through the simulations, changes in the number of BU cases in the human population were determined for different values of the clearance rates of the water bugs and M. ulcerans. Our main result is that the management of BU depends mostly on the environmental management, that is, clearance of the bacteria from the environment and reduction in shedding. This in turn will reduce the infection of the water bugs that transmit the infection to humans. As mentioned earlier, clearance of M. ulcerans and water bugs is not practically feasible, but non the less very informative. Research in malaria now looks at vector control, sterilisation, and genetic modification of the mosquito. Such an approach could be beneficial with regard to the control of water bugs.
This model presents the very few attempts to mathematically model BU. A lot of additional extensions can be made. The model can be transformed into a delay differential equation dynamical system to capture treatment delays that are often fatal to BU victims. Social interventions such as educational campaigns can be included in the model to capture various campaigns and initiatives to stop the disease. Finally, this model can be used to suggest the type of data that should be collected as research on the ulcer intensifies.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
E. Bonyah and I. Dontwi acknowledge, with thanks, the support of the Department of Mathematics, Kwame Nkrumah University of Science and Technology. F. Nyabadza acknowledges, with gratitude, the support from the National Research Foundation, South Africa, and Post Graduate and International Office, Stellenbosch University, for the production of this paper.