# A Multiscale Model for the World’s First Parasitic Disease Targeted for Eradication: Guinea Worm Disease

Authors: Rendani Netshikweta; Winston Garira
Journal: Computational and Mathematical Methods in Medicine
 (2017)

---

## Abstract

Guinea worm disease (GWD) is both a neglected tropical disease and an environmentally driven infectious disease. Environmentally driven infectious diseases remain one of the biggest health threats for human welfare in developing countries and the threat is increased by the looming danger of climate change. In this paper we present a multiscale model of GWD that integrates the within-host scale and the between-host scale. The model is used to concurrently examine the interactions between the three organisms that are implicated in natural cases of GWD transmission, the copepod vector, the human host, and the protozoan worm parasite(Dracunculus medinensis), and identify their epidemiological roles. The results of the study (through sensitivity analysis of R0) show that the most efficient elimination strategy for GWD at between-host scale is to give highest priority to copepod vector control by killing the copepods in drinking water (the intermediate host) by applying chemical treatments (e.g., temephos, an organophosphate). This strategy should be complemented by health education to ensure that greater numbers of individuals and communities adopt behavioural practices such as voluntary reporting of GWD cases, prevention of GWD patients from entering drinking water bodies, regular use of water from safe water sources, and, in the absence of such water sources, filtering or boiling water before drinking. Taking into account the fact that there is no drug or vaccine for GWD (interventions which operate at within-host scale), the results of our study show that the development of a drug that kills female worms at within-host scale would have the highest impact at this scale domain with possible population level benefits that include prevention of morbidity and prevention of transmission.


---

## Body

## 1. Introduction

Guinea worm disease (GWD), sometimes known asDracunculiasis or dracontiasis [1], is a nematode infection transmitted to humans exclusively through contaminated drinking water. People become infected when they drink water contaminated with copepods or cyclopoids (tiny aquatic crustaceans) harbouring infective Dracunculus larvae also known asDracunculus medinensis. The larvae ofDracunculus medinensis are released into the stomach, when the copepods are digested by the effect of the gastric juice and get killed by the acid environment. Although the disease has low mortality, its morbidity is considerably high causing huge economic losses and devastating disabilities [2]. There is no vaccine or drug for the disease. Our ability to eliminate GWD rests partly on gaining better insights into the functioning of the immune system, especially its interaction with Guinea worm parasite and partly on development of drugs to treat the disease together with implementation of preventive measures. Currently, the only therapy for GWD is to physically extract the worm from the human body. Humans are the sole definitive host for GWD parasite. Efforts to eradicate the disease are focused on preventive measures which include the following:(a)
Parasite control in the physical water environment. This may involve chlorination of drinking water, or boiling the water before drinking, or applying a larvicide, all of which have the effect of killing the parasite and thereby reduce parasite population in the physical water environment.(b)
Parasite control within the human host. This involves physically extracting the worm from the human body by rolling it over an ordinary stick or matchstick [1, 3] and ensuring that the patient receives care by cleaning and bandaging the wound until all the worms are extracted from the patient. This process may take up to two months to complete, as the worm can grow up to a meter in length and only 1-2 centimeters can be removed per day [4, 5].(c)
Vector control. This consists of killing the copepods in water (the intermediate host) by applying a chemical called temephos, an organophosphate, to unsafe drinking water sources every month during the transmission season, thus reducing vector population and reducing the chances of individuals contracting the disease [2, 6, 7]. The adult vector may also be removed from drinking water by filtering the water using a nylon cloth or by boiling the water.(d)
Health education. This is disseminated through poster, radio and television broadcast, village criers and markets, face-to-face communication (social mobilization and house-to-house visits) by health workers and volunteers to ensure that greater numbers of individuals and communities adopt behavioural practices aimed at preventing transmission of GWD [8]. These behavioural practices include voluntary reporting of GWD cases, prevention of GWD patients from entering drinking water bodies, regular use of water from safe water sources, and, in the absence of such water sources, filtering or boiling water before drinking [6].(e)
Provision of safe water sources. This involves providing safe drinking water supplies through protecting hand-dug wells and sinking deep bore wells, improving existing surface water sources by constructing barriers to prevent humans from entering water, and filtering the water through sand-filters [4].To date, these preventive measures have reduced the incidence of GWD by over 99% [6], making GWD the most likely parasitic disease that will soon be eradicated without the use of any drug or vaccine. Most countries, including the whole of Asia, are now declared free from GWD and transmission of the disease is now limited to African countries, especially Sudan, Ghana, Mali, Niger, and Nigeria [8]. GWD is one of the neglected tropical diseases. It is also an environmentally driven infectious disease. Therefore, its transmission depends on the parasite’s survival in the environment and finding new hosts (humans and copepod vectors) in order to replicate and sustain parasite population. Because this process is complex, it has hampered eradication efforts. During the parasite’s movement through the environment to the human and copepod vector hosts, many environmental factors influence both the parasite’s population and the vector population.For infectious diseases, including environmentally driven infectious diseases such as GWD, mathematical models have a long history of being used to study their transmission and also to compare and evaluate the effectiveness and affordability of intervention strategies that can be used to control or eliminate them [9, 10]. Currently, the predominant focus of modelling of infectious diseases is centered on concepts of epidemiological modelling and immunological modelling being considered as separate disease processes even for the same infectious disease. In epidemiological or between-host modelling of infectious diseases, the focus is on studying of transmission of infectious diseases between hosts, be they animals or humans or even both in the case of multiple host infections. In the immunological or within-host modelling of infectious diseases, the focus is on studying the interaction of pathogen and the immune system together with other within-host processes in order to elucidate outcomes of infection within a single host [11, 12]. To the best of our knowledge there has been no mathematical model to study the multiscale nature of GWD transmission by integrating between-host scale and within-host scale disease processes. Such models are sometimes called immunoepidemiological models [13]. Most of the mathematical models that have been developed so far are focused on the study of GWD at the epidemiological scale [14–16]. The purpose of this study is to develop an immunoepidemiological model of GWD. Immunoepidemiological modelling of infectious diseases is the quantitative approach which assists in developing a systems approach to understanding infectious disease transmission dynamics with regard to the interdependences between epidemiological (between-host scale) and immunological (within-host scale) processes [17, 18]. The immunoepidemiological model of GWD presented in this paper is based on a modelling framework of the immunoepidemiology of environmentally driven infectious diseases developed recently by the authors [13]. This new and innovative immunoepidemiological modelling framework, while maintaining the limits of a mathematical model, offers a solid platform to bring the separate modelling efforts (immunological modelling and epidemiological modelling) that focus on different aspects of the disease processes together to cover a broad range of disease aspects and time-scales in an integrated systems approach. It bridges host, environmental, and parasitic disease phenomena using mathematical modelling of parasite-host-environment-vector interactions and epidemiology to illuminate the fundamental processes of disease transmission in changing environments. For GWD there are three distinct time-scales associated with its transmission cycle which are as follows.(i)
The epidemiological time-scale, which is associated with the infection between hosts (human and copepod vector hosts).(ii)
The within-host time-scale, which is related to the replication and developmental stages of Guinea worm parasite within an individual human host and the individual copepod vector host.(iii)
The environmental time-scale, which is associated with the abundance and survival of Guinea worm parasite population and vector population in the physical water environment.In order to try and integrate these different processes and the associated time-scales of GWD, the immunoepidemiological model of GWD presented here incorporates the actual parasite load of the human host and copepod vector, rather than simply tracking the total number of infected humans. It also incorporates the various stages of the parasite’s life cycle as well as the within-host effects such as the effect of gastric juice within an infected human host and describes how the life stages in the definitive human host, environment, and intermediate vector are interconnected with the parasite’s life cycle through contact, establishment, and parasite fecundity. The paper is organized as follows. In Section2 we present brief discussion of the life cycle of Guinea worm parasite and use this information to develop the immunoepidemiological model of GWD in the same section. In Sections 3, 4, and 5, we derive the analytical results associated with the immunoepidemiological model and show that the model is mathematically and epidemiologically well-posed. We also show the reciprocal influence between the within-host scale and between-host scale of GWD transmission dynamics. The results of the sensitivity analysis of the reproductive number are given in Section 6 while the numerical results of the model are presented in Section 7. The paper ends with conclusions in Section 8.


## 2. The Mathematical Model

We develop a multiscale model of Guinea worm disease that traces the parasite’s life cycle of Guinea worm disease. The life cycle of GWD involves three different environments: physical water environment, biological human host environment, and biological copepod host environment. For more details on the life cycle of GWD see the published works [6, 19]. We only give a brief description in this section. The transmission cycle of Guinea worm disease begins when the human individual drinks contaminated water with copepods that are infected with Guinea mature worm larvae (L3 larvae). After ingestion, gastric juice in the human stomach kills the infected copepods and mature worm larvae are released. Then the released mature worm larvae penetrate the human stomach and intestinal wall and move to abdominal tissues where they grow and mate. After mating the male worms die soon and fertilized female worms migrate towards the skin surface (usually on the lower limbs or feet). After a year of infection, the fertilized female worm makes a blister on the infected individual’s skin causing burning and itching, which forces an infected individual to immerse his or her feet into water (which is the only source of drinking water) to seek relief from pain. At that point the female worm emerges and releases thousands of worm eggs. The worm eggs then hatch Guinea worm larvae (L1 larvae stage) which are then consumed by copepods and take approximately two weeks to develop and become infective mature larvae (L3 larvae) within the copepods. Then ingestion of the infected copepods by human closes the life cycle. The multiscale model which we now present explicitly traces this life cycle ofDracunculus medinensis in three different environments, which are physical water environment, biological human environment, and biological copepod environment. The model flow diagram is shown in Figure 1.Figure 1
A conceptual diagram of the multiscale model of Guinea worm disease transmission dynamics.The full multiscale model presented in this paper is based on monitoring the dynamics of ten populations at any timet, which are susceptible humans SH(t) and infected humans IH(t) in the behavioural human environment; infected copepods IC in the human biological environment; mature Guinea worms WM(t) and fertilized female Guinea worms WF(t) in the biological human environment (within-host parasite dynamics); Guinea worm eggs EW(t) and Guinea worm larvae LW(t) in the physical water environment; susceptible copepods SE(t) and infected copepods IE(t) in the physical water environment; and gastric juice GJ(t) in the human biological environment. We make the following assumptions for the model:(a)
There is no vertical transmission of the disease.(b)
The transmission of the disease in the human population is only through drinking contaminated water with infected copepods,IE(t), harbouring infective free-living pathogens (first-stage larvae), LW(t), in the physical water environment.(c)
For an infected individual, more than one Guinea worm can emerge simultaneously or sequentially over the course of weeks, depending on the number and intensity of infection the preceding year.(d)
Humans do not develop temporary or permanent immunity.(e)
Copepods do not recover from infection.(f)
The total population of humans and copepods is constant.(g)
Except for the effects of gastric juice in the stomach, there is no immune response in the human host.(h)
Copepods die in the human stomach due to the effects of gastric juice at a rateαC before their larvae undergoes two molts in the copepod to become L3 larvae and therefore are nonviable and noninfectious larvae.From the model flow diagram presented in Figure1 and the assumptions that we have now made, we have the following system of ordinary differential equations as our multiscale model for GWD transmission dynamics:(1)1dSHtdt=ΛH-λHtSHt-μHSHt+αHIHt,2dIHtdt=λHtSHt-μH+δH+αHIHt,3dICtdt=λhtSht-μCGJtICt-αCICt,4dWMtdt=NCμCGJtICt-αM+μMWMt,5dWFtdt=αM2WMt-μF+αFWFt,6dGJtdt=G0+αJGJtICt-μJGJt,7dEWtdt=αFWFtIht-μW+αWEWt,8dLWtdt=NWαWEWt-μLLWt,9dSEtdt=ΛE-λEtSEt-μESEt,10dIEtdt=λEtSEt-μE+δEIEt,where(2)λHt=βHIEtP0+ϵIEt,λEt=βELWtL0+ϵLWt,Iht=IHt+1,Sht=SHt-1,λht=βHIEtP0+ϵIEtIHt+1.Equations (1) and (2) of the model system (1) describe the evolution with time of susceptible and infected human hosts, respectively. At any time t, new susceptible humans are recruited at a constant rate ΛH and we assume that the recruited humans are all susceptible. Susceptible individuals leave the susceptible class either through infection at rate λH(t)SH(t) by drinking contaminated water with infected copepods to join infected group or through natural death at a rate μH. The infected group is generated through infection when susceptible humans acquire the disease at a rate λH(t)SH(t) through drinking water contaminated with copepods infected withDracunculus medinensis. Infected humans leave the infected group either through recovery at a rate αH to join the susceptible group or through natural death at a rate μH, or through disease induced death at a rate δH. Equation (3) of the model system (1) represents the evolution with time of infected copepods within an infected human host. The infected copepods within a human host are generated following uptake of infected copepods in the physical water environment through drinking contaminated water. In the human population, this uptake of infected copepods, which harbour Guinea worm larvae, is the transmission of Guinea worm parasite from the physical water environment to susceptible humans who become infected humans. Following the methodology described in [13] for modelling reinfection (superinfection) for environmentally transmitted infectious disease systems (because GWD and schistosomiasis are both water-borne and vector-borne infections), we model the average rate at which a single susceptible human host uptakes the infected copepods in the physical water environment through drinking contaminated water and becomes an infected human host by the expression(3)λhtSht=λHtSHt-1IHt+1,where λH(t), SH(t), and IH(t) are as defined previously. This is because, in our case, we define such a single infection by a single transition(4)SHt,IHt,IEt⟶SHt-1,IHt+1,IEt.Therefore, the average number of infected copepods,IC(t), within a single infected human host increases at a mean rate λh(t)Sh(t) and decreases through death due to digestion by human gastric acid at a rate μC after their larvae undergo two molts in the copepod to become L3 larvae and release viable and infectious larvae or naturally at a rate αC before their larvae undergo two molts in the copepod to become L3 larvae and release nonviable and noninfectious larvae.Equations (4–6) of the model system (1) represent changes with time of the average population of mature worms WM(t), fertilized female worms WF(t), and the amount of gastric acid GJ(t) within a single infected human host, respectively. The average mature worm population WM(t) in a single infected human host is generated following the digestion of infected copepods in the human stomach by gastric acid and then mature worms are released. We assume that mature worms die naturally at a rate μM and they exit the human stomach to the abdominal tissues at a rate αM, where they grow and mate. The population of fertilized female worms, WF(t) within an infected human host, is generated following the developmental changes undergone by mature fertilized female worms. These developmental changes result in mature worms reaching sexual maturity and mating and all male worms die soon after mating. We assume that fertilized female worms die naturally at a rate μF and emerge out through an infected human individual’s skin (usually the lower limbs) to release Guinea worm eggs into a water source at a rate αF, when an infected human comes into contact with water. The average amount of gastric acid inside a human stomach is generated following copepod vector induced proliferation at a rate αJIC(t), which is proportional to the density of infected copepods within an infected human host. We assume that the amount of gastric acid is also increased by the spontaneous production of gastric acid by the human body at a rate G0 and diluted or degraded at a rate μJ. Equation (7) of model system (1) describes the evolution with time of the Guinea worm eggs EW(t) in the physical water environment. We note that the population of Guinea worm eggs increases when each infected human host excretes eggs at a rate αFWF(t). Therefore the rate at which infected humans contaminate the physical water environment by excreting Guinea worm eggs is modelled by αFWF(t)Ih(t). The last three equations of the model system (1) describe the evolution with time of Guinea worm larvae LW(t), susceptible copepods SE(t), and infected copepods IE(t) in the physical water environment, respectively. The population of Guinea worm larvae is generated through each egg hatching an average of NW worms larvae with eggs hatching at an average rate of αW. Therefore the total Guinea worm larvae in the physical water environment are modelled by NWαWEW(t). We assume that worm larvae in the physical water environment die naturally at a constant rate μL. Similar to human population, at any time t, new susceptible copepods are recruited at a constant ΛE. Susceptible copepods leave the susceptible group to join the infected copepods group through infection at a rate λE(t)SE(t) when they consume first-stage Guinea worm larvae in the physical water environment. We assume that the population of copepods die naturally at a constant rate μE and further, we also assume that infected copepods have an additional mortality rate δE due to infection. The model state variables are summarized in Table 1.Table 1
Description of the state variables of the model system (1).
State variable
Description
Initial value
S
H
(
t
)
The susceptible human population size in the behavioural human environment
2500
I
H
(
t
)
The infected human population size in the behavioural human environment
10
I
C
(
t
)
The infected copepod population size in the biological human environment
0
W
M
(
t
)
The mature worm population size in the biological human environment
0
W
F
(
t
)
The female worm population size in the biological human environment
0
G
J
(
t
)
Amount of gastric acid in the human stomach
1.5
S
E
(
t
)
The susceptible copepod population size in the physical water environment
10
5
I
E
(
t
)
The infected copepod population size in the physical water environment
0
E
W
(
t
)
The worm egg population size in the physical water environment
0
L
W
(
t
)
The worm larvae population size in the physical water environment
5000


## 3. Invariant Region of the Model

The model system (1) can be analysed in a region Ω⊂R+10 of biological interest. Now assume that all parameters and state variables for model system (1) are positive for all t>0 and further suppose that GJ is bounded above by G0/μJ. It can be shown that all solutions for the model system (1) with positive initial conditions remain bounded.LettingNH=SH+IH and adding (1) and (2) of model system (1) we obtain(5)dSHdt+dIHdt=dNHdt=ΛH-μHNH-δHIH≤ΛH-μHNH.This implies that(6)limt→∞sup⁡NHt≤ΛHμH.Similarly, lettingNE=SE+IE and adding (9) and (10) of model system (1) we obtain(7)dSEdt+dIEdt=dNEdt=ΛE-μHNE-δHIE≤ΛE-μHNE.This also implies that(8)limt→∞sup⁡NEt≤ΛEμE.Now considering the third equation of model system (1), given by(9)dICdt=λhSh-μCGJIC-δCIC=βHIESH-1P0+ϵIEIH+1-μCGJ+αCIC,we obtain(10)dICdt≤βHΛEΛH-μHP0μE+ϵΛEΛH+μH-1μJμCG0+αCμJIC.This implies that(11)limt→∞sup⁡ICt≤βHΛEΛH-μHP0μE+ϵΛEΛH+μHμJμCG0+δHμJ.Using (6), (8), and (11) similar expression can be derived for the remaining model variables. Hence, all feasible solutions of the model system (1) are positive and enter a region defined by(12)Ω=SH,IH,IC,WM,WF,GJ,EW,LW,SE,IE∈R+10:0≤SH+IH≤S1,0≤SE+IE≤S2,0≤IC≤S3,0≤WM≤S4,0≤WF≤S5,0≤GJ≤S6,0≤EW≤S7,0≤LW≤S8,which is positively invariant and attracting for allt>0, where(13)S1=ΛHμH,S2=ΛEμE,S3=μJμCG0+αCμJS9,S4=NCαCμCαM+μMμJμCG0+αCμJS9,S5=12αMαF+μFNCαCμCαM+μMμJμCG0+αCμJS9,S6=G0μJ,S7=12αFαW+μWαMαF+μFNCαCμCαM+μMμJμCG0+αCμJS9,S8=12NWαWμLαFαW+μWαMαF+μFNCαCμCαM+μMμJμCG0+αCμJS9,S9=βHΛCΛH-μHP0μC+ϵΛCΛH+μH.Therefore it is sufficient to consider solutions of the model system (1) in Ω, since all solutions starting in Ω remain there for all t≥0. Hence, the model system is mathematically and epidemiologically well-posed and it is sufficient to consider the dynamics of the flow generated by model system (1) in Ω whenever ΛH>μH. We shall assume in all that follows (unless stated otherwise) that ΛH>μH.


## 4. Determination of Disease-Free Equilibrium and Its Stability

To obtain the disease-free equilibrium point of system (1), we set the left-hand side of the equations equal to zero and further we assume that IH=IC=WH=WH=EW=LW=IE=0. This means that all the populations are free from the disease. Thus we get(14)E0=SH0,IH0,IC0,WM0,WF0,GJ0,EW0,LW0,SE0,IE0,=ΛHμH,0,0,0,0,G0μJ,0,0,ΛEμE,0,as the disease-free equilibrium of the model system (1).


### 4.1. The Basic Reproduction Number of the Model System (1)

The basic reproduction number of the system model (1) is calculated in this section using next generation operator approach described in [20]. Thus the model system (1) can also be written in the form(15)dXdt=fX,Y,Z,dYdt=gX,Y,Z,dZdt=hX,Y,Z,where(i)
X=(SH,SE,GJ) represents all compartments of individuals who are not infected,(ii)
Y=(IH,IC,WM,WF,EW) represents all compartments of infected individuals who are not capable of infecting others,(iii)
Z=(IE,LW) represents all compartments of infected individuals who are capable of infecting others.We also let the disease-free equilibrium of the model (1) be denoted by the following expression:(16)U¯0=ΛHμH,0,0,0,0,G0μJ,0,0,ΛEμE,0.Following [20], we let(17)g~X∗,Z=g~1X∗,Z,g~2X∗,Z,g~3X∗,Z,g~4X∗,Z,g~5X∗,Z,with(18)g~1X∗,Z=βHΛHZ1μHμH+δH+αHP0+ϵIE,g~2X∗,Z=βHΛH-μHμJμH+δH+αHZ1μCG0+μJαCM11,g~3X∗,Z=NCαCβHΛH-μHG0μH+δH+αHZ1μM+αMμCG0+μJαCM11,g~4X∗,Z=αMNCαCβHΛH-μHG0μH+δH+αHZ12μF+αFμM+αMμCG0+μJαCM11,g~5X∗,Z=αFαMNCαCβHΛH-μHG0Z12μW+αWμF+αFμM+αMμCG0+μJαCμHP0+ϵIE,where(19)M11=ϵβHΛHIE+μHμH+δH+μHP0+ϵIE.We deduce that(20)hX,Y,Z=h1X,Y,Z,h2X,Y,Z,with(21)h1X,Y,Z=λESE-μE+αEIE=βEΛEZ2μEL0+ϵZ2-μE+αEZ1,h2X,Y,Z=NWαWEW-μLLW=KZ1P0+ϵZ1-μLZ2,where(22)K=αFαMNWαWNCαCβHΛH-μHG02μW+αWμF+αFμM+αMμCG0+μJαCμH.A matrix(23)A=DZhX∗,g~X∗,0,0=-μE+αEKP0βEΛEμEL0-μLcan be written in the formA=M-D, so that(24)M=0KP0βEΛEμEL00,(25)D=μE+αE00μL.The basic reproductive number is the spectral radius (dominant eigenvalue) of the matrixT=MD-1. Hence, the basic reproduction number of the immumoepidemiological model (1) is expressed by the following quantity.(26)R0=12·αMαM+μM·αFαF+μF·NCμCG0μCG0+μJαC·βHΛH-μHP0μH·NWαWαW+μWμLβEΛEμEμE+δEL0=R0BR0W,with(27)R0B=βHΛH-μHP0μH·NWαWαW+μWμLβEΛEμEμE+δEL0,(28)R0W=12·αMαM+μM·αFαF+μF·NCμCG0μCG0+μJαC.The expression,R0B, in (27) represents GWD’s partial reproductive number associated with the between-host transmission of the disease while the expression, R0W, in (28) represents GWD’s partial reproductive number associated with the within-host transmission of the disease. From the above two expressions in (27) and (28), respectively, we therefore make the following deductions.(i)
The epidemiological (between-host) transmission parameters such as the rate at which susceptible humans come into contact with water contaminated with infected copepodsβH (through drinking contaminated water with infected copepods) and the rate at which susceptible copepods come into contact with Guinea worm larvae βE; the supply rate of susceptible humans ΛH and copepods ΛE (through birth); the rate at which worms emerge from infected humans to contaminate the physical water environment αF, by laying eggs every time infected humans come into contact with water sources; the rate at which eggs in physical water environment hatch to produce worm larvae NWαW all contribute to the transmission of Guinea worm disease. Therefore control measures such as reducing the rate at which infected human hosts visit water sources when an individual is infected, reducing contact rate between susceptible humans with contaminated water through educating the public, and treating water bodies with chemicals that kill worm eggs, worm larvae, and copepods may help to reduce the transmission risk of GWD.(ii)
The immunological (within-host) transmission parameters such as the rate at which infected copepods within an infected human host release mature wormsNCμC after digestion by human gastric juice; the rate at which mature worms become fertilized females worms αM/2; and the rate at which mature worms and females worms die all contribute to the transmission of Guinea worm disease. Therefore immune mechanisms that kill infected copepods and worms within infected human host and also treatment intend to kill both mature worms and fertilized female worm population may help to reduce the transmission risk of GWD.Therefore, both the epidemiological and immunological factors affect the transmission cycle of GWD in both humans and copepod population.


### 4.2. Local Stability of DFE

In this section we determine the local stability of DFE of the model system (1). We linearize equations of the model system (1) in order to obtain a Jacobian matrix. Then we evaluate the Jacobian matrix of the system at the disease-free equilibrium (DFE),(29)E0=ΛHμH,0,0,0,0,G0μJ,0,0,ΛEμE,0.The Jacobian matrix of the model system (1) evaluated at the disease-free equilibrium state (DFE) is given by(30)JE0=-μHαH0000000-A00-q00000000A000-q1000000A100NCμCG0μJ-q2000000000αM2-q30000000αJG0μJ00-μJ00000000αF0-q4000000000NWαW-μL000000000-βEΛEμEL0-μE00000000βEΛEμEL00-q5,where(31)q0=μH+δH+αH,q1=μCG0+μJαCμJ,q2=μM+αM,q3=μF+αF,q4=μW+αW,q5=μE+αE,A0=βHΛHP0μH,A1=βHΛH-μHμHP0.We consider stability of DFE by calculating the eigenvalues (λs) of the Jacobian matrix given by (30). The characteristic equation for the eigenvalues is given by(32)λ0λ6+π1λ5+π2λ4+π3λ3+π4λ2+π5λ+π6=0,where(33)λ0=-μH-λ-μE-λ-μJ-λ-q0-λ.It is clear from (32) that there are four negative eigenvalues (-μH, -μE, -μJ, and -q0). Now in order to make conclusions about the stability of the DFE, we use the Routh-Hurwitz criteria to determine the sign of the remaining eigenvalues of the polynomial(34)λ6+π1λ5+π2λ4+π3λ3+π4λ2+π5λ+π6=0,where(35)π1=q1+q2+q3+q4+q5+μL,π2=q1q2+q3q4+q1+q2+q3+q4q5+μL+q5μL+q1+q2q3+q4,π3=q1q2q3+q4+q3q4q1+q2+q1+q2q3+q4q5+μL+q1q2q5+μL+q3q4q5+μL+q5μLq1+q2+q3+q4,π4=q1q2q3q4+q3q4q1+q2q5+μL+q1q2q3+q4q5+μL+q1+q2q3+q4q1μL+q5q4q3μL+q1q2q5μL,π5=q1q2q3q4q5+μL+q3q4q1+q2q5μL+q1q2q3+q4q5μL,π6=q1q2q3q4q5μL1-R02.Using the Routh-Hurwitz stability criterion, the equilibrium state associated with the model system (1) is stable if and only if the determinants of all the Hurwitz matrices associated with the characteristic equation (34) are positive; that is,(36)Det⁡Hj>0;j=1,2,…,6,where(37)H1=π1;H2=π11π3π2;H3=π110π3π2π1π5π4π3;H4=π1100π3π2π11π5π4π3π20π6π5π4;H5=π11000π3π2π110π5π4π3π2π10π6π5π4π3000π6π5;H6=π110000π3π2π1100π5π4π3π2π110π6π5π4π3π2000π6π5π400000π6.The Routh-Hurwitz criterion applied to (37) requires that the following conditions (H1)–(H6) be satisfied, in order to guarantee the local stability of the disease-free equilibrium point of the model system (1).(H1)
π1>0.(H2)
π1π2-π3>0.(H3)
π1(π2π3+π5)>π1π4+π32.(H4)
π1[π2(π3(π4+π5)+π1π6)+(π1+π4)]>π1[π22π5+π3π6+π1π42]+π32π4+π52.(H5)
π6[π1(2π2π5+π3(π1π4-3π5-π3))+π33π6]+π5[π5(2π1π4+π2π3-π1π2π2+1+π4π1π2π3-π12π4-π32)]>0.(H6)
π62[π1(2π2π5+π3(π1π4-3π5-π3))+π33π6]+π5π6[π5(2π1π4+π2π3-π1π2π2+1+π4π1π2π3-π12π4-π32)]>0.From (37) we note that all the coefficients π1, π2, π3, π4, π5, and π6 of the polynomial P(λ) are greater than zero whenever R02<1. And we also noted that the conditions above are satisfied if and only if R02<1. Hence all the roots of the polynomial P(λ) either are negative or have negative real parts. The results are summarized in the following theorem.Theorem 1.
The disease-free equilibrium point of the model system (1) is locally asymptotically stable whenever R0<1.


### 4.3. Global Stability of DFE

To determine the global stability of DFE of the model system (1), we use Theorem 2 in [21] to establish that the disease-free equilibrium is globally asymptotically stable whenever R0<1 and unstable when R0>1. In this section, we list two conditions that if met, also guarantee the global asymptotic stability of the disease-free state. We write the model system (1) in the form(38)dXdt=FX,Z,dYdt=GX,Z,where(i)
X=(SH,SE,GJ) represents all uninfected components.(ii)
Z=(IH,IC,WM,WF,EW,LW,IE) represents all compartments of infected and infectious components.We let(39)U0=X∗,0=ΛHμH,0,0,ΛCμC,0,0,0denote the disease-free equilibrium (DFE) of the system. To guarantee global asymptotic stability of the disease-free equilibrium, conditions (H1) and (H2) below must be met [20].(H1)
dX/dt=F(X,0) is globally asymptotically stable,(H2)
G(X,Z)=AZ-G^(X,Z) and G^(X,Z)≥0 for X,Z∈R+10, where A=DZGX∗,0 is an M-matrix and R+10 is the region where the model makes biological sense.In our case(40)FX,0=ΛH-μHSHΛE-μESEG0-μJGJ.MatrixA is given by(41)A=-a000000βHΛHP0μH0-a10000βHΛH-μHμHP00NCμCG0μJ-a2000000αM2-a3000000αF-a4000000NWαW-μL000000βEΛEL0μE-a5,where(42)a0=μH+δH+αH,a1=1μJμC+αCμJ,a2=μM+αM,a3=μF+αF,a4=μW+αW,a5=μE+αE,G^X,Z=ΛHμHP0-SHP0+ϵIEβHIEΛH-μHμHP0-SH-1P0+ϵIEβHIE+ICμCGJ-G0μJ+αC1-μJ0000ΛEμEL0-SEL0+ϵLWβELW.Assume thatGJ=G0/μJ and μJ∈[0,1]. It is clear that G^(X,Z)≥0 for all (X,Z)∈R+10, since ΛH/μHP0≥SH/P0+ϵIE, ΛE/μEL0≥SE/L0+ϵLW, and (ΛH-μH)/μHP0≥(SH-1)/P0+ϵIC provided that ΛH>μH. It is also clear that A is an M-matrix, since the off diagonal elements of A are nonnegative. We state a theorem which summarizes the above result.Theorem 2.
The disease-free equilibrium of model system (1) is globally asymptotically stable if R0≤1 and the assumptions (H1) and (H2) are satisfied.


## 4.1. The Basic Reproduction Number of the Model System (1)

The basic reproduction number of the system model (1) is calculated in this section using next generation operator approach described in [20]. Thus the model system (1) can also be written in the form(15)dXdt=fX,Y,Z,dYdt=gX,Y,Z,dZdt=hX,Y,Z,where(i)
X=(SH,SE,GJ) represents all compartments of individuals who are not infected,(ii)
Y=(IH,IC,WM,WF,EW) represents all compartments of infected individuals who are not capable of infecting others,(iii)
Z=(IE,LW) represents all compartments of infected individuals who are capable of infecting others.We also let the disease-free equilibrium of the model (1) be denoted by the following expression:(16)U¯0=ΛHμH,0,0,0,0,G0μJ,0,0,ΛEμE,0.Following [20], we let(17)g~X∗,Z=g~1X∗,Z,g~2X∗,Z,g~3X∗,Z,g~4X∗,Z,g~5X∗,Z,with(18)g~1X∗,Z=βHΛHZ1μHμH+δH+αHP0+ϵIE,g~2X∗,Z=βHΛH-μHμJμH+δH+αHZ1μCG0+μJαCM11,g~3X∗,Z=NCαCβHΛH-μHG0μH+δH+αHZ1μM+αMμCG0+μJαCM11,g~4X∗,Z=αMNCαCβHΛH-μHG0μH+δH+αHZ12μF+αFμM+αMμCG0+μJαCM11,g~5X∗,Z=αFαMNCαCβHΛH-μHG0Z12μW+αWμF+αFμM+αMμCG0+μJαCμHP0+ϵIE,where(19)M11=ϵβHΛHIE+μHμH+δH+μHP0+ϵIE.We deduce that(20)hX,Y,Z=h1X,Y,Z,h2X,Y,Z,with(21)h1X,Y,Z=λESE-μE+αEIE=βEΛEZ2μEL0+ϵZ2-μE+αEZ1,h2X,Y,Z=NWαWEW-μLLW=KZ1P0+ϵZ1-μLZ2,where(22)K=αFαMNWαWNCαCβHΛH-μHG02μW+αWμF+αFμM+αMμCG0+μJαCμH.A matrix(23)A=DZhX∗,g~X∗,0,0=-μE+αEKP0βEΛEμEL0-μLcan be written in the formA=M-D, so that(24)M=0KP0βEΛEμEL00,(25)D=μE+αE00μL.The basic reproductive number is the spectral radius (dominant eigenvalue) of the matrixT=MD-1. Hence, the basic reproduction number of the immumoepidemiological model (1) is expressed by the following quantity.(26)R0=12·αMαM+μM·αFαF+μF·NCμCG0μCG0+μJαC·βHΛH-μHP0μH·NWαWαW+μWμLβEΛEμEμE+δEL0=R0BR0W,with(27)R0B=βHΛH-μHP0μH·NWαWαW+μWμLβEΛEμEμE+δEL0,(28)R0W=12·αMαM+μM·αFαF+μF·NCμCG0μCG0+μJαC.The expression,R0B, in (27) represents GWD’s partial reproductive number associated with the between-host transmission of the disease while the expression, R0W, in (28) represents GWD’s partial reproductive number associated with the within-host transmission of the disease. From the above two expressions in (27) and (28), respectively, we therefore make the following deductions.(i)
The epidemiological (between-host) transmission parameters such as the rate at which susceptible humans come into contact with water contaminated with infected copepodsβH (through drinking contaminated water with infected copepods) and the rate at which susceptible copepods come into contact with Guinea worm larvae βE; the supply rate of susceptible humans ΛH and copepods ΛE (through birth); the rate at which worms emerge from infected humans to contaminate the physical water environment αF, by laying eggs every time infected humans come into contact with water sources; the rate at which eggs in physical water environment hatch to produce worm larvae NWαW all contribute to the transmission of Guinea worm disease. Therefore control measures such as reducing the rate at which infected human hosts visit water sources when an individual is infected, reducing contact rate between susceptible humans with contaminated water through educating the public, and treating water bodies with chemicals that kill worm eggs, worm larvae, and copepods may help to reduce the transmission risk of GWD.(ii)
The immunological (within-host) transmission parameters such as the rate at which infected copepods within an infected human host release mature wormsNCμC after digestion by human gastric juice; the rate at which mature worms become fertilized females worms αM/2; and the rate at which mature worms and females worms die all contribute to the transmission of Guinea worm disease. Therefore immune mechanisms that kill infected copepods and worms within infected human host and also treatment intend to kill both mature worms and fertilized female worm population may help to reduce the transmission risk of GWD.Therefore, both the epidemiological and immunological factors affect the transmission cycle of GWD in both humans and copepod population.


## 4.2. Local Stability of DFE

In this section we determine the local stability of DFE of the model system (1). We linearize equations of the model system (1) in order to obtain a Jacobian matrix. Then we evaluate the Jacobian matrix of the system at the disease-free equilibrium (DFE),(29)E0=ΛHμH,0,0,0,0,G0μJ,0,0,ΛEμE,0.The Jacobian matrix of the model system (1) evaluated at the disease-free equilibrium state (DFE) is given by(30)JE0=-μHαH0000000-A00-q00000000A000-q1000000A100NCμCG0μJ-q2000000000αM2-q30000000αJG0μJ00-μJ00000000αF0-q4000000000NWαW-μL000000000-βEΛEμEL0-μE00000000βEΛEμEL00-q5,where(31)q0=μH+δH+αH,q1=μCG0+μJαCμJ,q2=μM+αM,q3=μF+αF,q4=μW+αW,q5=μE+αE,A0=βHΛHP0μH,A1=βHΛH-μHμHP0.We consider stability of DFE by calculating the eigenvalues (λs) of the Jacobian matrix given by (30). The characteristic equation for the eigenvalues is given by(32)λ0λ6+π1λ5+π2λ4+π3λ3+π4λ2+π5λ+π6=0,where(33)λ0=-μH-λ-μE-λ-μJ-λ-q0-λ.It is clear from (32) that there are four negative eigenvalues (-μH, -μE, -μJ, and -q0). Now in order to make conclusions about the stability of the DFE, we use the Routh-Hurwitz criteria to determine the sign of the remaining eigenvalues of the polynomial(34)λ6+π1λ5+π2λ4+π3λ3+π4λ2+π5λ+π6=0,where(35)π1=q1+q2+q3+q4+q5+μL,π2=q1q2+q3q4+q1+q2+q3+q4q5+μL+q5μL+q1+q2q3+q4,π3=q1q2q3+q4+q3q4q1+q2+q1+q2q3+q4q5+μL+q1q2q5+μL+q3q4q5+μL+q5μLq1+q2+q3+q4,π4=q1q2q3q4+q3q4q1+q2q5+μL+q1q2q3+q4q5+μL+q1+q2q3+q4q1μL+q5q4q3μL+q1q2q5μL,π5=q1q2q3q4q5+μL+q3q4q1+q2q5μL+q1q2q3+q4q5μL,π6=q1q2q3q4q5μL1-R02.Using the Routh-Hurwitz stability criterion, the equilibrium state associated with the model system (1) is stable if and only if the determinants of all the Hurwitz matrices associated with the characteristic equation (34) are positive; that is,(36)Det⁡Hj>0;j=1,2,…,6,where(37)H1=π1;H2=π11π3π2;H3=π110π3π2π1π5π4π3;H4=π1100π3π2π11π5π4π3π20π6π5π4;H5=π11000π3π2π110π5π4π3π2π10π6π5π4π3000π6π5;H6=π110000π3π2π1100π5π4π3π2π110π6π5π4π3π2000π6π5π400000π6.The Routh-Hurwitz criterion applied to (37) requires that the following conditions (H1)–(H6) be satisfied, in order to guarantee the local stability of the disease-free equilibrium point of the model system (1).(H1)
π1>0.(H2)
π1π2-π3>0.(H3)
π1(π2π3+π5)>π1π4+π32.(H4)
π1[π2(π3(π4+π5)+π1π6)+(π1+π4)]>π1[π22π5+π3π6+π1π42]+π32π4+π52.(H5)
π6[π1(2π2π5+π3(π1π4-3π5-π3))+π33π6]+π5[π5(2π1π4+π2π3-π1π2π2+1+π4π1π2π3-π12π4-π32)]>0.(H6)
π62[π1(2π2π5+π3(π1π4-3π5-π3))+π33π6]+π5π6[π5(2π1π4+π2π3-π1π2π2+1+π4π1π2π3-π12π4-π32)]>0.From (37) we note that all the coefficients π1, π2, π3, π4, π5, and π6 of the polynomial P(λ) are greater than zero whenever R02<1. And we also noted that the conditions above are satisfied if and only if R02<1. Hence all the roots of the polynomial P(λ) either are negative or have negative real parts. The results are summarized in the following theorem.Theorem 1.
The disease-free equilibrium point of the model system (1) is locally asymptotically stable whenever R0<1.


## 4.3. Global Stability of DFE

To determine the global stability of DFE of the model system (1), we use Theorem 2 in [21] to establish that the disease-free equilibrium is globally asymptotically stable whenever R0<1 and unstable when R0>1. In this section, we list two conditions that if met, also guarantee the global asymptotic stability of the disease-free state. We write the model system (1) in the form(38)dXdt=FX,Z,dYdt=GX,Z,where(i)
X=(SH,SE,GJ) represents all uninfected components.(ii)
Z=(IH,IC,WM,WF,EW,LW,IE) represents all compartments of infected and infectious components.We let(39)U0=X∗,0=ΛHμH,0,0,ΛCμC,0,0,0denote the disease-free equilibrium (DFE) of the system. To guarantee global asymptotic stability of the disease-free equilibrium, conditions (H1) and (H2) below must be met [20].(H1)
dX/dt=F(X,0) is globally asymptotically stable,(H2)
G(X,Z)=AZ-G^(X,Z) and G^(X,Z)≥0 for X,Z∈R+10, where A=DZGX∗,0 is an M-matrix and R+10 is the region where the model makes biological sense.In our case(40)FX,0=ΛH-μHSHΛE-μESEG0-μJGJ.MatrixA is given by(41)A=-a000000βHΛHP0μH0-a10000βHΛH-μHμHP00NCμCG0μJ-a2000000αM2-a3000000αF-a4000000NWαW-μL000000βEΛEL0μE-a5,where(42)a0=μH+δH+αH,a1=1μJμC+αCμJ,a2=μM+αM,a3=μF+αF,a4=μW+αW,a5=μE+αE,G^X,Z=ΛHμHP0-SHP0+ϵIEβHIEΛH-μHμHP0-SH-1P0+ϵIEβHIE+ICμCGJ-G0μJ+αC1-μJ0000ΛEμEL0-SEL0+ϵLWβELW.Assume thatGJ=G0/μJ and μJ∈[0,1]. It is clear that G^(X,Z)≥0 for all (X,Z)∈R+10, since ΛH/μHP0≥SH/P0+ϵIE, ΛE/μEL0≥SE/L0+ϵLW, and (ΛH-μH)/μHP0≥(SH-1)/P0+ϵIC provided that ΛH>μH. It is also clear that A is an M-matrix, since the off diagonal elements of A are nonnegative. We state a theorem which summarizes the above result.Theorem 2.
The disease-free equilibrium of model system (1) is globally asymptotically stable if R0≤1 and the assumptions (H1) and (H2) are satisfied.


## 5. The Endemic Equilibrium State and Its Stability

At the endemic equilibrium humans are infected by copepods that have been infected by first-stage larvae(LW). The endemic equilibrium point of the model system (1) given by(43)E^1=SH∗,IH∗,IC∗,WM∗,WF∗,GJ∗,EW∗,LW∗,SE∗,IE∗satisfies(44)0=ΛH-λH∗SH∗-μHSH∗+αHIH∗,0=λH∗SH∗-μH+δH+αHIH∗,0=λh∗Sh∗-μCGJ∗IC∗-αCIC∗,0=NCμCGJ∗IC∗-αM+μMWM∗,0=αM2WM∗-μF+αFWF∗,0=G0+αJGJ∗IC∗-μJGJ∗,0=αFWF∗Ih∗-μW+αWEW∗,0=NWαWEW∗-μLLW∗,0=ΛE-λE∗SE∗-μESE∗,0=λE∗SE∗-μE+δEIE∗,for allSH∗,IH∗,IC∗,WM∗,WF∗,GJ∗,EW∗,LW∗,SE∗,IE∗>0. We therefore obtain the following endemic values. The endemic value of susceptible humans is given by (45)SH∗=ΛH+αHIH∗λH∗+μH.From (45) we note that the susceptible human population at endemic equilibrium is proportional to the average time of stay in the susceptible class and the rate at which new susceptible individuals are entering the susceptible class either through birth or through infected individuals who recover from the disease. Individuals leave the susceptible class through either infection or death. The endemic value of infected humans is given by(46)IH∗=λH∗SH∗μH+δH+αH.We note from (46) that the population of infected humans at the endemic equilibrium point is proportional to the average time of stay in the infected class, the rate at which susceptible individuals become infected, and the density of susceptible individuals. The endemic value of infected copepods population within a single infected human at the equilibrium point is given by(47)IC∗=λH∗SH∗-1IH∗+1μCGJ∗+αC,where SH∗>1. From (47) we note that the average infected copepod population within a single infected human is proportional to the average life-span of infected copepods within a single infected human host and the rate of infection of a single susceptible individual to become infected. We also note that this expression provides a link between the dynamics of the infected copepods within-host and human population dynamics. The endemic value of mature worm population within a single infected human is given by(48)WM∗=NCμCGJ∗IC∗αM+μM.We note from (48) that the population of mature worms within a single infected human at endemic equilibrium point is proportional to the average life-span of mature worms and the rate at which mature worms are released after infected copepods within human host have been killed by human gastric juice. The endemic value of fertilized female worm population within a single infected human is given by(49)WF∗=12αMWM∗αF+μF.The average population of fertilized female worms within an infected human at endemic equilibrium point is equal to the average life-span of female worms and the rate at which mature worms become fertilized female worms. The endemic value of a single human gastric juice is given by(50)GJ∗=G0μJ-αJIC∗,where μJ>αJIC∗. The endemic value of Guinea worm eggs population in the physical water environment is given by(51)EW∗=αFWF∗IH∗+1αW+μW.We note from (51) that the worm egg population at equilibrium point is proportional to the average life-span of eggs, the rate at which each infected human host excretes Guinea worm eggs, and the total number of infected humans. The endemic value of Guinea worm larva population in the physical water environment is given by(52)LW∗=NWαWEW∗μL.We note from (52) that the larvae population at equilibrium point is proportional to the rate at which Guinea worm eggs hatch, the number of larvae generated by each egg, and the average life-span of larvae. The value of susceptible copepod population at equilibrium point is given by(53)SE∗=ΛEλE∗+μE.From (53) we note that susceptible copepod population at endemic equilibrium is proportional to the average time of stay in susceptible copepod class and the rate at which new susceptible copepods are entering the susceptible copepod class through birth. The endemic value of infected copepod population is given by(54)IE∗=λE∗SE∗δE+μE=λE∗ΛEλE∗+μEδE+μE.We note from (54) that infected copepod population at the endemic equilibrium point is proportional to the average time of stay in the infected copepod class, the rate at which susceptible copepods become infected, and the density of susceptible copepods. We also make the endemic equilibrium of the model system (1) given by expressions (45)–(54) depend on both within-host and between-host disease parameters.


### 5.1. Existence of the Endemic Equilibrium State

In this section we present some results concerning the existence of an endemic equilibrium solution for the model system (1). To determine the existence and uniqueness of the endemic equilibrium point (EEP) of the model system (1), we can easily express SH∗, IH∗, IC∗, WM∗, WF∗, EW∗, and LW∗ in terms of IE∗ in the form(55)SH∗IE∗=ΛHa1+a2IE∗+αHa0IE∗P0+ϵIE∗a1+a2IH∗βHIE∗+μHP0+ϵIE∗,IH∗IE∗=a0IE∗a1+a2IE∗,IC∗IE∗=IE∗βHΛH-μHZEa+ZEbβHIE∗P0+ϵIE∗μCHGJ∗+αCIH∗+1ZEc,WM∗IE∗=NCμCGJ∗IE∗βHΛH-μHZEa+ZEbβHIE∗αM+μMP0+ϵIE∗μCGJ∗+αCIH∗+1ZEc,WF∗IE∗=12αMNCμCGJ∗IE∗βHΛH-μHZEa+ZEbβHIE∗αF+μFαM+μMP0+ϵIE∗μCHGJ∗+αCIH∗+1ZEc,EW∗IE∗=αFαMNCμCGJ∗IE∗βHΛH-μHZEa+ZEbβHIE∗2αW+μWαF+μFαM+μMP0+ϵIE∗μCHGJ∗+αCZEc,LW∗IE∗=QEGJ∗μCGJ∗+αC·βHΛH-μHIE∗ZEa+ZEbβHIE∗2P0+ϵIE∗ZEc,where(56)ZEa=a1+a2IE∗P0+ϵIE∗,ZEb=ΛHαHa0P0+ϵIE∗-a1+a2IE∗2βH,ZEc=a1+a2IE∗βHIE∗+μHP0+ϵIE∗,QE=12·NCμCμL·NWαWμW+αW·αFμF+αF·αMμM+αM,a0=βHΛH,a1=μHP0μH+δH+αH,a2=βHμH+δHμHϵμH+δH+αH.Substituting the expressionλE=βELW/L0+ϵLW and LW∗=QEGJ∗/(μCGJ∗+αC)·(βHΛH-μHIE∗ZEa+ZEbβHIE∗2/P0+ϵIE∗ZEc) into (25) we get(57)IE∗hIE∗=IE∗γ3IE∗3+γ2IE∗2+γ1IE∗+γ0=0,where(58)γ3=GJ∗μCG0+μJαCβE+ϵμEL0μEP0μHϵβEΛEG0μCGJ∗+αCR02a2βH-μH+a2βH+μHϵμEL0-βHΛHαHa0>0,γ1=GJ∗μCG0+μJαCβE+ϵμEL0μEμE+δHP02μHa1+ΛHαHa0βEΛEG0μCGJ∗+αCR02+a11-βHμE+δEμCG0+μJαCG0R02+P0μEL0a2βH+μHϵ-L0μEP0μHμCG0+μJαCGJ∗G0μCGJ∗+αCa1ϵ+a2P0+ΛHG0ϵΛH-μH,γ2=BβHΛH-μHa1ϵ+a2P0+a2μEL0βH+ϵμEΛH-μH+AβHβH-ΛH-μHϵ,γ0=μEL0μHa1P01-GJ∗μCG0+μJδCμCGJ∗+αCG0R02-GJ∗μCG0+μJαCL0μEP02μHβHΛH-μHG0μCGJ∗+αCμHμH+δH+αH-ΛH2αHR02,A=QEGJ∗βEΛEμCGJ∗+αCμE+δE,B=QEGJ∗βE+ϵμEμCGJ∗+αC.We can easily note that (57) gives IE∗=0, which corresponds to the disease-free equilibrium and(59)hIE∗=γ3IE∗3+γ2IE∗2+γ1IE∗+γ0=0,which corresponds to the existence of endemic equilibria. Solving forIE∗ in h(IE∗)=0, the roots of h(IE∗)=0 are determined by using Descartes’s rule of sign. The various possibilities are tabulated in Table 2.Table 2
Number of possible positive roots ofhIE∗=0.
Cases
γ
3
γ
2
γ
1
γ
0
Number of sign changes
Number of possible real roots (endemic equilibrium)
1
+
+
+
+
0
0
2
+
+
+
-
1
1
3
+
+
-
+
2
0, 2
4
+
+
-
-
1
1
5
+
-
-
+
2
0, 2
6
+
-
-
-
1
1
7
+
-
+
+
2
0, 2
8
+
-
+
-
3
1, 3We summarize the results in Table2 in the following Theorem 3.Theorem 3.
The model system (1)(1)
has a unique endemic equilibrium whenever Cases  1, 2, 3, 4, 5, 6, 7, and 8 are satisfied and ifR0>1,(2)
could have more than one endemic equilibrium if Case  8 is satisfied andR0>1,(3)
could have two endemic equilibria if Cases  3, 5, and 7 are satisfied.We now employ the center manifold theory [22] to establish the local asymptotic stability of the endemic equilibrium of model system (1).


### 5.2. Local Stability of the Endemic Equilibrium

We determine the local asymptotic stability of the endemic steady state of the model system (1) by using the center manifold theory described in [22]. In our case, we use center manifold theory by making the following change of variables. Let SH=x1, IH=x2, IC=x3, WM=x4, WF=x5, GJ=x6, EW=x7, EW=x8, SE=x9, and IE=x10. We also use the vector notation x=(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10)T so that the model system (1) can be written in the form(60)dxdt=fx,β∗,where(61)f=f1,f2,f3,f4,f5,f6,f7,f8,f9,f10.Therefore, model system (1) can be rewritten as(62)x˙1=ΛH-λHx1-μHx1+αHx2,x˙2=λHx1-μH+δH+αHx2,x˙3=λHx1-1x2+1-μCx6+αCx3,x˙4=NCμCx6x3-αM+μMx4,x˙5=αM2x4-μF+αFx5,x˙6=G0+αJx6x3-μJx6,x˙7=αFx5x2+1-μW+αWx7,x˙8=NWαWx7-μLx8,x˙9=ΛE-λEx9-μEx9,x˙10=λEx9-μE+δEx10,where(63)λH=β∗x10P0+ϵx10,λE=kβ∗x8L0+ϵx8.The method involves evaluating the Jacobian matrix of system (62) at the disease-free equilibrium E0 denoted by J(E0). The Jacobian matrix associated with the system of (62) evaluated at the disease-free equilibrium (E0) is given by (64)JE0=-μHαH0000000-βHΛHμHP00b00000000βHΛHP0μH00b1000000βHΛH-μHμHP000NCμCG0μJb2000000000αM2b30000000αJG0μJ00-μJ00000000αF0b4000000000NWαW-μL000000000-βEΛEμEL0-μE00000000βEΛEμEL00b5,where(65)b0=-μH+δH+αH,b1=-μCG0+μJαCμJ,b2=-μM+αM,b3=-μF+αF,b4=-μW+αW,b5=-μE+αE.By using the similar approach from Section4.1, the basic reproductive number of model system (62) is(66)R0=12·αMαM+μM·αFαF+μF·NCμCG0μCG0+μJαC·βHΛH-μHP0μH·NWαWαW+μWμLβEΛEμEμE+δEL0.Now let us considerβE=kβH, regardless of whether k∈(0,1) or k≥1, and let βH=β∗. Taking β∗ as the bifurcation parameter and if we consider R0=1 and solve for β∗ in (66), we obtain(67)β∗=2L0μE+δEμEμW+αWμLμM+αMμF+αFμCG0+αCμJP0μHkαFαMNCμCG0NWαWΛH-μHΛE.Note that the linearized system of the transformed equations (62) with bifurcation point β∗ has a simple zero eigenvalue. Hence, the center manifold theory [22] can be used to analyse the dynamics of (62) near βH=β∗.In particular, Theorem  4.1 in Castillo-Chavez and Song [23], reproduced below as Theorem 4 for convenience, will be used to show the local asymptotic stability of the endemic equilibrium point of (62) (which is the same as the endemic equilibrium point of the original system (1), for βH=β∗).Theorem 4.
Consider the following general system of ordinary differential equations with parameterϕ:(68)dxdt=fx,ϕ,f:Rn×R⟶R,f:C2R2×R,where 0 is an equilibrium of the system, that is, f(0,ϕ)=0 for all ϕ, and assume that(A1)
A=Dxf(0,0)=((∂fi/∂xj)(0,0)) is a linearization matrix of the model system (68) around the equilibrium 0 with ϕ evaluated at 0. Zero is a simple eigenvalue of A, and other eigenvalues of A have negative real parts,(A2)
matrixA has a right eigenvector u and a left eigenvector v corresponding to the zero eigenvalue.
Letfk be the kth component of f and(69)a=∑k,i,j=1nukvivj∂2fk∂xi∂xj0,0,b=∑k,i=1nukvi∂2fk∂xi∂ϕ0,0.
The local dynamics of (68) around 0 are totally governed by a and b and are summarized as follows.(i)
a>0 and b>0. When ϕ<0 with ϕ≪1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when 0<ϕ≪1, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium.(ii)
a<0 and b<0. When ϕ<0 with |ϕ|≪1, 0 is unstable; when 0<ϕ≪1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium.(iii)
a>0 and b<0. When ϕ<0 with |ϕ|≪1, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0<ϕ≪1, 0 is stable and a positive unstable equilibrium appears.(iv)
a<0 and b>0. When ϕ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.In order to apply Theorem4, the following computations are necessary (it should be noted that we are using β∗ as the bifurcation parameter, in place of ϕ in Theorem 4).Eigenvectors of
J
β
∗. For the case when R0=1, it can be shown that the Jacobian matrix of (62) at βH=β∗ (denoted by Jβ∗) has a right eigenvector associated with the zero eigenvalue given by(70)u=u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12T,where(71)u1=β∗ΛHμH2P0αHμH+δH+αH-1,u2=β∗ΛHμH+δH+αHP0μH,u3=β∗ΛH-μHP0μHμJμCG0+μJαC,u4=NCαCG0μCG0+μJαC·β∗ΛH-μHP0μHμM+αM,u5=αM2μM+αMμF+αFNCμCG0μCG0+μJαC·β∗ΛH-μHP0μH,u6=αJG0μJμCG0+μJαC·β∗ΛH-μHP0μH,u7=αMαF2μM+αMμF+αFNCμCG0μCG0+μJαC·β∗ΛH-μHP0μH1μW+αW,u8=αMαF2μM+αMμF+αFNCμCG0μCG0+μJαC·β∗ΛH-μHP0μHNWαWμLμW+αW,u9=-αMαF2μM+αMμF+αFNCμCG0μCG0+μJαC·β∗2ΛH-μHP0μHNWαWμLμW+αW·kΛEL0μE2.u10=1.In addition, the left eigenvector of the Jacobian matrix in (64) associated with the zero eigenvalue at βH=β∗ is given by(72)v=v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12T,where(73)v1=0,v2=0v3=1,v4=β∗2ΛH-μHμHP0·αFαM2μF+αFμM+αM·NWαWμLμW+αW·kΛEμEμE+δEL0,v5=β∗2ΛH-μHμHP0·αFμF+αF·NWαWμLμW+αW·kΛEμEμE+δEL0,v6=0,v7=β∗2ΛH-μHμHP0·NWαWμLμW+αW·kΛEμEμE+δEL0,v8=β∗2ΛH-μHμHP0·1μL·kΛEμEμE+δEL0,v9=0,v10=β∗ΛH-μHμE+δEμHP0.Computation of Bifurcation Parameters
a
and
b. We evaluate the nonzero second-order mixed derivatives off with respect to the variables and β∗ in order to determine the signs of a and b. The sign of a is associated with the following nonvanishing partial derivatives of f:(74)∂2f1∂x102=2ϵβ∗ΛHP02μH,∂2f2∂x102=-2ϵβ∗ΛHP02μH,∂2f3∂x102=-2ϵβ∗ΛH-μHP02μH,∂2f9∂x82=2ϵkβ∗ΛEL02μE,∂2f10∂x82=-2ϵkβ∗ΛEL02μE.The sign ofb is associated with the following nonvanishing partial derivatives off:(75)∂2f1∂x10∂β∗=-ΛHμHP0,∂2f2∂x10∂β∗=ΛHμHP0,∂2f3∂x10∂β∗=ΛH-μHμHP0,∂2f9∂x8∂β∗=-kΛEμEL0,∂2f10∂x8∂β∗=kΛEμEL0.Substituting expressions (71), (73), and (74) into (69), we get(76)a=u1v102∂2f1∂x102+u2v102∂2f2∂x102+u3v102∂2f3∂x102+u9v82∂2f9∂x82+u10v82∂2f10∂x82=u1v1022ϵβ∗ΛHP02μH+u2v102-2ϵβ∗ΛHP02μH+u3v102-2ϵβ∗ΛH-μHP02μH+u9v822ϵkβ∗ΛEL02μE+u10v82-2ϵkβ∗ΛEL02μE=2ϵβ∗ΛHP02μH·v102u1-u2-u3v1022ϵβ∗ΛH-μHP02μH+2ϵkβ∗ΛEL02μE·v82u9-u10<0since (u1-u2)<0, (u9-u10)<0, u3>0, and v10>0.Similarly, substituting expressions (71) and (73) and (75) into (69), we get(77)b=u1v10∂2f1∂x10∂β∗+u2v10∂2f2∂x10∂β∗+u3v10∂2f3∂x8∂β∗+u9v8∂2f9∂x10∂β∗+u10v8∂2f10∂x8∂β∗=v10ΛHP0μH·u2-ΛHP0μH·u1+ΛH-μHP0μH·u3+kΛEL0μE·v8u10-u9=ΛHP0μHv10u2-u1+ΛH-μHP0μHv10u3+kΛEL0μEv8u10-u9>0,since(u2-u1)>0, (u10-u9)>0, u3>0, and v10>0.Thus,a<0 and b>0. Using Theorem 4, item (iv), we have established the following result which only holds for R0>1 but close to 1.Theorem 5.
The endemic equilibrium guaranteed by Theorem3 is locally asymptotically stable for R0>1 near 1.


## 5.1. Existence of the Endemic Equilibrium State

In this section we present some results concerning the existence of an endemic equilibrium solution for the model system (1). To determine the existence and uniqueness of the endemic equilibrium point (EEP) of the model system (1), we can easily express SH∗, IH∗, IC∗, WM∗, WF∗, EW∗, and LW∗ in terms of IE∗ in the form(55)SH∗IE∗=ΛHa1+a2IE∗+αHa0IE∗P0+ϵIE∗a1+a2IH∗βHIE∗+μHP0+ϵIE∗,IH∗IE∗=a0IE∗a1+a2IE∗,IC∗IE∗=IE∗βHΛH-μHZEa+ZEbβHIE∗P0+ϵIE∗μCHGJ∗+αCIH∗+1ZEc,WM∗IE∗=NCμCGJ∗IE∗βHΛH-μHZEa+ZEbβHIE∗αM+μMP0+ϵIE∗μCGJ∗+αCIH∗+1ZEc,WF∗IE∗=12αMNCμCGJ∗IE∗βHΛH-μHZEa+ZEbβHIE∗αF+μFαM+μMP0+ϵIE∗μCHGJ∗+αCIH∗+1ZEc,EW∗IE∗=αFαMNCμCGJ∗IE∗βHΛH-μHZEa+ZEbβHIE∗2αW+μWαF+μFαM+μMP0+ϵIE∗μCHGJ∗+αCZEc,LW∗IE∗=QEGJ∗μCGJ∗+αC·βHΛH-μHIE∗ZEa+ZEbβHIE∗2P0+ϵIE∗ZEc,where(56)ZEa=a1+a2IE∗P0+ϵIE∗,ZEb=ΛHαHa0P0+ϵIE∗-a1+a2IE∗2βH,ZEc=a1+a2IE∗βHIE∗+μHP0+ϵIE∗,QE=12·NCμCμL·NWαWμW+αW·αFμF+αF·αMμM+αM,a0=βHΛH,a1=μHP0μH+δH+αH,a2=βHμH+δHμHϵμH+δH+αH.Substituting the expressionλE=βELW/L0+ϵLW and LW∗=QEGJ∗/(μCGJ∗+αC)·(βHΛH-μHIE∗ZEa+ZEbβHIE∗2/P0+ϵIE∗ZEc) into (25) we get(57)IE∗hIE∗=IE∗γ3IE∗3+γ2IE∗2+γ1IE∗+γ0=0,where(58)γ3=GJ∗μCG0+μJαCβE+ϵμEL0μEP0μHϵβEΛEG0μCGJ∗+αCR02a2βH-μH+a2βH+μHϵμEL0-βHΛHαHa0>0,γ1=GJ∗μCG0+μJαCβE+ϵμEL0μEμE+δHP02μHa1+ΛHαHa0βEΛEG0μCGJ∗+αCR02+a11-βHμE+δEμCG0+μJαCG0R02+P0μEL0a2βH+μHϵ-L0μEP0μHμCG0+μJαCGJ∗G0μCGJ∗+αCa1ϵ+a2P0+ΛHG0ϵΛH-μH,γ2=BβHΛH-μHa1ϵ+a2P0+a2μEL0βH+ϵμEΛH-μH+AβHβH-ΛH-μHϵ,γ0=μEL0μHa1P01-GJ∗μCG0+μJδCμCGJ∗+αCG0R02-GJ∗μCG0+μJαCL0μEP02μHβHΛH-μHG0μCGJ∗+αCμHμH+δH+αH-ΛH2αHR02,A=QEGJ∗βEΛEμCGJ∗+αCμE+δE,B=QEGJ∗βE+ϵμEμCGJ∗+αC.We can easily note that (57) gives IE∗=0, which corresponds to the disease-free equilibrium and(59)hIE∗=γ3IE∗3+γ2IE∗2+γ1IE∗+γ0=0,which corresponds to the existence of endemic equilibria. Solving forIE∗ in h(IE∗)=0, the roots of h(IE∗)=0 are determined by using Descartes’s rule of sign. The various possibilities are tabulated in Table 2.Table 2
Number of possible positive roots ofhIE∗=0.
Cases
γ
3
γ
2
γ
1
γ
0
Number of sign changes
Number of possible real roots (endemic equilibrium)
1
+
+
+
+
0
0
2
+
+
+
-
1
1
3
+
+
-
+
2
0, 2
4
+
+
-
-
1
1
5
+
-
-
+
2
0, 2
6
+
-
-
-
1
1
7
+
-
+
+
2
0, 2
8
+
-
+
-
3
1, 3We summarize the results in Table2 in the following Theorem 3.Theorem 3.
The model system (1)(1)
has a unique endemic equilibrium whenever Cases  1, 2, 3, 4, 5, 6, 7, and 8 are satisfied and ifR0>1,(2)
could have more than one endemic equilibrium if Case  8 is satisfied andR0>1,(3)
could have two endemic equilibria if Cases  3, 5, and 7 are satisfied.We now employ the center manifold theory [22] to establish the local asymptotic stability of the endemic equilibrium of model system (1).


## 5.2. Local Stability of the Endemic Equilibrium

We determine the local asymptotic stability of the endemic steady state of the model system (1) by using the center manifold theory described in [22]. In our case, we use center manifold theory by making the following change of variables. Let SH=x1, IH=x2, IC=x3, WM=x4, WF=x5, GJ=x6, EW=x7, EW=x8, SE=x9, and IE=x10. We also use the vector notation x=(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10)T so that the model system (1) can be written in the form(60)dxdt=fx,β∗,where(61)f=f1,f2,f3,f4,f5,f6,f7,f8,f9,f10.Therefore, model system (1) can be rewritten as(62)x˙1=ΛH-λHx1-μHx1+αHx2,x˙2=λHx1-μH+δH+αHx2,x˙3=λHx1-1x2+1-μCx6+αCx3,x˙4=NCμCx6x3-αM+μMx4,x˙5=αM2x4-μF+αFx5,x˙6=G0+αJx6x3-μJx6,x˙7=αFx5x2+1-μW+αWx7,x˙8=NWαWx7-μLx8,x˙9=ΛE-λEx9-μEx9,x˙10=λEx9-μE+δEx10,where(63)λH=β∗x10P0+ϵx10,λE=kβ∗x8L0+ϵx8.The method involves evaluating the Jacobian matrix of system (62) at the disease-free equilibrium E0 denoted by J(E0). The Jacobian matrix associated with the system of (62) evaluated at the disease-free equilibrium (E0) is given by (64)JE0=-μHαH0000000-βHΛHμHP00b00000000βHΛHP0μH00b1000000βHΛH-μHμHP000NCμCG0μJb2000000000αM2b30000000αJG0μJ00-μJ00000000αF0b4000000000NWαW-μL000000000-βEΛEμEL0-μE00000000βEΛEμEL00b5,where(65)b0=-μH+δH+αH,b1=-μCG0+μJαCμJ,b2=-μM+αM,b3=-μF+αF,b4=-μW+αW,b5=-μE+αE.By using the similar approach from Section4.1, the basic reproductive number of model system (62) is(66)R0=12·αMαM+μM·αFαF+μF·NCμCG0μCG0+μJαC·βHΛH-μHP0μH·NWαWαW+μWμLβEΛEμEμE+δEL0.Now let us considerβE=kβH, regardless of whether k∈(0,1) or k≥1, and let βH=β∗. Taking β∗ as the bifurcation parameter and if we consider R0=1 and solve for β∗ in (66), we obtain(67)β∗=2L0μE+δEμEμW+αWμLμM+αMμF+αFμCG0+αCμJP0μHkαFαMNCμCG0NWαWΛH-μHΛE.Note that the linearized system of the transformed equations (62) with bifurcation point β∗ has a simple zero eigenvalue. Hence, the center manifold theory [22] can be used to analyse the dynamics of (62) near βH=β∗.In particular, Theorem  4.1 in Castillo-Chavez and Song [23], reproduced below as Theorem 4 for convenience, will be used to show the local asymptotic stability of the endemic equilibrium point of (62) (which is the same as the endemic equilibrium point of the original system (1), for βH=β∗).Theorem 4.
Consider the following general system of ordinary differential equations with parameterϕ:(68)dxdt=fx,ϕ,f:Rn×R⟶R,f:C2R2×R,where 0 is an equilibrium of the system, that is, f(0,ϕ)=0 for all ϕ, and assume that(A1)
A=Dxf(0,0)=((∂fi/∂xj)(0,0)) is a linearization matrix of the model system (68) around the equilibrium 0 with ϕ evaluated at 0. Zero is a simple eigenvalue of A, and other eigenvalues of A have negative real parts,(A2)
matrixA has a right eigenvector u and a left eigenvector v corresponding to the zero eigenvalue.
Letfk be the kth component of f and(69)a=∑k,i,j=1nukvivj∂2fk∂xi∂xj0,0,b=∑k,i=1nukvi∂2fk∂xi∂ϕ0,0.
The local dynamics of (68) around 0 are totally governed by a and b and are summarized as follows.(i)
a>0 and b>0. When ϕ<0 with ϕ≪1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when 0<ϕ≪1, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium.(ii)
a<0 and b<0. When ϕ<0 with |ϕ|≪1, 0 is unstable; when 0<ϕ≪1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium.(iii)
a>0 and b<0. When ϕ<0 with |ϕ|≪1, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0<ϕ≪1, 0 is stable and a positive unstable equilibrium appears.(iv)
a<0 and b>0. When ϕ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.In order to apply Theorem4, the following computations are necessary (it should be noted that we are using β∗ as the bifurcation parameter, in place of ϕ in Theorem 4).Eigenvectors of
J
β
∗. For the case when R0=1, it can be shown that the Jacobian matrix of (62) at βH=β∗ (denoted by Jβ∗) has a right eigenvector associated with the zero eigenvalue given by(70)u=u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12T,where(71)u1=β∗ΛHμH2P0αHμH+δH+αH-1,u2=β∗ΛHμH+δH+αHP0μH,u3=β∗ΛH-μHP0μHμJμCG0+μJαC,u4=NCαCG0μCG0+μJαC·β∗ΛH-μHP0μHμM+αM,u5=αM2μM+αMμF+αFNCμCG0μCG0+μJαC·β∗ΛH-μHP0μH,u6=αJG0μJμCG0+μJαC·β∗ΛH-μHP0μH,u7=αMαF2μM+αMμF+αFNCμCG0μCG0+μJαC·β∗ΛH-μHP0μH1μW+αW,u8=αMαF2μM+αMμF+αFNCμCG0μCG0+μJαC·β∗ΛH-μHP0μHNWαWμLμW+αW,u9=-αMαF2μM+αMμF+αFNCμCG0μCG0+μJαC·β∗2ΛH-μHP0μHNWαWμLμW+αW·kΛEL0μE2.u10=1.In addition, the left eigenvector of the Jacobian matrix in (64) associated with the zero eigenvalue at βH=β∗ is given by(72)v=v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12T,where(73)v1=0,v2=0v3=1,v4=β∗2ΛH-μHμHP0·αFαM2μF+αFμM+αM·NWαWμLμW+αW·kΛEμEμE+δEL0,v5=β∗2ΛH-μHμHP0·αFμF+αF·NWαWμLμW+αW·kΛEμEμE+δEL0,v6=0,v7=β∗2ΛH-μHμHP0·NWαWμLμW+αW·kΛEμEμE+δEL0,v8=β∗2ΛH-μHμHP0·1μL·kΛEμEμE+δEL0,v9=0,v10=β∗ΛH-μHμE+δEμHP0.Computation of Bifurcation Parameters
a
and
b. We evaluate the nonzero second-order mixed derivatives off with respect to the variables and β∗ in order to determine the signs of a and b. The sign of a is associated with the following nonvanishing partial derivatives of f:(74)∂2f1∂x102=2ϵβ∗ΛHP02μH,∂2f2∂x102=-2ϵβ∗ΛHP02μH,∂2f3∂x102=-2ϵβ∗ΛH-μHP02μH,∂2f9∂x82=2ϵkβ∗ΛEL02μE,∂2f10∂x82=-2ϵkβ∗ΛEL02μE.The sign ofb is associated with the following nonvanishing partial derivatives off:(75)∂2f1∂x10∂β∗=-ΛHμHP0,∂2f2∂x10∂β∗=ΛHμHP0,∂2f3∂x10∂β∗=ΛH-μHμHP0,∂2f9∂x8∂β∗=-kΛEμEL0,∂2f10∂x8∂β∗=kΛEμEL0.Substituting expressions (71), (73), and (74) into (69), we get(76)a=u1v102∂2f1∂x102+u2v102∂2f2∂x102+u3v102∂2f3∂x102+u9v82∂2f9∂x82+u10v82∂2f10∂x82=u1v1022ϵβ∗ΛHP02μH+u2v102-2ϵβ∗ΛHP02μH+u3v102-2ϵβ∗ΛH-μHP02μH+u9v822ϵkβ∗ΛEL02μE+u10v82-2ϵkβ∗ΛEL02μE=2ϵβ∗ΛHP02μH·v102u1-u2-u3v1022ϵβ∗ΛH-μHP02μH+2ϵkβ∗ΛEL02μE·v82u9-u10<0since (u1-u2)<0, (u9-u10)<0, u3>0, and v10>0.Similarly, substituting expressions (71) and (73) and (75) into (69), we get(77)b=u1v10∂2f1∂x10∂β∗+u2v10∂2f2∂x10∂β∗+u3v10∂2f3∂x8∂β∗+u9v8∂2f9∂x10∂β∗+u10v8∂2f10∂x8∂β∗=v10ΛHP0μH·u2-ΛHP0μH·u1+ΛH-μHP0μH·u3+kΛEL0μE·v8u10-u9=ΛHP0μHv10u2-u1+ΛH-μHP0μHv10u3+kΛEL0μEv8u10-u9>0,since(u2-u1)>0, (u10-u9)>0, u3>0, and v10>0.Thus,a<0 and b>0. Using Theorem 4, item (iv), we have established the following result which only holds for R0>1 but close to 1.Theorem 5.
The endemic equilibrium guaranteed by Theorem3 is locally asymptotically stable for R0>1 near 1.


## 6. Sensitivity Analysis

In this section we carry out sensitivity analysis to evaluate the relative change in basic reproduction number(R0) when the within-host and between-host parameters as well as the environmental parameters of the model system (1) change. We used the normalized forward sensitivity index of the basic reproduction number, R0 of the model system (1) to each of the model parameters. The normalized forward sensitivity index of a variable to a parameter is typically defined as “the ratio of the relative change in the variable to the relative change in the parameter” [24]. In this case, if we let R0 be a differentiable function of the parameter u, then the normalized forward sensitivity index of R0 at u is defined as(78)ΥuR0=∂R0∂u×uR0,where the quotient u/R0 is introduced to normalize the coefficient by removing the effect of units [25]. For example, the sensitivity index of R0 with respect to the human infection rate βH is given by(79)ΥβHR0=∂R0∂βH×βHR0=0.5.It can be easily noted that the sensitivity index ofR0 with respect to the parameter βH does not depend on any of the parameter values. The indices of worm larvae death rate within a host and copepods death rate in the physical environment are, respectively, given by(80)ΥμFR0=-12μFμF+αF=-0.5,ΥμER0=-122μE+δEμE+δE=-0.9991.Using (78)–(80) similar expressions can be derived for the remaining parameters. The resulting sensitivity indices of R0 to the different model parameters are shown in Table 3. We see from (78)–(80) that the index of parameter βH is positive and indexes of both parameters μL and μE are negative. The sign of the index value indicates whether the parameter increases the reproduction number or reduces the reproduction number. Therefore increasing human infection rate βH reduces R0 and also increasing μL or μE reduces R0. Based on the results shown in Table 3, we observe that the reproduction number R0 is sensitive to the changes of both the within-host and between-host parameters as well as the environmental parameters (parameters which can be modified by environmental conditions which impact on survival and reproduction of the parasite and vector populations). More specifically, we deduce the following results for the between-host scale:(i)
The reproductive number is most sensitive to the changes of parameterμE and the natural death rate of copepods in the physical water environment. This implies that interventions focused on vector control have highest impact on GWD control. Since ΥμER0=-0.9991, increasing μE by 10% decreases the reproduction number by 9.991%. Therefore increasing the death rate of copepods by using chemical such as ABATE or temephos will eventually reduce the transmission of Guinea worm disease.(ii)
The reproductive number also shows significant sensitivity toβH and βE since ΥβHR0=ΥβER0=0.5. This implies that reducing human infection rates βH and βE by 10% reduces R0 by 5% for each of these parameters. Therefore, health education to ensure that greater numbers of individuals and communities adopt behavioural practices such as voluntary reporting of GWD cases, prevention of GWD patients from entering drinking water bodies, regular use of water from safe water sources, and, in the absence of such water sources, filtering or boiling water before drinking aimed at preventing transmission of GWD would have high impact in complementing vector control in elimination of GWD.(iii)
Similarly,ΥαFR0=0.4639. This implies that reducing the rate at which eggs are excreted in the physical water environment, αF, by 10% reduces R0 by 4.639%. Therefore educating people about GWD (i.e., teaching people not to immerse their infected feet into the drinking water when the fertilized female worm is emerging out from their feet or to always filter contaminated water before drinking the water) will reduce the transmission of the disease.Table 3
Sensitivity indices of model reproduction numberR0 to parameters for model system (1), evaluated at the parameters values presented in Tables 4–6.
Parameter
Description
Sensitivity index with positive sign
Sensitivity index with negative sign
μ
E
Natural decay rate of copepods in the water environment
−0.9991
α
C
Natural decay rate of copepods within human host
−0.4853
μ
C
Release rate of mature worms within human host
−0.4853
Λ
H
Human birth rate
+0.50013
β
H
Human infection rate
+0.5
N
W
Fecundity rate of worm larvae in the environment
+0.5
μ
L
Natural decay rate of Guinea worm larvae
−0.5
L
0
Larvae saturation constant
−0.5
N
C
Fecundity rate of mature worm
+0.5
Λ
E
Copepods birth rate
+0.5
β
E
Copepods infection rate
+0.5
P
0
Copepods saturated constant
−0.5
α
H
Human recovery
−0.4998
μ
F
Natural decay rate of fertilized female worms
−0.4639
α
F
Migration rate of fertilized female worms to surface of host’s skin
+0.4639
μ
W
Natural decay rate of worm eggs in the water environment
−0.4545
α
W
Worm egg hatching rate
+0.4545
μ
M
Natural decay rate of mature worms within human host
−0.25
α
M
Migration rate of mature worms to subcutaneous tissues
+0.25
μ
J
Dilution/degradation rate of gastric juice
−0.0147
G
0
Supply rate of gastric juice from the source of the body
+0.0147
δ
E
Induced decay rate of copepods in the water environment
−0.000894Table 4
Human host parameter values used in simulations.
Parameter
Description
Initial values
Units
Source
Λ
H
Human birth rate
0.1013
Peopleday-1
[14]
β
H
Human infection rate
0.1055
Copepodday-1
Estimated
μ
H
Human natural death rate
2.548
×
1
0
-
5
D
a
y
-
1
[15, 16]
α
H
Human recovery rate
0.03
D
a
y
-
1
Estimated
δ
H
GWD induced death rate
4
×
1
0
-
8
D
a
y
-
1
EstimatedTable 5
Within-host parameter values of the model system (1).
Parameter
Description
Initial values
Units
Source
N
C
Fecundity rate of mature worms
700
People
Estimated
μ
C
Decay rate of copepods within a human host due to gastric juice
0.99
Copepodday-1
Estimated
α
C
Natural death rate of copepods within a human host
0.001
D
a
y
-
1
Estimated
μ
M
Natural decay rate of mature worms within a human host
0.9
D
a
y
-
1
Estimated
α
M
Migration rate of mature worms to subcutaneous tissues
0.9
D
a
y
-
1
Estimated
μ
F
Natural death rate of fertilized female worms within a human host
0.9
D
a
y
-
1
Estimated
α
F
Migration rate of fertilized female worms to surface of skin
0.07
D
a
y
-
1
Estimated
μ
J
Dilution/degradation rate of gastric juice
0.05
D
a
y
-
1
Estimated
α
J
Proliferation rate of gastric juice due to infection
0.4
D
a
y
-
1
Estimated
G
0
Supply rate of gastric juice from within a human body
1.5
D
a
y
-
1
EstimatedTable 6
Free-living pathogens and their associated environmental parameter values used in simulations.
Parameter
Description
Initial value
Units
Source
Λ
E
Copepods birth rate
0.75
Copepodday-1
Estimated
β
E
Copepods infection rate
0.7
Larvaeday-1
Estimated
μ
E
Natural decay rate of copepods
0.005
D
a
y
-
1
[15, 16]
δ
E
Disease induced death rate of copepods
9
×
1
0
-
6
D
a
y
-
1
Estimated
P
0
Copepods saturation constant
20 0000
D
a
y
-
1
[15, 16]
μ
W
Natural decay rate of Guinea worm eggs
0.333
D
a
y
-
1
[15, 16]
α
W
Hatching rate of worm eggs
0.009
D
a
y
-
1
Estimated
N
W
Number of Guinea worm larvae hatched
300
Larvae egg-1day-1
Estimated
μ
L
Natural decay rate of Guinea worm larvae
0.0333
D
a
y
-
1
[15, 16]
L
0
Larvae saturation constant
5000000
D
a
y
-
1
[15, 16]
ϵ
Limitation growth rate
0.0991
D
a
y
-
1
EstimatedFurther, we also deduce the following results for the within-host scale:(i)
The development of a drug that would kill mature worms within human host would have significant benefits at within-host. However, the drug would have even higher impact if it would kill fertilized female worms.(ii)
The development of interventions that would increase the supply rate of gastric juice would have no benefits in the control of GWD.Therefore, the lack of drugs to treat GWD has delayed progress in eliminating GWD.


## 7. Numerical Analysis

The behaviour of model system (1) was investigated using numerical simulations using a Python program version V 2.6 on the Linux operation system (Ubuntu 14.04). The program uses a package odeint function in the scipy.integrate for solving a system of differential equations. The behaviour of the system model (1) was simulated in order to illustrate the analytical results we obtained in this paper. We used parameter values presented in Tables 4–6. Some of the parameter values used in the numerical simulations are from published literature while others were estimated as values of some parameters are generally not reported in literature. The initial conditions used for simulations are given by SH(0)=2500, IH(0)=10, IC(0)=0, GJ(0)=1.50, WM(0)=0, WF(0)=0, SE(0)=100000, IE(0)=0, EW(0)=0, and LW(0)=50000.Figure2 illustrates the solution profile of the population of (a) infected humans, (b) infected copepods in the physical water environment, (c) worm eggs in the physical water environment, and (d) worm larvae in the physical water environment, for different values of the infection rate of humans βH: βH=0.1055, βH=0.55, and βH=0.9. The numerical results show that higher rates of infection at the human population level result in increased population of parasites (worm eggs and worm larvae) in the physical water environment and a noticeable increase in infected copepod population in the physical water environment. Therefore, human behavioural changes which reduce contact with contaminated water bodies through drinking contaminated water reduce transmission of the disease at both human and copepod population level.Figure 2
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs in the physical water environment, and (d) population of Guinea worm larvae in the physical water environment, for different values of the infection rate of humans βH: βH=0.1055, βH=0.55, and βH=0.9.
(a)
(b)
(c)
(d)Figure3 illustrates the solution profile of the population of (a) infected humans, (b) infected copepods in the physical water environment, (c) worm eggs in the physical water environment, and (d) worm larvae in the physical water environment, for different values of natural death rate of copepod population in the physical water environment μE: μE=0.005, μE=0.05, and μE=0.5. The results show that environmental conditions which increase death of copepods affect transmission of the disease in the human population. Increased death of copepod population reduces transmission risk of the disease at humans population; therefore any mechanisms which enhance the killing of copepod population in the physical water environment reduces transmission risk of GWD within disease endemic communities.Figure 3
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs in the physical water environment, and (d) population of Guinea worm larvae in the physical water environment, for different values of natural death rate of copepods μE: μE=0.005, μE=0.05, and μE=0.5.
(a)
(b)
(c)
(d)Figure4 shows graphs of numerical solutions of model system (1) showing propagation of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs in the physical water environment, and (d) population of Guinea worm larvae in the physical water environment, for different values of natural death rate of Guinea worm eggs in the physical water environment μW: μW=0.005, μW=0.5, and μW=0.9. The results show that the environmental conditions which enhance death of worm eggs affect transmission of GWD in the human population. Increased death of worm egg population reduces transmission risk of the disease at human population level. Therefore any mechanisms which enhance the killing of worm egg population in the physical water environment reduce transmission risk of the disease within GWD endemic communities.Figure 4
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs in the physical water environment, and (d) population of Guinea worm larvae in the physical water environment, for different values of natural death rate of Guinea worm eggs in the physical water environment μW: μW=0.005, μW=0.5, and μW=0.9.
(a)
(b)
(c)
(d)Figure5 shows graphs of numerical solutions of model system (1) showing propagation of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs, and (d) population of Guinea worm larvae in the physical water environment, for different values of natural death rate of Guinea worm larvae in the physical water environment μL: μL=0.005, μL=0.5, and μL=0.9. The results show that the environmental conditions which increase death of worm larvae reduce transmission of GWD in the human population. Increased death of worm larvae population reduces transmission risk of the disease at human population level. Therefore any interventions which enhance the killing of worm larvae population in the physical water environment reduce transmission risk of GWD within the human population.Figure 5
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs, and (d) population of Guinea worm larvae in the physical water environment, for different values of natural death rate of Guinea worm larvae in the physical water environment μL: μL=0.005, μL=0.5, and μL=0.9.
(a)
(b)
(c)
(d)Figure6 shows graphs of numerical solution of model system (1) showing propagation of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs in the physical water environment, and (d) population of Guinea worm larvae in the physical water environment, for different values of natural death rate of mature worms inside a single infected human host μM: μM=0.9, μM=0.09, and μM=0.009. The results show that the within-host process of death of mature worms affects transmission of the disease in the human population level. Increased death of mature worm population within an infected human host reduces transmission risk of the disease at human population level. Therefore any interventions which enhance the killing of mature worm population inside an infected human host reduce transmission risk of the disease within communities.Figure 6
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs in the physical water environment, and (d) population of Guinea worm larvae in the physical water environment, for different values of natural death rate of mature worms inside a single infected human host μM: μM=0.9, μM=0.09, and μM=0.009.
(a)
(b)
(c)
(d)Figure7 shows graphs of numerical solution of model system (1) showing propagation of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs, and (d) population of Guinea worm larvae in the physical water environment, for different values of natural death rate of fertilized female worms inside a single infected human host μF: μF=0.9, μF=0.09, and μF=0.009. The results show that the within-host processes which increase the death of fertilized female worms can be a potent control measure for GWD. Increased death of fertilized female worms within infected human hosts reduces transmission risk of the disease at human population level. Therefore any interventions which enhance the killing of fertilized female worm population inside an infected human host reduce transmission risk of GWD within a community.Figure 7
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs, and (d) population of Guinea worm larvae in the physical water environment, for different values of natural death rate of fertilized female worm within a single infected human host μF: μF=0.9, μF=0.09, and μF=0.009.
(a)
(b)
(c)
(d)Figure8 illustrates the solution profiles of the population of (a) infected humans (IH), (b) infected copepods (IE) in the physical water environment, (c) worm eggs in the physical water environment, and (d) worm larvae in the physical water environment, for different values of the rate of worm larvae fecundity NW: NW=30, NW=300, and NW=30000. The results show that an increase of worm larvae produced per day by worm eggs increases the transmission of the disease. Therefore, any interventions which reduce worm larvae fecundity in the physical water environment reduce the transmission risk of the disease in the community.Figure 8
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs in the physical water environment, and (d) population of Guinea worm larvae in the physical water environment, for different values of Guinea worm larvae fecundity NW: NW=30, NW=300, and NW=30000.
(a)
(b)
(c)
(d)Figure9 shows graphs of numerical solution of model system (1) showing the propagation of the population of (a) infected humans (IH), (b) population of infected copepods (IE), (c) worm eggs in the physical water environment, and (d) worm larvae in the physical water environment, for different values of the rate at which an emerging fertilized female worm from a single infected human host excretes number of eggs into the physical water environment αF: αF=0.007, αF=0.07, and αF=0.7. The results show that higher rate of excretion of worm eggs by each infected human host results in increased population of parasites (worm eggs and worm larvae) in the physical water environment and a noticeable increase in infected copepods. Therefore, improvements in individual sanitation (which reduce contamination of water source with human eggs) are good for the community because they reduce the risk of the disease transmission in the community.Figure 9
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of infected humans (IH), (b) population of infected copepods (IE), (c) population of Guinea worm eggs in the physical water environment, and (d) population of Guinea worm larvae in the physical water environment, for different values of the rate at which an emerging fertilized female worm from a single infected human host excretes eggs into the physical water environment αF: αF=0.007, αF=0.07, and αF=0.7.
(a)
(b)
(c)
(d)Figure10 demonstrates numerical solutions showing the propagation of the population of (a) mature worm within infected human host and (b) population of fertilized female worm within infected human host, for different values of the infection rate of humans βH: βH=0.1055, βH=0.01055, βH=0.001055, and βH=0.55. The results show the influence of between-host disease process on within-host disease process of Guinea worm disease. Here, as transmission rate of GWD in the community increases, the within-host infection intensity of the disease also increases. The numerical results demonstrate that the transmission of the disease at the population level influences the dynamics within an infected individual. Therefore, human behavioural changes (such as filtering water before drinking) which reduce contact with infected copepods reduce infection intensity at individual level. Equally, good sanitation by community which reduces contamination of water bodies reduces the intensity of infection of humans at individual level.Figure 10
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of mature worm within infected human host and (b) population of fertilized female worm within infected human host, for different values of the infection rate of humans βH: βH=0.1055, βH=0.01055, βH=0.001055, and βH=0.55.
(a)
(b)Figure11 illustrates the graphs of numerical solutions showing the propagation of the population of (a) mature worm within infected human host and (b) population of fertilized female worm within infected human host, for different values of the natural death rate of copepods in the physical water environment μE: μE=0.005, μE=0.05, and μE=0.5. The results demonstrate the potency of public health interventions intended to reduce copepods population (such as killing copepods using temephos or boiling the water) on the infection intensity within an infected individual.Figure 11
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of mature worm within infected human host and (b) population of fertilized female worm within infected human host, for different values of the natural death rate of copepods μE: μE=0.005, μE=0.05, and μE=0.5.
(a)
(b)Figure12 illustrates the graphs of numerical solutions showing the propagation of the population of (a) mature worm within infected human host and (b) population of fertilized female worm within infected human host, for different values of the natural death rate of worm larvae in the physical water environment μL: μL=0.005, μL=0.05, and μL=0.5. The results demonstrate the influence of public health interventions intended to reduce worm larvae population in the physical water environment (such as destroying worm larvae using chemicals or boiling the water) on the infection intensity within an infected individual. Overall, the numerical results verify the following aspects about GWD transmission dynamics.(i)
Higher rates of infection at the human population level result in increased population of parasites (worm eggs and worm larvae) in the physical water environment and a noticeable increase in infected copepod population in the physical water environment.(ii)
Interventions which increase death of copepods through enhanced killing of copepod population in the physical water environment reduce transmission risk of GWD within a disease endemic communities.(iii)
Interventions which enhance death of worm eggs through enhanced killing of worm egg population in the physical water environment reduce transmission risk of the disease within GWD endemic communities.(iv)
Health interventions which increase death of worm larvae in the physical water environment reduce transmission risk of GWD within the human population.(v)
Within-host scale interventions which increase death of mature worm population inside an infected human host reduce transmission risk of the disease within communities.(vi)
Within-host scale interventions which increase the death of fertilized female worms can also be a potent control measure for GWD through reduced transmission risk of GWD within a community.(vii)
An increase in worm fecundity with an infected human host has considerable impact on the transmission of GWD.(viii)
Higher rate of contamination of the physical water environment through excretion of worm eggs by each infected human host results in increased population of parasites (worm eggs and worm larvae) in the physical water environment and a noticeable increase in infected copepods.(ix)
Transmission of the GWD shows reciprocal influence of within-host scale interventions (medical interventions) and between-host scale interventions (public health interventions). Therefore, human behavioural changes (such as filtering water before drinking) which reduce contact with infected copepods reduce infection intensity at individual level. Equally, good sanitation by community which reduces contamination of water bodies reduces the intensity of infection of humans at individual level.Figure 12
Graphs of numerical solutions of model system (1) showing the evolution with time of (a) population of mature worm within infected human host and (b) population of fertilized female worm within infected human host, for different values of natural death rate of worm larvae in the physical water environment μL: μL=0.005, μL=0.05, and μL=0.5.
(a)
(b)


## 8. Conclusions

Guinea worm disease, like most neglected parasitic diseases, urgently needs renewed attention and sustainable interventions in Africa. The limited scientific knowledge about GWD represents a challenge to the successful elimination of the disease. In this paper, we have sought to identify a broad range of within-host and between-host processes that should be better understood if GWD is to be eliminated. A multiscale model of GWD transmission dynamics is presented. The multiscale model is shown to be mathematically and epidemiologically well-posed. Sensitivity analyses of the basic reproduction number to the variation of model parameters were carried out. The sensitivity results of the reproduction number show that between-host model parameters (such as infection rate of human hostβH and supply rate of humans ΛH); within-host model parameters (such as excretion rate of eggs αF into the physical water environment by each infected human host, fecundity rate of mature worm NC, decay rate of mature worms μM, migration rate of mature worms to the subcutaneous tissue αM, and decay rate of fertilized worms μF); and environmental model parameters (such as the production rate of larvae per egg worm per day αW, fecundity of larvae NW generated by eggs, death rate of egg worms μW, larva worms μL in the physical water environment, supply rate of copepods ΛE, and decay rate of copepods μE) all are responsible for the transmission dynamics of Guinea worm disease within the community. Therefore reducing the infection rate of human, excretion rate of eggs into physical water, and the population of parasites (worm eggs and worm larvae) as well as population of vector host (copepods) in the physical water environment could eventually contribute in eradicating GWD completely from the community.



---

*Source: 1473287-2017-07-20.xml*