# A Theoretical Model for the Transmission Dynamics of the Buruli Ulcer with Saturated Treatment

Authors: Ebenezer Bonyah; Isaac Dontwi; Farai Nyabadza
Journal: Computational and Mathematical Methods in Medicine
 (2014)

---

## 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.


---

## Body

## 1. Introduction

The Buruli ulcer disease (BU) is a rapidly emerging, neglected tropical disease caused byMycobacterium 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 theM. 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 killM. 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 withM. 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 Section2, 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 theM. ulcerans. Our model is thus a coupled system of two submodels. The submodel of the human population is an (
S
H
,
I
H
,
T
H
,
R
H
) type model, with S
H denoting the susceptible humans, I
H those infected with the BU, T
H those in treatment, and R
H the recovered. The total human population is given by
(1)
N
H
=
S
H
+
I
H
+
T
H
+
R
H
.
The submodel of the water bugs andM. ulcerans has three compartments. The population of water bugs is comprised of susceptible water bugs S
W and the infected water bugs I
W. The total water bugs population is given by
(2)
N
W
=
S
W
+
I
W
.
The third compartment, D, is that ofM. ulcerans in the environment whose carrying capacity is K
d. 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
(3)
f
(
I
H
)
=
σ
I
H
1
+
I
H
,
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 andM. ulcerans has been used to model cholera epidemics [20] and the references cited therein. We note that if BU cases are few, then f
(
I
H
)
≈
σ
I
H, 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 f
(
I
H
)
≈
σ a constant. So for very large values of I
H of the uptake of BU patients into treatment become constant, thus reaching a saturation level. The parameters β
1 and β
2 are the effective contact rates of susceptible humans with the water bugs and the environment, respectively. Here β
1 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, β
2 is the product of density ofM. ulcerans per human host and the probability that a contact will result in an infection. The parameter K
50 gives the concentration ofM. ulcerans in the environment that yield 50
% chance of infection with BU.Figure 1
A schematic diagram for the model.The dynamics of the susceptible population for which new susceptible populations enter at a rate ofμ
H
N
H 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 (
μ
H
) being the same:
(4)
d
S
H
d
t
=
μ
H
N
H
+
θ
R
H
-
Λ
S
H
-
μ
H
S
H
,
where Λ
=
β
1
I
W
/
N
H
+
β
2
D
/
(
K
50
+
D
) and f
(
I
H
) is a function that models saturation in the treatment of BU.For the population infected with the BU, we have(5)
d
I
H
d
t
=
Λ
S
H
-
f
(
I
H
)
-
μ
H
I
H
.
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 f
(
I
H
). The last term represents the natural mortality of infected humans.Equation (6),
(6)
d
T
H
d
t
=
f
(
I
H
)
-
(
μ
H
+
γ
)
I
H
,
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 μ
H 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:(7)
d
R
H
d
t
=
γ
I
H
-
(
μ
H
+
θ
)
R
H
.
The first term denotes those who recover at a per capita rate γ and the second term, with rates μ
H and θ, respectively, represents the natural mortality and loss of immunity.The equations for the submodel of water bugs are(8)
d
S
W
d
t
=
μ
W
N
W
-
β
3
S
W
D
K
d
-
μ
W
S
W
,
(9)
d
I
W
d
t
=
β
3
S
W
D
K
d
-
μ
W
I
W
.
Equation (8) tracks susceptible water bugs. The first term is the recruitment of water bugs at a rate of μ
N
W. The second and third term model the infection rate of water bugs byM. ulcerans at the rate of β
3 and the natural mortality of the water bugs at a rate μ
W, 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 μ
W, from the environment.The dynamics ofM. ulcerans in the environment are modelled by
(10)
d
D
d
t
=
α
I
W
-
μ
d
D
K
d
.The first term models the shedding ofM. ulcerans by infected water bugs into the environment and the second term represents the removal ofM. ulcerans from the environment at the rate μ
d.System (4)–(10) is subject to the following initial conditions:
(11)
S
H
(
0
)
=
S
H
0
>
0
,
I
H
(
0
)
=
I
H
0
>
0
,
T
H
(
0
)
=
T
H
0
>
0
,
R
H
(
0
)
=
R
H
0
=
0
,
S
W
(
0
)
=
S
W
0
>
0
,
I
W
(
0
)
=
I
W
0
,
D
(
0
)
=
D
0
>
0
.It is easier to analyse the models (4)–(10) in dimensionless form. Using the following substitutions:
(12)
s
h
=
S
H
N
H
,
i
h
=
I
H
N
H
,
τ
h
=
T
H
N
H
,
r
h
=
R
H
N
H
,
s
w
=
S
W
N
W
,
i
w
=
I
W
N
W
,
x
=
D
K
d
,
m
1
=
N
W
N
H
,
and given that s
h
+
i
h
+
τ
h
+
r
h
=
1,  s
w
+
i
w
=
1 and 0
≤
x
≤
1, system (4)–(10) when decomposed into its subsystems becomes
(13)
d
s
h
d
t
=
(
μ
H
+
θ
)
(
1
-
s
h
)
-
θ
(
i
h
+
τ
h
)
-
Λ
~
s
h
,
d
i
h
d
t
=
Λ
~
s
h
-
σ
i
h
1
+
N
H
i
h
-
μ
H
i
h
,
d
τ
h
d
t
=
σ
i
h
1
+
N
H
i
h
-
(
μ
H
+
γ
)
τ
h
,
(14)
d
i
w
d
t
=
β
3
(
1
-
i
w
)
x
-
μ
W
i
w
,
d
x
d
t
=
α
~
i
w
-
μ
d
x
,
where α
~
=
α
N
W
/
K
d,  Λ
~
=
β
1
m
1
(
N
W
,
N
H
)
i
w
+
β
2
x
/
(
K
~
+
x
), and K
~
=
K
50
/
K
d. Given that the total number of bites made by the water bugs must equal the number of bites received by the humans, m
1
(
N
W
,
N
H
) is a constant; see [23].


## 2.1. Description

The transmission dynamics of the BU involve three populations: that of humans, water bugs, and theM. ulcerans. Our model is thus a coupled system of two submodels. The submodel of the human population is an (
S
H
,
I
H
,
T
H
,
R
H
) type model, with S
H denoting the susceptible humans, I
H those infected with the BU, T
H those in treatment, and R
H the recovered. The total human population is given by
(1)
N
H
=
S
H
+
I
H
+
T
H
+
R
H
.
The submodel of the water bugs andM. ulcerans has three compartments. The population of water bugs is comprised of susceptible water bugs S
W and the infected water bugs I
W. The total water bugs population is given by
(2)
N
W
=
S
W
+
I
W
.
The third compartment, D, is that ofM. ulcerans in the environment whose carrying capacity is K
d. 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
(3)
f
(
I
H
)
=
σ
I
H
1
+
I
H
,
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 andM. ulcerans has been used to model cholera epidemics [20] and the references cited therein. We note that if BU cases are few, then f
(
I
H
)
≈
σ
I
H, 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 f
(
I
H
)
≈
σ a constant. So for very large values of I
H of the uptake of BU patients into treatment become constant, thus reaching a saturation level. The parameters β
1 and β
2 are the effective contact rates of susceptible humans with the water bugs and the environment, respectively. Here β
1 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, β
2 is the product of density ofM. ulcerans per human host and the probability that a contact will result in an infection. The parameter K
50 gives the concentration ofM. ulcerans in the environment that yield 50
% chance of infection with BU.Figure 1
A schematic diagram for the model.The dynamics of the susceptible population for which new susceptible populations enter at a rate ofμ
H
N
H 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 (
μ
H
) being the same:
(4)
d
S
H
d
t
=
μ
H
N
H
+
θ
R
H
-
Λ
S
H
-
μ
H
S
H
,
where Λ
=
β
1
I
W
/
N
H
+
β
2
D
/
(
K
50
+
D
) and f
(
I
H
) is a function that models saturation in the treatment of BU.For the population infected with the BU, we have(5)
d
I
H
d
t
=
Λ
S
H
-
f
(
I
H
)
-
μ
H
I
H
.
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 f
(
I
H
). The last term represents the natural mortality of infected humans.Equation (6),
(6)
d
T
H
d
t
=
f
(
I
H
)
-
(
μ
H
+
γ
)
I
H
,
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 μ
H 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:(7)
d
R
H
d
t
=
γ
I
H
-
(
μ
H
+
θ
)
R
H
.
The first term denotes those who recover at a per capita rate γ and the second term, with rates μ
H and θ, respectively, represents the natural mortality and loss of immunity.The equations for the submodel of water bugs are(8)
d
S
W
d
t
=
μ
W
N
W
-
β
3
S
W
D
K
d
-
μ
W
S
W
,
(9)
d
I
W
d
t
=
β
3
S
W
D
K
d
-
μ
W
I
W
.
Equation (8) tracks susceptible water bugs. The first term is the recruitment of water bugs at a rate of μ
N
W. The second and third term model the infection rate of water bugs byM. ulcerans at the rate of β
3 and the natural mortality of the water bugs at a rate μ
W, 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 μ
W, from the environment.The dynamics ofM. ulcerans in the environment are modelled by
(10)
d
D
d
t
=
α
I
W
-
μ
d
D
K
d
.The first term models the shedding ofM. ulcerans by infected water bugs into the environment and the second term represents the removal ofM. ulcerans from the environment at the rate μ
d.System (4)–(10) is subject to the following initial conditions:
(11)
S
H
(
0
)
=
S
H
0
>
0
,
I
H
(
0
)
=
I
H
0
>
0
,
T
H
(
0
)
=
T
H
0
>
0
,
R
H
(
0
)
=
R
H
0
=
0
,
S
W
(
0
)
=
S
W
0
>
0
,
I
W
(
0
)
=
I
W
0
,
D
(
0
)
=
D
0
>
0
.It is easier to analyse the models (4)–(10) in dimensionless form. Using the following substitutions:
(12)
s
h
=
S
H
N
H
,
i
h
=
I
H
N
H
,
τ
h
=
T
H
N
H
,
r
h
=
R
H
N
H
,
s
w
=
S
W
N
W
,
i
w
=
I
W
N
W
,
x
=
D
K
d
,
m
1
=
N
W
N
H
,
and given that s
h
+
i
h
+
τ
h
+
r
h
=
1,  s
w
+
i
w
=
1 and 0
≤
x
≤
1, system (4)–(10) when decomposed into its subsystems becomes
(13)
d
s
h
d
t
=
(
μ
H
+
θ
)
(
1
-
s
h
)
-
θ
(
i
h
+
τ
h
)
-
Λ
~
s
h
,
d
i
h
d
t
=
Λ
~
s
h
-
σ
i
h
1
+
N
H
i
h
-
μ
H
i
h
,
d
τ
h
d
t
=
σ
i
h
1
+
N
H
i
h
-
(
μ
H
+
γ
)
τ
h
,
(14)
d
i
w
d
t
=
β
3
(
1
-
i
w
)
x
-
μ
W
i
w
,
d
x
d
t
=
α
~
i
w
-
μ
d
x
,
where α
~
=
α
N
W
/
K
d,  Λ
~
=
β
1
m
1
(
N
W
,
N
H
)
i
w
+
β
2
x
/
(
K
~
+
x
), and K
~
=
K
50
/
K
d. Given that the total number of bites made by the water bugs must equal the number of bites received by the humans, m
1
(
N
W
,
N
H
) 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 ofM. ulcerans in the environment, the model parameters and variables are nonnegative. The biologically feasible region for the systems (13)-(14) is in R
+
5 and is represented by the set
(15)
Γ
=
{
0
≤
i
w
≤
1
,
0
≤
x
≤
α
~
μ
d
(
s
h
,
i
h
,
τ
h
,
i
w
,
x
)
∈
R
+
5
∣
0
≤
s
h
+
i
h
+
τ
h
≤
1
,
mim
0
≤
i
w
≤
1
,
0
≤
x
≤
α
~
μ
d
}
,
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 d
x
/
d
t
=
α
~
i
w
-
μ
d
x
≤
α
~
-
μ
d
x, we have x
≤
α
~
/
μ
d. The solutions of systems (13)-(14) starting in Γ remain in Γ for all t
>
0. 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 t
∈
[
0
,
∞
). 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 t
>
0. We have the following proposition.Proposition 1.
For positive initial conditions of systems (13)-(14), the solutions s
h
(
t
),  i
h
(
t
),  τ
h
(
t
),  i
w
(
t
), and x
(
t
) are nonnegative for all t
>
0.Proof.
Assume that(16)
t
^
=
sup
⁡
{
t
>
0
:
s
h
>
0
,
i
h
>
0
,
τ
h
>
0
,
i
w
>
0
,
x
>
0
}
∈
(
0
,
t
]
.
Thus t
^
>
0, and it follows directly from the first equation of the subsystem (13) that
(17)
d
s
h
d
t
≤
(
μ
H
+
θ
)
-
[
(
μ
H
+
θ
)
+
Λ
]
s
h
.
This is a first order differential equation that can easily be solved using an integrating factor. For a nonconstant force of infection Λ, we have
(18)
s
h
(
t
^
)
≤
s
h
(
0
)
exp
⁡
[
-
(
(
μ
H
+
θ
)
t
^
+
∫
0
t
^
‍
Λ
(
s
)
d
s
)
]
+
exp
⁡
[
-
(
(
μ
H
+
θ
)
t
^
+
∫
0
t
^
‍
Λ
(
s
)
d
s
)
]
×
[
∫
0
t
^
‍
(
μ
H
+
θ
)
e
(
(
μ
H
+
θ
)
t
^
+
∫
0
t
^
‍
Λ
(
l
)
d
l
)
d
t
^
]
.
Since the right-hand side of (18) is always positive, the solution s
h
(
t
) will always be positive. If Λ is constant, this result still holds.
From the second equation of subsystem (13),
(19)
d
i
h
d
t
≥
-
(
μ
H
+
σ
)
i
h
≥
i
h
(
0
)
exp
⁡
[
-
(
μ
H
+
σ
)
t
]
>
0
.
The third equation of subsystem (13) yields
(20)
d
τ
h
d
t
≥
-
(
μ
H
+
γ
)
τ
h
≥
τ
h
(
0
)
exp
⁡
[
-
(
μ
H
+
γ
)
t
]
>
0
.
Similarly, we can show that i
w
(
t
)
>
0 and x
(
t
)
>
0 for all t
>
0 and this completes the proof.


### 3.3. Environmental Dynamics

The subsystem (14) represents the dynamics of water bugs andM. ulcerans in the environment. From the second equation, we have
(21)
x
*
=
α
~
i
w
*
μ
d
,
i
w
*
=
0
or
i
w
*
=
1
-
1
R
T
,
where
(22)
R
T
=
α
~
β
3
μ
d
μ
W
.
In this case x
*
=
(
α
~
/
μ
d
)
(
1
-
(
1
/
R
T
)
).The casei
w
*
=
0 yields the infection free equilibrium point of the environmental dynamics submodel given by
(23)
E
0
=
(
0,0
)
.
The submodel also has an endemic equilibrium given by
(24)
E
1
=
(
α
~
μ
W
(
R
T
-
1
)
,
μ
d
μ
W
(
R
T
-
1
)
)
.Remark 2.
It is important to note that theR
T is our model reproduction number for the BU epidemic in the presence of treatment driven by the dynamics of the water bug andM. 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 sheddedM. ulcerans in the environment. R
T 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 andM. 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 theM. ulcerans into the environment and the effective contact rate between the water bugs andM. ulcerans. This implies that the control and management of the ulcer largely depend on environmental management.


#### 3.3.1. Stability ofE
0

Theorem 3.
The infection free equilibriumE
0 is globally stable when R
T
<
1 and unstable otherwise.Proof.
We propose a Lyapunov function of the form(25)
V
(
t
)
=
i
w
+
β
3
μ
d
x
.
The time derivative of (25) is
(26)
V
˙
=
d
i
w
d
t
+
β
3
μ
d
d
x
d
t
≤
μ
W
(
R
T
-
1
)
i
w
.
When R
T
≤
1,  V
˙ is negative and semidefinite, with equality at the infection free equilibrium and/or at R
T
=
1. So the largest compact invariant set in Γ such that V
/
d
t
≤
0 when R
T
≤
1 is the singleton E
0. Therefore, by the LaSalle Invariance Principle [26], the infection free equilibrium point E
0 is globally asymptotically stable if R
T
<
1 and unstable otherwise.


#### 3.3.2. Stability ofE
1

Theorem 4.
The endemic steady stateE
1 of the subsystem (14) is locally asymptotically stable if R
T
>
1.Proof.
The Jacobian matrix of system (14) at the equilibrium point E
1 is given by
(27)
J
E
1
=
(
-
μ
W
β
3
α
~
-
μ
d
)
.
Given that the trace of J
E
1 is negative and the determinant is negative if R
T
>
1, we can thus conclude that the unique endemic equilibrium is locally asymptotically stable whenever R
T
>
1.Theorem 5.
IfR
T
>
1, then the unique endemic equilibrium E
1 is globally stable in the interior of Γ.Proof.
We now prove the global stability of endemic steady stateE
1 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 (
F
(
i
w
,
x
)
,
G
(
i
w
,
x
)
), we can construct a Dulac function
(28)
B
(
i
w
,
x
)
=
1
β
3
i
w
x
,
i
w
>
0
,
x
>
0
.
We will thus have
(29)
∂
(
F
B
)
∂
i
w
+
∂
(
G
B
)
∂
x
=
-
(
1
i
w
2
+
α
~
β
3
x
2
)
<
0
.
Thus, subsystem (14) does not have a limit cycle in Γ. From Theorem 4, if R
T
>
1, then E
1 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 andM. ulcerans impact the human population. The overall goal is to mitigate the influence of theM. ulceranson the human population. We can actually evaluate the force of infection so that
(30)
Λ
~
=
(
R
T
-
1
)
μ
W
(
m
1
β
1
μ
d
+
α
~
β
2
K
~
+
α
~
(
R
T
-
1
)
μ
W
)
.
This means that the analysis of submodel (13) is subject to R
T
>
1. 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 R
T.Figure 2
The plot of the force of infection as a function ofR
T. The force of infection increases linearly with the reproduction number. The human population is at risk only if R
T
>
1.Using the second equation of system (13), we can evaluate s
h
* so that
(31)
s
h
*
=
(
[
σ
i
h
*
+
μ
h
i
h
(
1
+
N
H
i
h
*
)
]
[
K
~
+
α
~
μ
W
(
R
T
-
1
)
]
)
×
(
(
1
+
N
H
i
h
*
)
[
m
1
β
1
μ
d
μ
W
(
R
T
-
1
)
m
m
m
m
m
m
m
m
m
m
m
×
{
K
~
+
α
~
μ
W
(
R
T
-
1
)
}
m
m
m
m
m
m
m
m
m
m
m
+
α
~
β
2
μ
W
(
R
T
-
1
)
]
(
1
+
N
H
i
h
*
)
[
m
1
β
1
μ
d
μ
W
(
R
T
-
1
)
)
-
1
.
From the third equation of (13), we have
(32)
τ
h
*
=
σ
i
h
*
(
1
+
N
H
i
h
*
)
(
γ
+
μ
H
)
.Substitutings
h
* and τ
h
* in the first equation of (13) at the steady state yields a quadratic equation in i
h
* given by
(33)
a
i
h
*
2
+
b
i
h
*
+
c
=
0
,
where
(34)
a
=
N
H
(
μ
H
+
γ
)
(
μ
H
+
θ
)
×
(
α
~
β
2
μ
W
(
R
T
-
1
)
+
(
K
~
+
α
~
(
R
T
-
1
)
μ
W
)
mmm
×
[
μ
H
+
m
1
β
1
μ
d
μ
W
(
R
T
-
1
)
]
α
~
β
2
μ
W
(
R
T
-
1
)
+
(
K
~
+
α
~
(
R
T
-
1
)
μ
W
)
)
,
b
=
K
~
γ
θ
σ
+
K
~
μ
H
[
θ
σ
+
γ
(
θ
+
σ
)
mmmimmmmm
+
μ
H
(
μ
H
+
γ
+
θ
+
σ
)
]
+
(
R
p
-
1
)
[
α
~
(
μ
H
+
γ
)
(
μ
H
+
θ
)
(
μ
H
+
σ
)
mmmmmimm
+
α
~
β
2
(
γ
θ
+
(
γ
+
θ
)
σ
-
N
H
(
μ
H
+
γ
)
miimmmmmmimmm
×

(
μ
H
+
θ
)
+
μ
H

(
μ
H
+
γ
+
θ
+
σ
)
)
mmmmmmim
+
K
~
m
1
β
1
μ
d
{
γ
θ
+

(
γ
+
θ
)
σ
-
N
H

(
μ
H
+
γ
)
mmmmmmmmmmimmmmi
×
(
μ
H
+
θ
)
+
μ
H
mmmmimimmmmmmmmm
×

(
μ
H
+
γ
+
θ
+
σ
)
}
]
μ
W
-
α
~
m
1
(
R
p
-
1
)
2
β
1
μ
d
(
-
(
θ
σ
+
γ
(
θ
+
σ
)
)
mmmmmmmmmmmmm
+
N
H
(
γ
+
μ
H
)
(
θ
+
μ
H
)
mmmmmmmmmmmmm
-
μ
H
(
γ
+
θ
+
σ
+
μ
H
)
)
μ
W
2
,
c
=
-
μ
W
(
μ
H
+
γ
)
(
μ
H
+
θ
)
×
[
α
~
β
2
+
m
1
β
1
μ
d
(
K
~
+
α
~
μ
W
(
R
T
-
1
)
)
]
(
R
T
-
1
)
.
Clearly our model has two possible steady states given by
(35)
E
2
a
=
(
s
h
*
,
i
h
*
+
,
τ
h
*
)
,
E
2
b
=
(
s
h
*
,
i
h
*
-
,
τ
h
*
)
,
where i
h
*
± are roots of the quadratic equation (33). We note that if R
T
>
1, we have a
>
0 and c
<
0. By Descartes’ rule of signs, irrespective of the sign of b, the quadratic equation (33) has one positive root; the endemic equilibrium E
2
a
=
E
2. We thus have the following result.Theorem 6.
System (13) has a unique endemic equilibrium E
2 whenever R
T
>
1.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 E
0
h
=
(
1,0
,
0
). The existence of E
0
h is thus subject to the water bugs and the environment being free ofM. ulcerans.


#### 3.4.1. Local Stability ofE
0
h

Theorem 8.
The disease free equilibriumE
0
h whenever it exists is locally asymptotically stable if R
T
<
1 and unstable otherwise.Proof.
WhenR
T
<
1, then either there are no infections in the water bugs or they are simply carriers. So E
0
h exists. The Jacobian matrix of system (13) at the disease free equilibrium point E
0
h is given by
(36)
J
E
0
h
=
(
-
(
μ
H
+
θ
)
-
θ
-
θ
0
-
(
μ
H
+
σ
)
0
0
σ
-
(
μ
H
+
γ
)
)
.
The eigenvalues of J
E
0
h are λ
1
=
-
(
μ
H
+
θ
), λ
2
=
-
(
μ
H
+
σ
), and λ
3
=
-
(
μ
H
+
γ
). We can thus conclude that the disease free equilibrium is locally asymptotically stable whenever R
T
<
1.


#### 3.4.2. Local Stability ofE
2

Theorem 9.
The unique endemic equilibrium pointE
2 is locally asymptotically stable for R
T
>
1.Proof.
The Jacobian matrix at the endemic steady stateE
2 is given by
(37)
J
E
2
=
(
-
(
μ
H
+
θ
)
-
Λ
~
-
θ
-
θ
Λ
~
-
μ
H
-
σ
(
1
+
i
h
*
)
2
0
0
σ
(
1
+
i
h
*
)
2
-
(
μ
H
+
γ
)
)
.
If we let ψ
=
σ
/
(
1
+
i
h
*
)
2, then the eigenvalues of J
E
2 are given by the solutions of the characteristic polynomial
(38)
ϑ
3
+
η
1
ϑ
2
+
η
2
ϑ
+
η
3
=
0
,
where
(39)
η
1
=
(
μ
H
+
γ
)
+
(
μ
H
+
θ
)
+
(
μ
H
+
ψ
)
+
Λ
~
,
η
2
=
(
μ
H
+
θ
)
(
μ
H
+
γ
Λ
~
)
+
(
μ
H
+
γ
)
(
μ
H
+
ψ
+
Λ
~
)
+
(
μ
H
+
θ
)
(
μ
H
+
ψ
)
+
Λ
~
ψ
,
η
3
=
θ
Λ
~
(
γ
+
ψ
)
+
γ
(
θ
+
Λ
~
)
ψ
mmm
+
μ
H
(
+
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
Λ
~
ψ
+
θ
(
Λ
~
+
ψ
)
+
γ
(
θ
+
Λ
~
+
ψ
)
mmmmmim
+
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
)
.
Using the Routh-Hurwitz criterion, we note thatη
1
>
0,  η
2
>
0 and η
3
>
0. The evaluation of η
1
η
2
-
η
3 yields
(40)
(
θ
+
Λ
~
)
(
θ
+
ψ
)
(
Λ
~
+
ψ
)
+
γ
2
(
θ
+
Λ
~
+
ψ
)
+
γ
(
θ
+
Λ
~
+
ψ
)
2
+
2
μ
H
(
+
4
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
γ
2
+
θ
2
+
Λ
~
2
+
3
Λ
~
ψ
+
ψ
2
+
3
θ
(
Λ
~
+
ψ
)
mmmmmm
+
3
γ
(
θ
+
Λ
~
+
ψ
)
+
4
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
)
>
0
.
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 E
2 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.
IfR
T
>
1, system (13) is uniformly persistent in Γ
^, the interior of Γ.
The existence ofE
0
h, only if R
T
>
1, guarantees uniform persistence [31]. System (13) is said to be uniformly persistent if there exists a positive constant c such that any solution (
s
h
(
t
)
,
i
h
(
t
)
,
τ
h
(
t
)
) with initial conditions (
s
h
(
0
)
,
i
h
(
0
)
,
τ
h
(
0
)
)
∈
Γ
^ satisfies
(41)
lim
⁡
inf
⁡
t
→
∞
s
h
(
t
)
>
c
,
lim
⁡
inf
⁡
t
→
∞
i
h
(
t
)
>
c
,
lim
⁡
inf
⁡
t
→
∞
τ
h
(
t
)
>
c
.The proof of uniform persistence can be done using uniform persistence results in [31, 32].Theorem 11.
IfΛ
~
>
γ, the endemic equilibrium point E
2 of system (13) is globally asymptotically stable when R
T
>
1.Proof.
Using the arguments in [28], system (13) satisfies assumptions H
(
1
) and H
(
2
) in Γ
^. Let x
=
(
s
h
,
i
h
,
τ
h
) and f
(
x
) be the vector field of system (13). The Jacobian matrix corresponding to system (13) is
(42)
J
(
s
h
,
i
h
,
τ
h
)
=
(
-
(
θ
+
Λ
~
+
μ
H
)
-
θ
-
θ
Λ
~
-
(
σ
(
1
+
N
H
i
h
)
2
+
μ
H
)
0
0
σ
(
1
+
N
H
i
h
)
2
-
(
μ
H
+
γ
)
)
.
The second additive compound matrix J
(
s
h
,
i
h
,
τ
h
)
[
2
] is given by(43)
J
(
s
h
,
i
h
,
τ
h
)
[2]
=
(
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
0
θ
σ
(
1
+
N
H
i
h
)
2
-
(
θ
+
Λ
~
+
2
μ
H
+
γ
)
0
0
Λ
~
-
(
2
μ
H
+
γ
+
σ
(
1
+
N
H
i
h
)
2
)
)
.We let the matrix function P take the form
(44)
P
(
s
h
,
i
h
,
τ
h
)
=
diag
⁡
{
i
h
τ
h
,
i
h
τ
h
,
i
h
τ
h
}
.
We thus have
(45)
P
f
P
-
1
=
diag
⁡
{
i
h
′
i
h
-
τ
h
′
τ
h
,
i
h
′
i
h
-
τ
h
′
τ
h
,
i
h
′
i
h
-
τ
h
′
τ
h
}
,
where P
f is the diagonal element matrix derivative of P with respect to time and(46)
P
J
[
2
]
P
-
1
=
(
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
0
θ
σ
(
1
+
N
H
i
h
)
2
-
(
θ
+
Λ
~
+
2
μ
H
+
γ
)
-
θ
0
Λ
~
-
(
2
μ
H
+
γ
+
σ
(
1
+
N
H
i
h
)
2
)
)
,where ′ represents the derivative with respect to time.
The matrixQ
=
P
f
P
-
1
+
P
J
[
2
]
P
-
1 can be written as a block matrix so that
(47)
Q
=
(
Q
11
Q
12
Q
21
Q
22
)
,
where(48)
Q
11
=
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
+
i
h
′
i
h
-
τ
h
′
τ
h
,
Q
12
=
(
0
θ
)
,
Q
21
=
(
σ
(
1
+
N
H
i
h
)
2
0
)
,
Q
22
=
(
-
(
θ
+
Λ
~
+
2
μ
H
+
γ
)
+
i
h
′
i
h
-
τ
h
′
τ
h
-
θ
Λ
~
-
(
2
μ
H
+
γ
+
σ
(
1
+
N
H
i
h
)
2
)
+
i
h
′
i
h
-
τ
h
′
τ
h
)
.
Let(
x
,
y
,
z
) denote the vectors in R
3 and let the norm in R
3 be defined by
(49)
|
(
x
,
y
,
z
)
|
=
max
⁡
{
|
x
|
,
|
y
+
z
|
}
.
Also letL denote the Lozinskii
ˇ measure with respect to this norm. Following [33] we have
(50)
L
(
Q
)
≤
sup
⁡
{
g
1
,
g
2
}
≡
sup
⁡
{
L
1
(
Q
11
)
+
|
Q
12
|
,
L
1
(
Q
22
)
+
|
Q
21
|
}
,
where |
Q
12
| and |
Q
21
| are the matrix norms with respect to the vector norm L
1 and L
1 is the Lozinskii
ˇ measure with respect to the L
1 norm.
In fact(51)
L
1
(
Q
11
)
=
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
+
i
h
′
i
h
-
τ
h
′
τ
h
,
|
Q
12
|
=
θ
,
|
Q
21
|
=
σ
(
1
+
N
H
i
h
)
2
,
L
1
(
Q
22
)
=
-
(
θ
+
2
μ
H
+
γ
)
+
i
h
′
i
h
-
τ
h
′
τ
h
.
We now have
(52)
g
1
=
i
h
′
i
h
-
[
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
-
τ
h
′
τ
h
,
g
2
=
i
h
′
i
h
-
(
θ
+
2
μ
H
+
γ
)
+
σ
(
1
+
N
H
i
h
)
2
-
τ
h
′
τ
h
.
The third equation of (13) gives
(53)
τ
h
′
τ
h
=
(
σ
(
1
+
N
H
i
h
)
)
(
i
h
τ
h
)
-
(
μ
H
+
γ
)
.
Substituting (53) into (52) yields
(54)
g
1
=
i
h
′
i
h
-
[
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
-
τ
h
′
τ
h
=
i
h
′
i
h
-
[
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
-
{
(
σ
(
1
+
N
H
i
h
)
)
(
i
h
τ
h
)
-
(
μ
H
+
γ
)
}
≤
i
h
′
i
h
-
{
Λ
~
-
γ
+
μ
H
+
σ
(
1
+
N
H
c
)
2
+
(
σ
(
1
+
N
H
c
)
2
)
}
,
g
2
=
i
h
′
i
h
-
(
θ
+
2
μ
H
+
γ
)
+
σ
(
1
+
N
H
i
h
)
2
-
τ
h
′
τ
h
=
i
h
′
i
h
-
(
θ
+
μ
H
)
+
(
σ
(
1
+
N
H
i
h
)
)
(
1
(
1
+
N
H
i
h
)
-
i
h
τ
h
)
≤
i
h
′
i
h
-
[
(
θ
+
μ
H
)
-
(
σ
(
1
+
N
H
c
)
)
(
1
(
1
+
N
H
c
)
-
1
)
]
,
where c is the constant of uniform persistence. The inequalities follow Theorem 10.
If we impose the conditionΛ
~
>
γ, then
(55)
L
(
Q
)
≤
sup
⁡
{
g
1
,
g
2
}
=
i
h
′
i
h
-
ω
,
where ω
=
min
⁡
{
ω
1
,
ω
2
} with
(56)
ω
1
=
Λ
~
-
γ
+
μ
H
+
σ
(
1
+
N
H
c
)
2
+
(
σ
(
1
+
N
H
c
)
2
)
,
ω
2
=
(
θ
+
μ
H
)
-
(
σ
(
1
+
N
H
c
)
)
(
1
(
1
+
N
H
c
)
-
1
)
.
Hence
(57)
1
t
∫
0
t
L
(
Q
)
d
s
≤
1
t
log
⁡
i
h
(
t
)
i
h
(
0
)
-
ω
.
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].


## 3.1. Basic Properties

Since the model monitors changes in the populations of humans and water bugs and the density ofM. ulcerans in the environment, the model parameters and variables are nonnegative. The biologically feasible region for the systems (13)-(14) is in R
+
5 and is represented by the set
(15)
Γ
=
{
0
≤
i
w
≤
1
,
0
≤
x
≤
α
~
μ
d
(
s
h
,
i
h
,
τ
h
,
i
w
,
x
)
∈
R
+
5
∣
0
≤
s
h
+
i
h
+
τ
h
≤
1
,
mim
0
≤
i
w
≤
1
,
0
≤
x
≤
α
~
μ
d
}
,
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 d
x
/
d
t
=
α
~
i
w
-
μ
d
x
≤
α
~
-
μ
d
x, we have x
≤
α
~
/
μ
d. The solutions of systems (13)-(14) starting in Γ remain in Γ for all t
>
0. 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 t
∈
[
0
,
∞
). 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 t
>
0. We have the following proposition.Proposition 1.
For positive initial conditions of systems (13)-(14), the solutions s
h
(
t
),  i
h
(
t
),  τ
h
(
t
),  i
w
(
t
), and x
(
t
) are nonnegative for all t
>
0.Proof.
Assume that(16)
t
^
=
sup
⁡
{
t
>
0
:
s
h
>
0
,
i
h
>
0
,
τ
h
>
0
,
i
w
>
0
,
x
>
0
}
∈
(
0
,
t
]
.
Thus t
^
>
0, and it follows directly from the first equation of the subsystem (13) that
(17)
d
s
h
d
t
≤
(
μ
H
+
θ
)
-
[
(
μ
H
+
θ
)
+
Λ
]
s
h
.
This is a first order differential equation that can easily be solved using an integrating factor. For a nonconstant force of infection Λ, we have
(18)
s
h
(
t
^
)
≤
s
h
(
0
)
exp
⁡
[
-
(
(
μ
H
+
θ
)
t
^
+
∫
0
t
^
‍
Λ
(
s
)
d
s
)
]
+
exp
⁡
[
-
(
(
μ
H
+
θ
)
t
^
+
∫
0
t
^
‍
Λ
(
s
)
d
s
)
]
×
[
∫
0
t
^
‍
(
μ
H
+
θ
)
e
(
(
μ
H
+
θ
)
t
^
+
∫
0
t
^
‍
Λ
(
l
)
d
l
)
d
t
^
]
.
Since the right-hand side of (18) is always positive, the solution s
h
(
t
) will always be positive. If Λ is constant, this result still holds.
From the second equation of subsystem (13),
(19)
d
i
h
d
t
≥
-
(
μ
H
+
σ
)
i
h
≥
i
h
(
0
)
exp
⁡
[
-
(
μ
H
+
σ
)
t
]
>
0
.
The third equation of subsystem (13) yields
(20)
d
τ
h
d
t
≥
-
(
μ
H
+
γ
)
τ
h
≥
τ
h
(
0
)
exp
⁡
[
-
(
μ
H
+
γ
)
t
]
>
0
.
Similarly, we can show that i
w
(
t
)
>
0 and x
(
t
)
>
0 for all t
>
0 and this completes the proof.


## 3.3. Environmental Dynamics

The subsystem (14) represents the dynamics of water bugs andM. ulcerans in the environment. From the second equation, we have
(21)
x
*
=
α
~
i
w
*
μ
d
,
i
w
*
=
0
or
i
w
*
=
1
-
1
R
T
,
where
(22)
R
T
=
α
~
β
3
μ
d
μ
W
.
In this case x
*
=
(
α
~
/
μ
d
)
(
1
-
(
1
/
R
T
)
).The casei
w
*
=
0 yields the infection free equilibrium point of the environmental dynamics submodel given by
(23)
E
0
=
(
0,0
)
.
The submodel also has an endemic equilibrium given by
(24)
E
1
=
(
α
~
μ
W
(
R
T
-
1
)
,
μ
d
μ
W
(
R
T
-
1
)
)
.Remark 2.
It is important to note that theR
T is our model reproduction number for the BU epidemic in the presence of treatment driven by the dynamics of the water bug andM. 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 sheddedM. ulcerans in the environment. R
T 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 andM. 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 theM. ulcerans into the environment and the effective contact rate between the water bugs andM. ulcerans. This implies that the control and management of the ulcer largely depend on environmental management.


### 3.3.1. Stability ofE
0

Theorem 3.
The infection free equilibriumE
0 is globally stable when R
T
<
1 and unstable otherwise.Proof.
We propose a Lyapunov function of the form(25)
V
(
t
)
=
i
w
+
β
3
μ
d
x
.
The time derivative of (25) is
(26)
V
˙
=
d
i
w
d
t
+
β
3
μ
d
d
x
d
t
≤
μ
W
(
R
T
-
1
)
i
w
.
When R
T
≤
1,  V
˙ is negative and semidefinite, with equality at the infection free equilibrium and/or at R
T
=
1. So the largest compact invariant set in Γ such that V
/
d
t
≤
0 when R
T
≤
1 is the singleton E
0. Therefore, by the LaSalle Invariance Principle [26], the infection free equilibrium point E
0 is globally asymptotically stable if R
T
<
1 and unstable otherwise.


### 3.3.2. Stability ofE
1

Theorem 4.
The endemic steady stateE
1 of the subsystem (14) is locally asymptotically stable if R
T
>
1.Proof.
The Jacobian matrix of system (14) at the equilibrium point E
1 is given by
(27)
J
E
1
=
(
-
μ
W
β
3
α
~
-
μ
d
)
.
Given that the trace of J
E
1 is negative and the determinant is negative if R
T
>
1, we can thus conclude that the unique endemic equilibrium is locally asymptotically stable whenever R
T
>
1.Theorem 5.
IfR
T
>
1, then the unique endemic equilibrium E
1 is globally stable in the interior of Γ.Proof.
We now prove the global stability of endemic steady stateE
1 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 (
F
(
i
w
,
x
)
,
G
(
i
w
,
x
)
), we can construct a Dulac function
(28)
B
(
i
w
,
x
)
=
1
β
3
i
w
x
,
i
w
>
0
,
x
>
0
.
We will thus have
(29)
∂
(
F
B
)
∂
i
w
+
∂
(
G
B
)
∂
x
=
-
(
1
i
w
2
+
α
~
β
3
x
2
)
<
0
.
Thus, subsystem (14) does not have a limit cycle in Γ. From Theorem 4, if R
T
>
1, then E
1 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.3.1. Stability ofE
0

Theorem 3.
The infection free equilibriumE
0 is globally stable when R
T
<
1 and unstable otherwise.Proof.
We propose a Lyapunov function of the form(25)
V
(
t
)
=
i
w
+
β
3
μ
d
x
.
The time derivative of (25) is
(26)
V
˙
=
d
i
w
d
t
+
β
3
μ
d
d
x
d
t
≤
μ
W
(
R
T
-
1
)
i
w
.
When R
T
≤
1,  V
˙ is negative and semidefinite, with equality at the infection free equilibrium and/or at R
T
=
1. So the largest compact invariant set in Γ such that V
/
d
t
≤
0 when R
T
≤
1 is the singleton E
0. Therefore, by the LaSalle Invariance Principle [26], the infection free equilibrium point E
0 is globally asymptotically stable if R
T
<
1 and unstable otherwise.


## 3.3.2. Stability ofE
1

Theorem 4.
The endemic steady stateE
1 of the subsystem (14) is locally asymptotically stable if R
T
>
1.Proof.
The Jacobian matrix of system (14) at the equilibrium point E
1 is given by
(27)
J
E
1
=
(
-
μ
W
β
3
α
~
-
μ
d
)
.
Given that the trace of J
E
1 is negative and the determinant is negative if R
T
>
1, we can thus conclude that the unique endemic equilibrium is locally asymptotically stable whenever R
T
>
1.Theorem 5.
IfR
T
>
1, then the unique endemic equilibrium E
1 is globally stable in the interior of Γ.Proof.
We now prove the global stability of endemic steady stateE
1 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 (
F
(
i
w
,
x
)
,
G
(
i
w
,
x
)
), we can construct a Dulac function
(28)
B
(
i
w
,
x
)
=
1
β
3
i
w
x
,
i
w
>
0
,
x
>
0
.
We will thus have
(29)
∂
(
F
B
)
∂
i
w
+
∂
(
G
B
)
∂
x
=
-
(
1
i
w
2
+
α
~
β
3
x
2
)
<
0
.
Thus, subsystem (14) does not have a limit cycle in Γ. From Theorem 4, if R
T
>
1, then E
1 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 andM. ulcerans impact the human population. The overall goal is to mitigate the influence of theM. ulceranson the human population. We can actually evaluate the force of infection so that
(30)
Λ
~
=
(
R
T
-
1
)
μ
W
(
m
1
β
1
μ
d
+
α
~
β
2
K
~
+
α
~
(
R
T
-
1
)
μ
W
)
.
This means that the analysis of submodel (13) is subject to R
T
>
1. 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 R
T.Figure 2
The plot of the force of infection as a function ofR
T. The force of infection increases linearly with the reproduction number. The human population is at risk only if R
T
>
1.Using the second equation of system (13), we can evaluate s
h
* so that
(31)
s
h
*
=
(
[
σ
i
h
*
+
μ
h
i
h
(
1
+
N
H
i
h
*
)
]
[
K
~
+
α
~
μ
W
(
R
T
-
1
)
]
)
×
(
(
1
+
N
H
i
h
*
)
[
m
1
β
1
μ
d
μ
W
(
R
T
-
1
)
m
m
m
m
m
m
m
m
m
m
m
×
{
K
~
+
α
~
μ
W
(
R
T
-
1
)
}
m
m
m
m
m
m
m
m
m
m
m
+
α
~
β
2
μ
W
(
R
T
-
1
)
]
(
1
+
N
H
i
h
*
)
[
m
1
β
1
μ
d
μ
W
(
R
T
-
1
)
)
-
1
.
From the third equation of (13), we have
(32)
τ
h
*
=
σ
i
h
*
(
1
+
N
H
i
h
*
)
(
γ
+
μ
H
)
.Substitutings
h
* and τ
h
* in the first equation of (13) at the steady state yields a quadratic equation in i
h
* given by
(33)
a
i
h
*
2
+
b
i
h
*
+
c
=
0
,
where
(34)
a
=
N
H
(
μ
H
+
γ
)
(
μ
H
+
θ
)
×
(
α
~
β
2
μ
W
(
R
T
-
1
)
+
(
K
~
+
α
~
(
R
T
-
1
)
μ
W
)
mmm
×
[
μ
H
+
m
1
β
1
μ
d
μ
W
(
R
T
-
1
)
]
α
~
β
2
μ
W
(
R
T
-
1
)
+
(
K
~
+
α
~
(
R
T
-
1
)
μ
W
)
)
,
b
=
K
~
γ
θ
σ
+
K
~
μ
H
[
θ
σ
+
γ
(
θ
+
σ
)
mmmimmmmm
+
μ
H
(
μ
H
+
γ
+
θ
+
σ
)
]
+
(
R
p
-
1
)
[
α
~
(
μ
H
+
γ
)
(
μ
H
+
θ
)
(
μ
H
+
σ
)
mmmmmimm
+
α
~
β
2
(
γ
θ
+
(
γ
+
θ
)
σ
-
N
H
(
μ
H
+
γ
)
miimmmmmmimmm
×

(
μ
H
+
θ
)
+
μ
H

(
μ
H
+
γ
+
θ
+
σ
)
)
mmmmmmim
+
K
~
m
1
β
1
μ
d
{
γ
θ
+

(
γ
+
θ
)
σ
-
N
H

(
μ
H
+
γ
)
mmmmmmmmmmimmmmi
×
(
μ
H
+
θ
)
+
μ
H
mmmmimimmmmmmmmm
×

(
μ
H
+
γ
+
θ
+
σ
)
}
]
μ
W
-
α
~
m
1
(
R
p
-
1
)
2
β
1
μ
d
(
-
(
θ
σ
+
γ
(
θ
+
σ
)
)
mmmmmmmmmmmmm
+
N
H
(
γ
+
μ
H
)
(
θ
+
μ
H
)
mmmmmmmmmmmmm
-
μ
H
(
γ
+
θ
+
σ
+
μ
H
)
)
μ
W
2
,
c
=
-
μ
W
(
μ
H
+
γ
)
(
μ
H
+
θ
)
×
[
α
~
β
2
+
m
1
β
1
μ
d
(
K
~
+
α
~
μ
W
(
R
T
-
1
)
)
]
(
R
T
-
1
)
.
Clearly our model has two possible steady states given by
(35)
E
2
a
=
(
s
h
*
,
i
h
*
+
,
τ
h
*
)
,
E
2
b
=
(
s
h
*
,
i
h
*
-
,
τ
h
*
)
,
where i
h
*
± are roots of the quadratic equation (33). We note that if R
T
>
1, we have a
>
0 and c
<
0. By Descartes’ rule of signs, irrespective of the sign of b, the quadratic equation (33) has one positive root; the endemic equilibrium E
2
a
=
E
2. We thus have the following result.Theorem 6.
System (13) has a unique endemic equilibrium E
2 whenever R
T
>
1.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 E
0
h
=
(
1,0
,
0
). The existence of E
0
h is thus subject to the water bugs and the environment being free ofM. ulcerans.


### 3.4.1. Local Stability ofE
0
h

Theorem 8.
The disease free equilibriumE
0
h whenever it exists is locally asymptotically stable if R
T
<
1 and unstable otherwise.Proof.
WhenR
T
<
1, then either there are no infections in the water bugs or they are simply carriers. So E
0
h exists. The Jacobian matrix of system (13) at the disease free equilibrium point E
0
h is given by
(36)
J
E
0
h
=
(
-
(
μ
H
+
θ
)
-
θ
-
θ
0
-
(
μ
H
+
σ
)
0
0
σ
-
(
μ
H
+
γ
)
)
.
The eigenvalues of J
E
0
h are λ
1
=
-
(
μ
H
+
θ
), λ
2
=
-
(
μ
H
+
σ
), and λ
3
=
-
(
μ
H
+
γ
). We can thus conclude that the disease free equilibrium is locally asymptotically stable whenever R
T
<
1.


### 3.4.2. Local Stability ofE
2

Theorem 9.
The unique endemic equilibrium pointE
2 is locally asymptotically stable for R
T
>
1.Proof.
The Jacobian matrix at the endemic steady stateE
2 is given by
(37)
J
E
2
=
(
-
(
μ
H
+
θ
)
-
Λ
~
-
θ
-
θ
Λ
~
-
μ
H
-
σ
(
1
+
i
h
*
)
2
0
0
σ
(
1
+
i
h
*
)
2
-
(
μ
H
+
γ
)
)
.
If we let ψ
=
σ
/
(
1
+
i
h
*
)
2, then the eigenvalues of J
E
2 are given by the solutions of the characteristic polynomial
(38)
ϑ
3
+
η
1
ϑ
2
+
η
2
ϑ
+
η
3
=
0
,
where
(39)
η
1
=
(
μ
H
+
γ
)
+
(
μ
H
+
θ
)
+
(
μ
H
+
ψ
)
+
Λ
~
,
η
2
=
(
μ
H
+
θ
)
(
μ
H
+
γ
Λ
~
)
+
(
μ
H
+
γ
)
(
μ
H
+
ψ
+
Λ
~
)
+
(
μ
H
+
θ
)
(
μ
H
+
ψ
)
+
Λ
~
ψ
,
η
3
=
θ
Λ
~
(
γ
+
ψ
)
+
γ
(
θ
+
Λ
~
)
ψ
mmm
+
μ
H
(
+
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
Λ
~
ψ
+
θ
(
Λ
~
+
ψ
)
+
γ
(
θ
+
Λ
~
+
ψ
)
mmmmmim
+
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
)
.
Using the Routh-Hurwitz criterion, we note thatη
1
>
0,  η
2
>
0 and η
3
>
0. The evaluation of η
1
η
2
-
η
3 yields
(40)
(
θ
+
Λ
~
)
(
θ
+
ψ
)
(
Λ
~
+
ψ
)
+
γ
2
(
θ
+
Λ
~
+
ψ
)
+
γ
(
θ
+
Λ
~
+
ψ
)
2
+
2
μ
H
(
+
4
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
γ
2
+
θ
2
+
Λ
~
2
+
3
Λ
~
ψ
+
ψ
2
+
3
θ
(
Λ
~
+
ψ
)
mmmmmm
+
3
γ
(
θ
+
Λ
~
+
ψ
)
+
4
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
)
>
0
.
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 E
2 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.
IfR
T
>
1, system (13) is uniformly persistent in Γ
^, the interior of Γ.
The existence ofE
0
h, only if R
T
>
1, guarantees uniform persistence [31]. System (13) is said to be uniformly persistent if there exists a positive constant c such that any solution (
s
h
(
t
)
,
i
h
(
t
)
,
τ
h
(
t
)
) with initial conditions (
s
h
(
0
)
,
i
h
(
0
)
,
τ
h
(
0
)
)
∈
Γ
^ satisfies
(41)
lim
⁡
inf
⁡
t
→
∞
s
h
(
t
)
>
c
,
lim
⁡
inf
⁡
t
→
∞
i
h
(
t
)
>
c
,
lim
⁡
inf
⁡
t
→
∞
τ
h
(
t
)
>
c
.The proof of uniform persistence can be done using uniform persistence results in [31, 32].Theorem 11.
IfΛ
~
>
γ, the endemic equilibrium point E
2 of system (13) is globally asymptotically stable when R
T
>
1.Proof.
Using the arguments in [28], system (13) satisfies assumptions H
(
1
) and H
(
2
) in Γ
^. Let x
=
(
s
h
,
i
h
,
τ
h
) and f
(
x
) be the vector field of system (13). The Jacobian matrix corresponding to system (13) is
(42)
J
(
s
h
,
i
h
,
τ
h
)
=
(
-
(
θ
+
Λ
~
+
μ
H
)
-
θ
-
θ
Λ
~
-
(
σ
(
1
+
N
H
i
h
)
2
+
μ
H
)
0
0
σ
(
1
+
N
H
i
h
)
2
-
(
μ
H
+
γ
)
)
.
The second additive compound matrix J
(
s
h
,
i
h
,
τ
h
)
[
2
] is given by(43)
J
(
s
h
,
i
h
,
τ
h
)
[2]
=
(
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
0
θ
σ
(
1
+
N
H
i
h
)
2
-
(
θ
+
Λ
~
+
2
μ
H
+
γ
)
0
0
Λ
~
-
(
2
μ
H
+
γ
+
σ
(
1
+
N
H
i
h
)
2
)
)
.We let the matrix function P take the form
(44)
P
(
s
h
,
i
h
,
τ
h
)
=
diag
⁡
{
i
h
τ
h
,
i
h
τ
h
,
i
h
τ
h
}
.
We thus have
(45)
P
f
P
-
1
=
diag
⁡
{
i
h
′
i
h
-
τ
h
′
τ
h
,
i
h
′
i
h
-
τ
h
′
τ
h
,
i
h
′
i
h
-
τ
h
′
τ
h
}
,
where P
f is the diagonal element matrix derivative of P with respect to time and(46)
P
J
[
2
]
P
-
1
=
(
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
0
θ
σ
(
1
+
N
H
i
h
)
2
-
(
θ
+
Λ
~
+
2
μ
H
+
γ
)
-
θ
0
Λ
~
-
(
2
μ
H
+
γ
+
σ
(
1
+
N
H
i
h
)
2
)
)
,where ′ represents the derivative with respect to time.
The matrixQ
=
P
f
P
-
1
+
P
J
[
2
]
P
-
1 can be written as a block matrix so that
(47)
Q
=
(
Q
11
Q
12
Q
21
Q
22
)
,
where(48)
Q
11
=
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
+
i
h
′
i
h
-
τ
h
′
τ
h
,
Q
12
=
(
0
θ
)
,
Q
21
=
(
σ
(
1
+
N
H
i
h
)
2
0
)
,
Q
22
=
(
-
(
θ
+
Λ
~
+
2
μ
H
+
γ
)
+
i
h
′
i
h
-
τ
h
′
τ
h
-
θ
Λ
~
-
(
2
μ
H
+
γ
+
σ
(
1
+
N
H
i
h
)
2
)
+
i
h
′
i
h
-
τ
h
′
τ
h
)
.
Let(
x
,
y
,
z
) denote the vectors in R
3 and let the norm in R
3 be defined by
(49)
|
(
x
,
y
,
z
)
|
=
max
⁡
{
|
x
|
,
|
y
+
z
|
}
.
Also letL denote the Lozinskii
ˇ measure with respect to this norm. Following [33] we have
(50)
L
(
Q
)
≤
sup
⁡
{
g
1
,
g
2
}
≡
sup
⁡
{
L
1
(
Q
11
)
+
|
Q
12
|
,
L
1
(
Q
22
)
+
|
Q
21
|
}
,
where |
Q
12
| and |
Q
21
| are the matrix norms with respect to the vector norm L
1 and L
1 is the Lozinskii
ˇ measure with respect to the L
1 norm.
In fact(51)
L
1
(
Q
11
)
=
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
+
i
h
′
i
h
-
τ
h
′
τ
h
,
|
Q
12
|
=
θ
,
|
Q
21
|
=
σ
(
1
+
N
H
i
h
)
2
,
L
1
(
Q
22
)
=
-
(
θ
+
2
μ
H
+
γ
)
+
i
h
′
i
h
-
τ
h
′
τ
h
.
We now have
(52)
g
1
=
i
h
′
i
h
-
[
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
-
τ
h
′
τ
h
,
g
2
=
i
h
′
i
h
-
(
θ
+
2
μ
H
+
γ
)
+
σ
(
1
+
N
H
i
h
)
2
-
τ
h
′
τ
h
.
The third equation of (13) gives
(53)
τ
h
′
τ
h
=
(
σ
(
1
+
N
H
i
h
)
)
(
i
h
τ
h
)
-
(
μ
H
+
γ
)
.
Substituting (53) into (52) yields
(54)
g
1
=
i
h
′
i
h
-
[
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
-
τ
h
′
τ
h
=
i
h
′
i
h
-
[
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
-
{
(
σ
(
1
+
N
H
i
h
)
)
(
i
h
τ
h
)
-
(
μ
H
+
γ
)
}
≤
i
h
′
i
h
-
{
Λ
~
-
γ
+
μ
H
+
σ
(
1
+
N
H
c
)
2
+
(
σ
(
1
+
N
H
c
)
2
)
}
,
g
2
=
i
h
′
i
h
-
(
θ
+
2
μ
H
+
γ
)
+
σ
(
1
+
N
H
i
h
)
2
-
τ
h
′
τ
h
=
i
h
′
i
h
-
(
θ
+
μ
H
)
+
(
σ
(
1
+
N
H
i
h
)
)
(
1
(
1
+
N
H
i
h
)
-
i
h
τ
h
)
≤
i
h
′
i
h
-
[
(
θ
+
μ
H
)
-
(
σ
(
1
+
N
H
c
)
)
(
1
(
1
+
N
H
c
)
-
1
)
]
,
where c is the constant of uniform persistence. The inequalities follow Theorem 10.
If we impose the conditionΛ
~
>
γ, then
(55)
L
(
Q
)
≤
sup
⁡
{
g
1
,
g
2
}
=
i
h
′
i
h
-
ω
,
where ω
=
min
⁡
{
ω
1
,
ω
2
} with
(56)
ω
1
=
Λ
~
-
γ
+
μ
H
+
σ
(
1
+
N
H
c
)
2
+
(
σ
(
1
+
N
H
c
)
2
)
,
ω
2
=
(
θ
+
μ
H
)
-
(
σ
(
1
+
N
H
c
)
)
(
1
(
1
+
N
H
c
)
-
1
)
.
Hence
(57)
1
t
∫
0
t
L
(
Q
)
d
s
≤
1
t
log
⁡
i
h
(
t
)
i
h
(
0
)
-
ω
.
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].


## 3.4.1. Local Stability ofE
0
h

Theorem 8.
The disease free equilibriumE
0
h whenever it exists is locally asymptotically stable if R
T
<
1 and unstable otherwise.Proof.
WhenR
T
<
1, then either there are no infections in the water bugs or they are simply carriers. So E
0
h exists. The Jacobian matrix of system (13) at the disease free equilibrium point E
0
h is given by
(36)
J
E
0
h
=
(
-
(
μ
H
+
θ
)
-
θ
-
θ
0
-
(
μ
H
+
σ
)
0
0
σ
-
(
μ
H
+
γ
)
)
.
The eigenvalues of J
E
0
h are λ
1
=
-
(
μ
H
+
θ
), λ
2
=
-
(
μ
H
+
σ
), and λ
3
=
-
(
μ
H
+
γ
). We can thus conclude that the disease free equilibrium is locally asymptotically stable whenever R
T
<
1.


## 3.4.2. Local Stability ofE
2

Theorem 9.
The unique endemic equilibrium pointE
2 is locally asymptotically stable for R
T
>
1.Proof.
The Jacobian matrix at the endemic steady stateE
2 is given by
(37)
J
E
2
=
(
-
(
μ
H
+
θ
)
-
Λ
~
-
θ
-
θ
Λ
~
-
μ
H
-
σ
(
1
+
i
h
*
)
2
0
0
σ
(
1
+
i
h
*
)
2
-
(
μ
H
+
γ
)
)
.
If we let ψ
=
σ
/
(
1
+
i
h
*
)
2, then the eigenvalues of J
E
2 are given by the solutions of the characteristic polynomial
(38)
ϑ
3
+
η
1
ϑ
2
+
η
2
ϑ
+
η
3
=
0
,
where
(39)
η
1
=
(
μ
H
+
γ
)
+
(
μ
H
+
θ
)
+
(
μ
H
+
ψ
)
+
Λ
~
,
η
2
=
(
μ
H
+
θ
)
(
μ
H
+
γ
Λ
~
)
+
(
μ
H
+
γ
)
(
μ
H
+
ψ
+
Λ
~
)
+
(
μ
H
+
θ
)
(
μ
H
+
ψ
)
+
Λ
~
ψ
,
η
3
=
θ
Λ
~
(
γ
+
ψ
)
+
γ
(
θ
+
Λ
~
)
ψ
mmm
+
μ
H
(
+
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
Λ
~
ψ
+
θ
(
Λ
~
+
ψ
)
+
γ
(
θ
+
Λ
~
+
ψ
)
mmmmmim
+
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
)
.
Using the Routh-Hurwitz criterion, we note thatη
1
>
0,  η
2
>
0 and η
3
>
0. The evaluation of η
1
η
2
-
η
3 yields
(40)
(
θ
+
Λ
~
)
(
θ
+
ψ
)
(
Λ
~
+
ψ
)
+
γ
2
(
θ
+
Λ
~
+
ψ
)
+
γ
(
θ
+
Λ
~
+
ψ
)
2
+
2
μ
H
(
+
4
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
γ
2
+
θ
2
+
Λ
~
2
+
3
Λ
~
ψ
+
ψ
2
+
3
θ
(
Λ
~
+
ψ
)
mmmmmm
+
3
γ
(
θ
+
Λ
~
+
ψ
)
+
4
μ
H
(
γ
+
θ
+
Λ
~
+
ψ
+
μ
H
)
)
>
0
.
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 E
2 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.
IfR
T
>
1, system (13) is uniformly persistent in Γ
^, the interior of Γ.
The existence ofE
0
h, only if R
T
>
1, guarantees uniform persistence [31]. System (13) is said to be uniformly persistent if there exists a positive constant c such that any solution (
s
h
(
t
)
,
i
h
(
t
)
,
τ
h
(
t
)
) with initial conditions (
s
h
(
0
)
,
i
h
(
0
)
,
τ
h
(
0
)
)
∈
Γ
^ satisfies
(41)
lim
⁡
inf
⁡
t
→
∞
s
h
(
t
)
>
c
,
lim
⁡
inf
⁡
t
→
∞
i
h
(
t
)
>
c
,
lim
⁡
inf
⁡
t
→
∞
τ
h
(
t
)
>
c
.The proof of uniform persistence can be done using uniform persistence results in [31, 32].Theorem 11.
IfΛ
~
>
γ, the endemic equilibrium point E
2 of system (13) is globally asymptotically stable when R
T
>
1.Proof.
Using the arguments in [28], system (13) satisfies assumptions H
(
1
) and H
(
2
) in Γ
^. Let x
=
(
s
h
,
i
h
,
τ
h
) and f
(
x
) be the vector field of system (13). The Jacobian matrix corresponding to system (13) is
(42)
J
(
s
h
,
i
h
,
τ
h
)
=
(
-
(
θ
+
Λ
~
+
μ
H
)
-
θ
-
θ
Λ
~
-
(
σ
(
1
+
N
H
i
h
)
2
+
μ
H
)
0
0
σ
(
1
+
N
H
i
h
)
2
-
(
μ
H
+
γ
)
)
.
The second additive compound matrix J
(
s
h
,
i
h
,
τ
h
)
[
2
] is given by(43)
J
(
s
h
,
i
h
,
τ
h
)
[2]
=
(
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
0
θ
σ
(
1
+
N
H
i
h
)
2
-
(
θ
+
Λ
~
+
2
μ
H
+
γ
)
0
0
Λ
~
-
(
2
μ
H
+
γ
+
σ
(
1
+
N
H
i
h
)
2
)
)
.We let the matrix function P take the form
(44)
P
(
s
h
,
i
h
,
τ
h
)
=
diag
⁡
{
i
h
τ
h
,
i
h
τ
h
,
i
h
τ
h
}
.
We thus have
(45)
P
f
P
-
1
=
diag
⁡
{
i
h
′
i
h
-
τ
h
′
τ
h
,
i
h
′
i
h
-
τ
h
′
τ
h
,
i
h
′
i
h
-
τ
h
′
τ
h
}
,
where P
f is the diagonal element matrix derivative of P with respect to time and(46)
P
J
[
2
]
P
-
1
=
(
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
0
θ
σ
(
1
+
N
H
i
h
)
2
-
(
θ
+
Λ
~
+
2
μ
H
+
γ
)
-
θ
0
Λ
~
-
(
2
μ
H
+
γ
+
σ
(
1
+
N
H
i
h
)
2
)
)
,where ′ represents the derivative with respect to time.
The matrixQ
=
P
f
P
-
1
+
P
J
[
2
]
P
-
1 can be written as a block matrix so that
(47)
Q
=
(
Q
11
Q
12
Q
21
Q
22
)
,
where(48)
Q
11
=
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
+
i
h
′
i
h
-
τ
h
′
τ
h
,
Q
12
=
(
0
θ
)
,
Q
21
=
(
σ
(
1
+
N
H
i
h
)
2
0
)
,
Q
22
=
(
-
(
θ
+
Λ
~
+
2
μ
H
+
γ
)
+
i
h
′
i
h
-
τ
h
′
τ
h
-
θ
Λ
~
-
(
2
μ
H
+
γ
+
σ
(
1
+
N
H
i
h
)
2
)
+
i
h
′
i
h
-
τ
h
′
τ
h
)
.
Let(
x
,
y
,
z
) denote the vectors in R
3 and let the norm in R
3 be defined by
(49)
|
(
x
,
y
,
z
)
|
=
max
⁡
{
|
x
|
,
|
y
+
z
|
}
.
Also letL denote the Lozinskii
ˇ measure with respect to this norm. Following [33] we have
(50)
L
(
Q
)
≤
sup
⁡
{
g
1
,
g
2
}
≡
sup
⁡
{
L
1
(
Q
11
)
+
|
Q
12
|
,
L
1
(
Q
22
)
+
|
Q
21
|
}
,
where |
Q
12
| and |
Q
21
| are the matrix norms with respect to the vector norm L
1 and L
1 is the Lozinskii
ˇ measure with respect to the L
1 norm.
In fact(51)
L
1
(
Q
11
)
=
-
[
θ
+
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
+
i
h
′
i
h
-
τ
h
′
τ
h
,
|
Q
12
|
=
θ
,
|
Q
21
|
=
σ
(
1
+
N
H
i
h
)
2
,
L
1
(
Q
22
)
=
-
(
θ
+
2
μ
H
+
γ
)
+
i
h
′
i
h
-
τ
h
′
τ
h
.
We now have
(52)
g
1
=
i
h
′
i
h
-
[
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
-
τ
h
′
τ
h
,
g
2
=
i
h
′
i
h
-
(
θ
+
2
μ
H
+
γ
)
+
σ
(
1
+
N
H
i
h
)
2
-
τ
h
′
τ
h
.
The third equation of (13) gives
(53)
τ
h
′
τ
h
=
(
σ
(
1
+
N
H
i
h
)
)
(
i
h
τ
h
)
-
(
μ
H
+
γ
)
.
Substituting (53) into (52) yields
(54)
g
1
=
i
h
′
i
h
-
[
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
-
τ
h
′
τ
h
=
i
h
′
i
h
-
[
Λ
~
+
2
μ
H
+
σ
(
1
+
N
H
i
h
)
2
]
-
{
(
σ
(
1
+
N
H
i
h
)
)
(
i
h
τ
h
)
-
(
μ
H
+
γ
)
}
≤
i
h
′
i
h
-
{
Λ
~
-
γ
+
μ
H
+
σ
(
1
+
N
H
c
)
2
+
(
σ
(
1
+
N
H
c
)
2
)
}
,
g
2
=
i
h
′
i
h
-
(
θ
+
2
μ
H
+
γ
)
+
σ
(
1
+
N
H
i
h
)
2
-
τ
h
′
τ
h
=
i
h
′
i
h
-
(
θ
+
μ
H
)
+
(
σ
(
1
+
N
H
i
h
)
)
(
1
(
1
+
N
H
i
h
)
-
i
h
τ
h
)
≤
i
h
′
i
h
-
[
(
θ
+
μ
H
)
-
(
σ
(
1
+
N
H
c
)
)
(
1
(
1
+
N
H
c
)
-
1
)
]
,
where c is the constant of uniform persistence. The inequalities follow Theorem 10.
If we impose the conditionΛ
~
>
γ, then
(55)
L
(
Q
)
≤
sup
⁡
{
g
1
,
g
2
}
=
i
h
′
i
h
-
ω
,
where ω
=
min
⁡
{
ω
1
,
ω
2
} with
(56)
ω
1
=
Λ
~
-
γ
+
μ
H
+
σ
(
1
+
N
H
c
)
2
+
(
σ
(
1
+
N
H
c
)
2
)
,
ω
2
=
(
θ
+
μ
H
)
-
(
σ
(
1
+
N
H
c
)
)
(
1
(
1
+
N
H
c
)
-
1
)
.
Hence
(57)
1
t
∫
0
t
L
(
Q
)
d
s
≤
1
t
log
⁡
i
h
(
t
)
i
h
(
0
)
-
ω
.
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μ
H, 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 μ
H
=
0.0166 per year or equivalently 4.5
×
1
0
-
5 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 1.6
×
1
0
-
5 to 0.5 per day [35]. The rate of loss of immunity θ for vector borne diseases ranges between 0 and 1.1
×
1
0
-
2 per day [35]. Although the mortality rate of the water bugs is not known, it is assumed to be 0.15 per day [18]. We assume that we have more water bugs than humans so that m
1
>
1. 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 with1000 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.The scatter plots for the parametersα, β
3, μ
d, and μ
W.
(a)
(b)
(c)
(d)Figures3(a) and 3(b) depict parameters with a positive correlation with the reproduction number. They show a monotonic increase of R
T as α and β
3 increase. This means that, to curtail the epidemic, the reduction in the shedding rate and infection of water bugs byM. 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 theM. 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 Figure4. 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 ofM. ulcerans into the environment will significantly slow the epidemic.Figure 4
The tornado plots for the four parameters in the model reproduction number.


### 4.3. Simulation Results

To validate our mathematical analysis results, we plot phase diagrams forR
T 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.The phase diagrams forR
T
=
0.8889 (a) and R
T
=
5.3333 (b).
(a)
(b)Figure6 shows a three-dimensional phase diagram for the human population dynamics. The existence of the endemic equilibrium, when R
T
>
1, 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.Figure 6
A phase diagram for the human population showing the endemic steady state. For a randomly chosen set of initial conditions, all trajectories tend to an endemic equilibrium for the following parameter values:μ
H
=
0.02, θ
=
0.04, Λ
=
0.07, σ
=
0.4, and γ
=
0.7.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μ
d, 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 μ
d. 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 μ
w 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.Fraction of the infected human population forR
T
=
1.6492 for the parameters μ
d and μ
W. The parameter values used for the constant parameters are μ
H
=
0.00045,  m
1
=
10,  θ
=
0.011,  γ
=
0.000016,  β
3
=
0.09,  σ
=
0.08,  K
~
=
0.4000,  N
H
=
100000,  α
=
0.00615,  β
1
=
0.00001, and β
2
=
0.0000002.
(a)
(b)We also explore the role played by directM. 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.Figure 8
The proportion of infected humans for the given values ofβ
2 and the following parameter values: μ
H
=
0.00045;  m
1
=
10;  θ
=
0.011;  γ
=
0.000016;  β
3
=
0.09;  μ
d
=
0.02;  σ
=
0.14;  K
=
0.4000;  μ
W
=
0.025;  N
H
=
100000;  α
=
0.006; and β
1
=
0.00001.The shedding rate ofM. ulcerans into the environment and the treatment rate of humans are important to consider. In Figure 9(a) we observe that increasing the amount ofM. 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 phase diagram for the infected water bugs andM. ulcerans in the environment for the same parameters presented in Figure 7 with μ
d
=
0.02 and μ
W
=
0.025.
(a)
(b)


## 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μ
H, 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 μ
H
=
0.0166 per year or equivalently 4.5
×
1
0
-
5 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 1.6
×
1
0
-
5 to 0.5 per day [35]. The rate of loss of immunity θ for vector borne diseases ranges between 0 and 1.1
×
1
0
-
2 per day [35]. Although the mortality rate of the water bugs is not known, it is assumed to be 0.15 per day [18]. We assume that we have more water bugs than humans so that m
1
>
1. 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 with1000 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.The scatter plots for the parametersα, β
3, μ
d, and μ
W.
(a)
(b)
(c)
(d)Figures3(a) and 3(b) depict parameters with a positive correlation with the reproduction number. They show a monotonic increase of R
T as α and β
3 increase. This means that, to curtail the epidemic, the reduction in the shedding rate and infection of water bugs byM. 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 theM. 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 Figure4. 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 ofM. ulcerans into the environment will significantly slow the epidemic.Figure 4
The tornado plots for the four parameters in the model reproduction number.


## 4.3. Simulation Results

To validate our mathematical analysis results, we plot phase diagrams forR
T 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.The phase diagrams forR
T
=
0.8889 (a) and R
T
=
5.3333 (b).
(a)
(b)Figure6 shows a three-dimensional phase diagram for the human population dynamics. The existence of the endemic equilibrium, when R
T
>
1, 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.Figure 6
A phase diagram for the human population showing the endemic steady state. For a randomly chosen set of initial conditions, all trajectories tend to an endemic equilibrium for the following parameter values:μ
H
=
0.02, θ
=
0.04, Λ
=
0.07, σ
=
0.4, and γ
=
0.7.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μ
d, 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 μ
d. 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 μ
w 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.Fraction of the infected human population forR
T
=
1.6492 for the parameters μ
d and μ
W. The parameter values used for the constant parameters are μ
H
=
0.00045,  m
1
=
10,  θ
=
0.011,  γ
=
0.000016,  β
3
=
0.09,  σ
=
0.08,  K
~
=
0.4000,  N
H
=
100000,  α
=
0.00615,  β
1
=
0.00001, and β
2
=
0.0000002.
(a)
(b)We also explore the role played by directM. 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.Figure 8
The proportion of infected humans for the given values ofβ
2 and the following parameter values: μ
H
=
0.00045;  m
1
=
10;  θ
=
0.011;  γ
=
0.000016;  β
3
=
0.09;  μ
d
=
0.02;  σ
=
0.14;  K
=
0.4000;  μ
W
=
0.025;  N
H
=
100000;  α
=
0.006; and β
1
=
0.00001.The shedding rate ofM. ulcerans into the environment and the treatment rate of humans are important to consider. In Figure 9(a) we observe that increasing the amount ofM. 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 phase diagram for the infected water bugs andM. ulcerans in the environment for the same parameters presented in Figure 7 with μ
d
=
0.02 and μ
W
=
0.025.
(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 numberR
T. 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 andM. 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 ofM. 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.



---

*Source: 576039-2014-08-21.xml*