Modeling Schistosomiasis and HIV/AIDS Codynamics

Abstract

We formulate a mathematical model for the cointeraction of schistosomiasis and HIV/AIDS in order to assess their synergistic relationship in the presence of therapeutic measures. Comprehensive mathematical techniques are used to analyze the model steady states. The disease-free equilibrium is shown to be locally asymptotically stable when the associated disease threshold parameter known as the basic reproduction number for the model is less than unity. Centre manifold theory is used to show that the schistosomiasis-only and HIV/AIDS-only endemic equilibria are locally asymptotically stable when the associated reproduction numbers are greater than unity. The impact of schistosomiasis and its treatment on the dynamics of HIV/AIDS is also investigated. To illustrate the analytical results, numerical simulations using a set of reasonable parameter values are provided, and the results suggest that schistosomiasis treatment will always have a positive impact on the control of HIV/AIDS.

1. Introduction

Schistosomiasis, also known as bilharzia after Theodor Bilharz who first identified the parasite in Egypt in 1851, is a disease caused by blood flukes [1]. It affects millions of people worldwide, especially in South America, the Middle East, and Southeast Asia where it remains a public health problem and poses a threat to 600 million people in more than 76 countries [1]. The disease is often associated with water resource development projects, such as dams and irrigation schemes, where the snail intermediate hosts of the parasite breed [2]. Human schistosomiasis (which has a relatively low mortality rate, but a high morbidity rate) is a family of diseases primarily caused by three species of the genus Schistosoma or flat worms. The adult worms inhabit the blood vessels lining either the intestine or bladder, depending on the species of the worm [3]. The highest number of human schistosomiasis infections is caused by S. haematobium, which has a predilection for the blood vessels around the bladder and causes urinary disease [4]. Schistosomiasis is the second most prevalent neglected tropical diseases after hookworm (192 million cases), accounting for 93% of the world's number of cases and possibly associated with increased horizontal transmission of HIV/AIDS [5].

On the other hand, the number of people living with HIV worldwide continued to grow in 2008, reaching an estimated 33.4 million, which is more than 20% higher than the number in 2000, and the prevalence was roughly threefold higher than in 1990 [6]. The HIV virus, by holding the immune system hostage, has opened many gates for pathological interactions with other diseases [7]. Schistosomiasis and HIV infections have major effects on the host immune response, and coinfection (of the two diseases which may increase the complexity of treatment for people living with HIV and may contribute to poorer medical outcomes) is widespread [8]. While schistosomiasis infections are caused by diverse species from three phyla, HIV is essentially a single entity. There is some evidence that schistosomiasis infection provides some benefit in some instances like the atopic disease [9, 10], and the inflammatory pathology of autoimmune disease [11–13]. For bacterial and viral infections, impaired control of replication and elimination may lead to a detrimental outcome [14–17]. That HIV infection is detrimental to the immune response to many pathogens is quite clear and poor regulation of immune system in advanced HIV infection is illustrated by an increased incidence of hypersensitive drug reactions [18, 19]. Studies that examine the codynamics of HIV and schistosomiasis infections have shown a significant association between HIV and the presence of S. haematobium eggs in the genital samples, supporting the argument that schistosomiasis infection enhances HIV susceptibility when genital lesions are present [20]. Host-parasite interactions such as schistosomiasis, where inflammatory responses have persisted through evolution, perhaps due to selective advantage for parasite egg excretion, may be more detrimental with regard to HIV infection [21].

Although the negative impact of the synergetic interactions between HIV and schistosomiasis has shown to be a public health burden, only few statistical or mathematical models have been used to explore the consequences of their joint dynamics at the population level. There are plenty of single disease dynamic models. A significant number focus on HIV/AIDS [22–26] or on the transmission dynamics of schistosomiasis [27–37]. Schistosomiasis model (24) considered in this study differs from those found in the literature in that we consider Schistosoma mansoni a human blood fluke which causes schistosomiasis and is the most widespread and the fresh water snail Biomphalaria glabrata serves as the main intermediate host, while the HIV/AIDS model (7) is an extension of the model by Murray [38] by including HIV therapy while neglecting the issue of seropositivity considered in [38]. Mathematical modeling assessing the impact of schistosomiasis on the transmission dynamics of HIV/AIDS is rare [39].

Quantifying by how much treatment of schistosomiasis affects HIV/AIDS dynamics will require an extensive sensitivity analysis with parameter values estimated from real and recent coinfection data. Nevertheless, this theoretical study provides a framework for the potential benefit of schistosomiasis treatment on the dynamics of HIV and highlights the fact that global public health challenges require comprehensive and multipronged approaches to dealing with coinfections [7], and current intervention efforts that focus on a single infection at a time may be losing substantial rewards of dealing synergistically and concurrently with multiple infectious diseases in one host. To the best of our knowledge, except for the study in [39] where the co-interaction of schistosomiasis and HIV without any form of treatment is investigated, this work is possibly the first to give a theoretical mathematical account of the impact of schistosomiasis on HIV dynamics in the presence of both schistosomiasis treatment and antiretroviral therapy at the population level.

The rest of the paper is structured as follows. In the next section, we present the schistosomiasis and HIV/AIDS coinfection model. In Section 3 we determine sufficient conditions for local stability of the disease-free and endemic equilibria and analyze the reproduction number for the two diseases separately while Section 4 provides a comprehensive analysis of the full model. Section 5 provides numerical results while Section 6 concludes the paper.

2. Model Description

The proposed model is an extension of an earlier study [39], which did not account for any intervention strategy. The schistosomiasis and HIV models will be coupled via the force of infection, and in the absence of any of the diseases (hence no coinfection), the two basic disease submodels can be decoupled from the general model (see Sections 3.1.4 and 3.3.1). The population of interest is divided into several compartments dictated by the epidemiological stages (disease status), namely, susceptibles 𝑆_𝐻(𝑡), who are not yet infected by either HIV or schistosomiasis, schistosomiasis-infected individuals 𝐼_𝐵(𝑡), HIV-infected individuals not yet displaying symptoms of AIDS 𝐼_𝐻(𝑡), individuals infected with HIV showing symptoms of AIDS 𝐴_𝐻(𝑡), individuals dually infected with schistosomiasis and HIV displaying symptoms of schistosomiasis only 𝐼_𝐻𝐵(𝑡), individuals dually infected with schistosomiasis and HIV displaying symptoms of schistosomiasis and AIDS 𝐴_𝐻𝐵(𝑡), treated individuals infected with HIV only, showing symptoms of AIDS 𝐴_HT𝐴(𝑡), and treated individuals dually infected with schistosomiasis and HIV displaying symptoms of schistosomiasis and AIDS 𝐴_𝐻𝑇𝐵(𝑡). Other important populations to consider in this model are the susceptible snails 𝑆_𝑠(𝑡), infected snails 𝐼_𝑠(𝑡), miracidia population 𝑀(𝑡), and the cercariae population 𝑃(𝑡). Individuals move from one class to the next as the disease progresses and/or through dual infection. We further make the following assumptions for the model. (i)There is no vertical transmission of both infections in humans. (ii)Infected snails do not reproduce due to castration by miracidia. (iii)Seasonal and weather variations do not affect snail populations and contact patterns. (iv)Susceptible humans become infected with schistosomiasis only through contact with free-living pathogen in infested waters.

At any time, new recruits enter the human and snail populations through birth/migration at constant rates Λ_𝐻 and Λ_𝑆, respectively. There is a constant natural death rate 𝜇_𝐻 in each human subclass. The force of infection associated with HIV infection, denoted by 𝜆_𝐻, is given by𝜆_𝐻(𝑡)=𝛽^𝐻𝑐𝐼_𝐻+𝐼_𝐻𝐵+𝜂𝐴_𝐻+𝐴_𝐻𝐵+𝜅𝐴_𝐻𝑇𝐴+𝐴_𝐻𝑇𝐵𝑁_𝐻,(1) with 𝛽_𝐻 being the probability of HIV transmission per sexual contact, 𝑐 is the effective contact rate for HIV infection to occur, and 𝜂>1 models the fact that individuals in the AIDS stage and not on antiretroviral therapy are more infectious since the viral load is correlated with infectiousness [42]. It is assumed that individuals on antiretroviral therapy transmit infection at the smallest rate 𝜅 (with 0<𝜅<1) because of the fact that these individuals have very small viral load. It has been estimated by an analysis of longitudinal cohort data that antiretroviral therapy reduces per-partnership infectivity by as much as 60% (so that 𝜅=0.4) [41]. Thus, the total human population 𝑁_𝐻(𝑡) is given by𝑁_𝐻(𝑡)=+𝑆_𝐻𝐼_B(𝑡)+𝐼_𝐻(𝑡)+𝐴_𝐻(𝑡)+𝐴_𝐻𝑇𝐴(𝑡)(𝑡)+𝐼_𝐻𝐵(𝑡)+𝐴_𝐻𝐵(𝑡)+𝐴_𝐻𝑇𝐵(𝑡).(2) Susceptible individuals acquire schistosomiasis following infection at a rate 𝜆_𝑃, where𝜆_𝑃=𝛽^𝑃𝑃(𝑡)𝑃₀+𝜖𝑃(𝑡),(3) with 𝛽_𝑃 being the maximum rate of exposure, 𝜖 is the limitation of the growth velocity of cercariae with the increase of cases, and 𝑃₀ is the half saturation constant. In the absence of the parasite, the functional response of individuals susceptible to the pathogen (schistosomiasis) is given by (𝜆_𝑃/𝛽_𝑃)[𝑆_𝐻(𝑡)+𝐼_𝐻(𝑡)+𝐴_𝐻(𝑡)+𝐴_𝐻𝑇𝐴(𝑡)], a modified Holling's type-II functional response (also known as the Michaelis-Menten function when 𝜖=1), the response refers to the change in the density of susceptibles per unit time per pathogen as the schistosomiasis susceptible population density changes. From the functional response, we note that at low parasite density, contacts are directly proportional to host density, but a maximum rate of contact is reached at very high densities (saturation incidence). Individuals infected with schistosomiasis have an additional disease-induced death rate 𝑑_𝐵. Similarly, susceptible and infected snails have a natural death rate 𝜇_𝑆, and the infected snails have an additional disease-induced death rate 𝑑_𝑆. The total snail population is given by 𝑁_𝑆(𝑡)=𝑆_𝑆(𝑡)+𝐼_𝑆(𝑡).

Considering a schistosomiasis-infected individual, a number (portion) 𝑁_𝐸 of eggs leave the body through excretion (faeces and urine) and find their way into the fresh water supply where they hatch into free swimming ciliated miracidium at a rate 𝛾 for individuals without AIDS. Given the weakened immune system of AIDS individuals, they tend to excrete more often, thus releasing more eggs which will hatch into miracidia at a rate 𝜎𝛾,𝜎>1. If the miracidium reaches a fresh water with snails of a suitable species, it penetrates at a rate 𝜆_𝑀, where𝜆_𝑀=𝛽^𝑀𝑀(𝑡)𝑀₀+𝜖𝑀(𝑡),(4) and transforms into a sporocyst otherwise, the miracidia die naturally at a rate 𝜇_𝑀. The infected snails release a second form of free swimming larva called a cercariae which is capable of infecting humans at rate 𝜃. Some cercariae also die naturally at a rate 𝜇_𝑃. Individuals infected with schistosomiasis are infected with HIV at a rate 𝛿𝜆_𝐻 with 𝛿>1 since infection by schistosomiasis creates wounds within the urethra as eggs are being released, which increases the likelihood of HIV infection per sexual contact. Individuals with HIV progress to the AIDS stage at a rate 𝜌. Individuals in the AIDS stage have an additional disease-induced death rate 𝑑_𝐴. We assume that antiretroviral therapy is given to AIDS individuals who are ill and have experienced AIDS-defining symptoms, or whose CD4+ T cell count is below 200/𝜇L, which is the recommended AIDS defining stage [42]. Thus, AIDS patients are assumed to get antiretroviral therapy at a constant rate 𝛼. Treated AIDS patients eventually succumb to AIDS-induced mortality at a reduced rate modeled by the parameter 𝜏(0<𝜏<1). Individuals treated for schistosomiasis are assumed to recover at a constant rate 𝜔, and 𝜔₁ denotes AIDS patients who have recovered from schistosomiasis but are on antiretroviral therapy since the latter is a life treatment. The model flowchart for the interaction of the two diseases is shown in Figure 1 and parameters described will assume values in Table 1.

Table 1 

Model parameters and their interpretations.

Figure 1 

Model flow diagram.

From the aforementioned model description and assumptions, we establish the following deterministic system of nonlinear differential equations Modelsystem⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩𝑑𝑆^𝐻𝑑𝑡=Λ^𝐻𝑑𝐼_𝐵+𝜔𝐼_𝐵−𝜆_𝐻+𝜆_𝑃𝑆_𝐻−𝜇_𝐻𝑆_𝐻,𝑑𝑡=𝜆^𝑃𝑑𝐼_𝐻𝑆_𝐻−𝛿𝜆_𝐻𝐼_𝐵−𝜇_𝐻+𝜔+𝑑_𝐵𝐼_𝐵,𝑑𝑡=𝜆^𝐻𝑑𝐴_𝐻𝑆_𝐻+𝜔𝐼_𝐻𝐵−𝜆_𝑃𝐼_𝐻−𝜇_𝐻+𝜌𝐼_𝐻,𝑑𝑡=𝜌𝐼^𝐻𝑑𝐴_𝐻𝑇𝐴+𝜔𝐴_𝐻𝐵−𝜆_𝑃𝐴_𝐻−𝜇_𝐻+𝛼+𝑑_𝐴𝐴_𝐻,𝑑𝑡=𝛼𝐴^𝐻+𝜔𝐴_𝐻𝑇𝐵−𝜆_𝑝𝐴_𝐻𝑇𝐴𝑑𝐼_𝐻𝐵−𝜇_𝐻+𝜏𝑑_𝐴𝐴_𝐻𝑇𝐴,𝑑𝑡=𝛿𝜆^𝐻𝑑𝐴_𝐻𝐵𝐼_𝐵+𝜆_𝑃𝐼_𝐻−𝜌+𝜔+𝜇_𝐻+𝑑_𝐵𝐼_𝐻𝐵,𝑑𝑡=𝜆^𝑃𝑑𝐴_𝐻𝑇𝐵𝐴_𝐻+𝜌𝐼_𝐻𝐵−𝜇_𝐻+𝜔+𝑑_𝐴+𝑑_𝐵𝐴_𝐻𝐵,𝑑𝑡=𝜆^𝑃𝐴_𝐻𝑇𝐴+𝛼𝐴_𝐻𝐵𝑑𝑀−𝜇_𝐻+𝜔+𝜏𝑑_𝐴+𝑑_𝐵𝐴_𝐻𝑇𝐵,𝑑𝑡=𝑁^𝐸𝑑𝑆_𝑆𝛾𝐼_𝐵+𝐼_𝐻𝐵+𝜎𝐴_𝐻𝐵+𝜎𝐴_𝐻𝑇𝐵−𝜇^𝑀𝑀,𝑑𝑡=Λ^𝑆𝑑𝐼_𝑆−𝜆_𝑀𝑆_𝑆−𝜇_𝑆𝑆_𝑆,𝑑𝑡=𝜆^𝑀𝑑𝑃𝑆_𝑆−𝜇_𝑆+𝑑_𝑆𝐼_𝑆,𝑑𝑡=𝜃𝐼^𝑆−𝜇_𝑃𝑃.(5)

2.1. Model Basic Properties

In this section, we study the basic properties of the solutions of model system (5), which are essential in the proofs of stability.

Lemma 1. The equations preserve positivity of solutions.

Proof. Considering the human population only, the vector field given by the right-hand side of (5) points inward on the boundary of ℝ⁸⁺⧵{0}. For example, if 𝐴_𝐻=0, then, 𝐴_𝐻′=𝜌𝐼_𝐻+𝜔𝐴_𝐻𝐵≥0. In an analogous manner, the same result can be shown for the other model components (variables). We shall use the human population to illustrate the boundedness of solutions for model system (5).

Lemma 2. Each nonnegative solution of model system (5) is bounded in 𝐿¹-norm.

Proof. Consider the human population only, and let 𝐿^1𝐻∈𝐿¹; then, the norm 𝐿^1𝐻 of each nonnegative solution in 𝑁_𝐻 is given by max{𝑁_𝐻(0),Λ_𝐻/𝜇_𝐻}. Thus, the norm 𝐿^1𝐻 satisfies the inequality 𝑁^𝐻≤Λ−𝜇_𝐻𝑁_𝐻. Solutions to the equation 𝑄^=Λ−𝜇𝑄 are monotone increasing and bounded by Λ/𝜇 if 𝑄(0)<Λ/𝜇. They are monotone decreasing and bounded above if 𝑄(0)≥Λ/𝜇. Since 𝑁^𝐻≤𝑄^, the claim follows and in a similar fashion, the remaining model variables can be shown to bounded.

Corollary 1. The region Φ=⎧⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎩𝑆^𝐻∈ℝ,𝐼_𝐵8+,𝐼_𝐻,𝐴_𝐻,𝐴_𝐻𝑇𝐴,𝐼_𝐻𝐵,𝐴_𝐻𝐵,𝐴_𝐻𝑇𝐵∶𝑁_𝐻≤Λ^𝐻𝜇_𝐻𝑀∈ℝ₊,∶𝑀≤𝛾Λ^𝐻𝑁_𝐸(1+𝜎)𝜇_𝑀𝜇_𝐻𝑆_𝑆,,𝐼_𝑆∈ℝ²⁺∶𝑁_𝑆≤Λ^𝑆𝜇_𝑆𝑃∈ℝ₊,∶𝑃≤𝜃Λ^𝑆𝜇_𝑃𝜇_𝑆(6) is invariant and attracting for system (5).

Theorem 1. For every nonzero, nonnegative initial value, solutions of model system (5) exist for all time 𝑡>0.

Proof. Local existence of solutions follows from standard arguments since the right-hand side of (5) is locally Lipschitz. Global existence follows from the a priori bounds.

3. Analysis of the Submodels

Before analyzing the full model system (5), it is essential to gain insights into the dynamics of the models for HIV only and schistosomiasis only.

3.1. HIV-Only Model

We now consider a model for HIV/AIDS only, obtained by setting 𝐼_𝐵=𝐼_𝐻𝐵=𝐴_𝐻𝐵=𝐴_𝐻𝑇𝐵=𝑀=𝑆_𝑆=𝐼_𝑆=𝑃=0, so that system in (5) reduces toHIV/AIDSonly⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩𝑑𝑆^𝐻𝑑𝑡=Λ^𝐻𝑑𝐼_𝐻−𝜆_𝐻+𝜇_𝐻𝑆_𝐻,𝑑𝑡=𝜆^𝐻𝑑𝐴_𝐻𝑆_𝐻−𝜌+𝜇_𝐻𝐼_𝐻,𝑑𝑡=𝜌𝐼^𝐻𝑑𝐴_𝐻𝑇𝐴−𝛼+𝑑_𝐴+𝜇_𝐻𝐴_𝐻,𝑑𝑡=𝛼𝐴^𝐻with,𝜆_𝐻−𝜏𝑑_𝐴+𝜇_𝐻𝐴_𝐻𝑇𝐴,=𝛽^𝐻𝑐𝐼_𝐻+𝜂𝐴_𝐻+𝜅𝐴_𝐻𝑇𝐴𝑁_𝐻𝑁_𝐻,=𝑆_𝐻+𝐼_𝐻+𝐴_𝐻+𝐴_𝐻𝑇𝐴.(7) For system (7), it can be shown that the region Φ_𝐻=𝑆_𝐻,𝐼_𝐻,𝐴_𝐻,𝐴_𝐻𝑇𝐴∈ℝ⁴⁺∶𝑁_𝐻≤Λ^𝐻𝜇_𝐻(8) is invariant and attracting. Thus, the dynamics of the HIV-only model will be considered in Φ_𝐻.

3.1.1. Disease-Free Equilibrium and Stability Analysis

Model system (7) has an evident disease-free given by𝒰_0𝐻=𝑆^0𝐻,𝐼^0𝐻,𝐴^0𝐻,𝐴^0𝐻_𝑇𝐴=Λ^𝐻𝜇_𝐻,0,0,0.(9) Following the next generation approach and the notation defined therein [43], matrices 𝐹 and 𝑉 for new infection terms and the remaining transfer terms are, respectively, given by𝐹=⎡⎢⎢⎣𝛽^𝐻𝑉=⎡⎢⎢⎣𝜇_H00𝑐𝛽_𝐻00𝑐𝜂𝛽_𝐻00𝑐𝜅⎤⎥⎥⎦,0𝜇_𝐻0000+𝜌𝜇_𝐻0+𝜏𝑑_𝐴⎤⎥⎥⎦.(10) It follows from (10) that the reproduction number of the system (7) is given byℛ_𝐴=𝛽^𝐻𝑐𝜌𝜅𝛼+𝜇_𝐻+𝜏𝑑_𝐴𝜂𝜌+𝛼+𝜇_𝐻+𝑑_𝐴𝜇_𝐻+𝜌𝜇_𝐻+𝑑_𝐴𝜇_𝐻+𝛼+𝑑_𝐴.(11) The threshold quantity ℛ_𝐴 measures the average number of new secondary cases generated by a single individual in a population where the aforementioned HIV control measures are in place. An associated epidemiological threshold which is the basic reproductive number ℛ₀, obtained using the same technique of the next generation operator [43], by considering model system (7) in the absence of HIV intervention strategies, is given byℛ_0𝐴=𝛽^𝐻𝑐𝜇_𝐻+𝑑_𝐴+𝜂𝜌𝜇_𝐻+𝜌𝜇_𝐻+𝑑_𝐴.(12) This disease threshold quantity ℛ_0𝐴 measures the average number of new infections generated by a single infected individual in a completely susceptible population where there are no HIV intervention strategies. Using Theorem 2 in [43], the following result is established.

Lemma 3. The disease-free equilibrium 𝒰_0𝐻 of system (7) is locally asymptotically stable (LAS) if ℛ_𝐴<1 and unstable if ℛ_𝐴>1.

3.1.2. Sensitivity Analysis of HIV-Only-Induced Reproductive Number

To avoid repetition we refer the reader to a detailed analysis of the reproductive number for model system (7), in the work of Bhunu et al. [44].

3.1.3. Global Stability of HIV/AIDS Model

We claim the following result.

Lemma 4. The disease-free equilibrium (𝒰_0𝐻) of model system (7) is globally asymptotically stable (GAS) if ℛ_𝐴<1 and unstable if ℛ_𝐴>1.

Proof. The proof is based on using a comparison theorem [45]. Note that the equations of the infected components in system (7) can be written as ⎡⎢⎢⎢⎢⎢⎣𝑑𝐼^𝐻𝑑𝐴𝑑𝑡_𝐻𝑑𝐴𝑑𝑡_𝐻𝑇𝐴𝑑𝑡⎤⎥⎥⎥⎥⎥⎦=[𝐹−𝑉]⎡⎢⎢⎣𝐼^𝐻𝐴_𝐻𝐴_𝐻𝑇𝐴−𝛽_𝐻⎤⎥⎥⎦𝑐1−𝑆^𝐻𝑁_𝐻⎡⎢⎢⎣100𝜂00𝜅00⎤⎥⎥⎦⎡⎢⎢⎣𝐼^𝐻𝐴_𝐻𝐴_𝐻𝑇𝐴⎤⎥⎥⎦,(13) where 𝐹 and 𝑉, are as defined earlier in (10). Since 𝑆_𝐻≤𝑁_𝐻, (for all 𝑡≥0) in Φ_𝐻, it follows that ⎡⎢⎢⎢⎢⎢⎣𝑑𝐼^𝐻𝑑𝐴𝑑𝑡_𝐻𝑑𝐴𝑑𝑡_𝐻𝑇𝐴𝑑𝑡⎤⎥⎥⎥⎥⎥⎦≤[𝐹−𝑉]⎡⎢⎢⎣𝐼^𝐻𝐴_𝐻𝐴_𝐻𝑇𝐴⎤⎥⎥⎦.(14) Using the fact that the eigenvalues of the matrix 𝐹−𝑉 all have negative real parts, it follows that the linearized differential inequality system (14) is stable whenever ℛ_𝐴<1. Consequently, (𝐼_𝐻,𝐴_𝐻,𝐴_𝐻𝑇𝐴)→(0,0,0) as 𝑡→∞. Thus, by a comparison theorem [45] (𝐼_𝐻,𝐴_𝐻,𝐴_𝐻𝑇𝐴)→(0,0,0) as 𝑡→∞, and evaluating system (7) at 𝐼_H=𝐴_𝐻=𝐴_𝐻𝑇𝐴=0 gives 𝑆_𝐻→𝑆_𝐻0 for ℛ_𝐴<1. Hence, the DFE (𝒰_0𝐻) is GAS for ℛ_𝐴<1.

3.1.4. HIV-Only Equilibrium

Expressed in terms of the equilibrium value of the force of infection 𝜆^∗𝐻, this equilibrium is given by𝒰^∗1⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩𝑆^∗𝐻=Λ^𝐻𝜇_𝐻+𝜆^∗𝐻𝐼_∗𝐻,=Λ^𝐻𝜆^∗𝐻𝐴𝜇_𝐻+𝜆^∗𝐻𝜇_𝐻+𝜌_∗𝐻,=𝜌𝜆^∗𝐻Λ_𝐻𝐴_{∗𝐻𝑇𝐴}𝜇_𝐻+𝜆^∗𝐻𝜇_𝐻+𝜌𝜇_𝐻+𝛼+𝑑_𝐴,=𝛼𝜌𝜆^∗𝐻Λ_𝐻𝜇_𝐻+𝜆^∗𝐻𝜇_𝐻+𝑑_𝐴𝜇_𝐻+𝜌𝜇_𝐻+𝛼+𝑑_𝐴.(15) The local bifurcation analysis is based on the centre manifold approach [46] as described by Theorem 4.1 in [47], stated in the appendix for convenience (also see [43] for more details). To apply the said Theorem 10 in order to establish the local asymptotic stability of the endemic equilibrium, it is convenient to make the following change of variables: 𝑆_𝐻=𝑥₁, 𝐼_𝐻=𝑥₂, 𝐴_𝐻=𝑥₃, and 𝐴_𝐻𝑇𝐴=𝑥₄, so that 𝑁_𝐻=∑^4𝑛=1𝑥_𝑛. We now use the vector notation 𝑋=(𝑥₁,𝑥₂,𝑥₃,𝑥₄)^𝑇. Then, model system (7) can be written in the form 𝑑𝑋/𝑑𝑡=𝐹=(𝑓₁,𝑓₂,𝑓₃,𝑓₄)^𝑇, where𝑥^1(𝑡)=𝑓₁=Λ_𝐻−𝛽^𝐻𝑐𝑥₂+𝜂𝑥₃+𝜅𝑥₄∑^4𝑛=1𝑥_𝑛𝑥_2𝑥₁−𝜇_𝐻𝑥₁,(𝑡)=𝑓₂=𝛽^𝐻𝑐𝑥₂+𝜂𝑥₃+𝜅𝑥₄∑^4𝑛=1𝑥_𝑛𝑥_3𝑥₁−𝜇_𝐻+𝜌𝑥₂,𝑥_4(𝑡)=𝑓₃=𝜌𝑥₂−𝜇_𝐻+𝛼+𝑑_𝐴𝑥₃,(𝑡)=𝑓₄=𝛼𝑥₃−𝜇_𝐻+𝜏𝑑_𝐴𝑥₄.(16) The Jacobian matrix of system (16) at 𝒰₀ is given by𝐽=⎡⎢⎢⎢⎢⎣𝒰−𝜇_0𝐻𝐻0𝛽_𝐻−𝛽_𝐻𝑐−𝜂𝛽_𝐻𝑐−𝜅𝛽_𝐻𝑐0𝑐−𝜌𝜇_𝐻+𝜌−𝜇_𝐻𝜂𝛽_𝐻𝑐𝜅𝛽_𝐻𝑐00+𝛼𝛼+𝑑_𝐴−𝜇_𝐻0+𝜏𝑑_𝐴⎤⎥⎥⎥⎥⎦,(17)from which it can be shown that the HIV/AIDS-induced reproduction number isℛ_𝐴=𝛽^𝐻𝑐𝜅𝜌𝛼+𝜇_𝐻+𝜏𝑑_𝐴𝜂𝜌+𝛼+𝜇_𝐻+𝑑_𝐴𝜇_𝐻+𝜌𝜇_𝐻+𝑑_𝐴𝜇_𝐻+𝛼+𝑑_𝐴.(18) If 𝛽_𝐻 is taken as a bifurcation parameter and by solving for 𝛽_𝐻 when ℛ_𝐴=1, we obtain𝛽_𝐻=𝛽^∗𝐻=𝜇^𝐻+𝜌𝜇_𝐻+𝑑_𝐴𝜇_𝐻+𝛼+𝑑_𝐴𝑐𝜅𝜌𝛼+𝜇_𝐻+𝜏𝑑_𝐴𝜂𝜌+𝛼+𝜇_𝐻+𝑑_𝐴.(19) Note that the linearized system of the transformed model (16) with 𝛽_𝐻=𝛽^∗𝐻 has a simple zero eigenvalue, which allows the use of Castillo-Chavez and Song result [47] to analyze the dynamics of (16) near 𝛽_𝐻=𝛽^∗𝐻. It can be shown that the Jacobian of (16) at 𝛽_𝐻=𝛽^∗𝐻 has a right eigenvector associated with the zero eigenvalue given by 𝑢=[𝑢₁,𝑢₂,𝑢₃,𝑢₄]^𝑇, where𝑢₁=−𝛽^𝐻𝑐𝑢₂+𝜂𝑢₃+𝜅𝑢₄𝜇_𝐻𝑢,𝑢₂>0,₃=𝜌𝑢²𝛼+𝑑_𝐴+𝜇_𝐻,𝑢₄=𝛼𝑢³𝜇_𝐻+𝜏𝑑_𝐴.(20) The left eigenvector of 𝐽(𝒰_0𝐻) associated with the zero eigenvalue at 𝛽_𝐻=𝛽^∗𝐻 is given by 𝑣=[𝑣₁,𝑣₂,𝑣₃,𝑣₄]^𝑇, where𝑣₁=0,𝑣₂=𝜌𝑣³𝜇_𝐻+𝜌−𝛽^∗𝐻𝑐,𝑣³>0,𝑣₄=𝜅𝛽^𝐻𝑐𝑣₂𝜇_𝐻+𝜏𝑑_𝐴.(21)

Computation of the Bifurcation Parameters 𝑎 and 𝑏
The application of Theorem 10 (see the appendix) entails the computation of two parameters 𝑎 and 𝑏, say. After some little algebraic manipulations and rearrangements, it can be shown that 𝑎=−2𝛽^∗𝐻𝑐𝜇_𝐻𝑣₂Λ_𝐻𝑢₂+𝑢₃+𝑢₄𝑢₂+𝜂𝑢₃+𝜅𝑢₄<0.(22) Furthermore, 𝑏=𝑐𝑢₂+𝜂𝑢₃+𝜅𝑢₄𝑣₂>0.(23) This sign of 𝑏 may be expected in general for epidemic models because, in essence, using β as a bifurcation parameter often ensures 𝑏>0 [43]. Since 𝑎<0 (which excludes any possibility of multiple equilibria and hence backward bifurcation), model system (16) has a forward (or transcritical) bifurcation at ℛ_𝐴=1, and consequently, the local stability implies global stability. This result is summarized below.

Theorem 2. The endemic equilibrium 𝒰^∗1 is locally asymptotically stable for ℛ_𝐴>1.

3.2. Schistosomiasis-Only Model

In the absence of HIV/AIDS in the community (obtained by setting HIV/AIDS-related parameters to zero from system (5)) schistosomiasis-only model is given bySchistosomiasis-onlymodel⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩𝑑𝑆^𝐻𝑑𝑡=Λ^𝐻𝑑𝐼_𝐵+𝜔𝐼_𝐵−𝜆_𝑃+𝜇_𝐻𝑆_𝐻,𝑑𝑡=𝜆^𝐻𝑑𝑀𝐼_𝐵−𝜇_𝐻+𝜔+𝑑_𝐵𝐼_𝐵,𝑑𝑡=𝑁^𝐸𝑑𝑆_𝑆𝛾𝐼_𝐵−𝜇_𝑀𝑀,𝑑𝑡=Λ^𝑆𝑑𝐼_𝑆−𝜆_𝑀𝑆_𝑆−𝜇_𝑆𝑆_𝑆,𝑑𝑡=𝜆^𝑀𝑑𝑃𝑆_𝑆−𝜇_𝑆+𝑑_𝑆𝐼_𝑆,𝑑𝑡=𝜃𝐼^𝑆with,𝜆_𝑃−𝜇_𝑃𝑃,=𝛽^𝑃𝑃(𝑡)𝜆𝑃₀+𝜖𝑃(𝑡)_𝑀,=𝛽^𝑀𝑀(𝑡)𝑀₀+𝜖𝑀(𝑡).(24) For system (24), it can be shown that the regionΦ_𝐵=⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩𝑆^𝐻,𝐼_𝐵∈ℝ²⁺∶𝑁_𝐻≤Λ^𝐻𝜇_𝐻𝑀∈ℝ₊,∶𝑀≤𝛾Λ^𝐻𝑁_𝐸(1+𝜎)𝜇_𝑀𝜇_𝐻𝑆_𝑆,,𝐼_𝑆∈ℝ²⁺∶𝑁_𝑆≤Λ^𝑆𝜇_𝑆𝑃∈ℝ₊,∶𝑃≤𝜃Λ^𝑆𝜇_𝑃𝜇_𝑆(25) is invariant and attracting. Thus, the dynamics of schistosomiasis-only model will be considered in Φ_𝐵.

3.2.1. Disease-Free Equilibrium and Stability Analysis

Model system (24) has an evident disease-free given by𝒰_0𝐵=𝑆^0𝐻,𝐼^0𝐵,𝑀⁰,𝑆^0𝑆,𝐼^0𝑆,𝑃⁰=Λ^𝐻𝜇_𝐻,0,0,Λ^𝑆𝜇_𝑆,0,0.(26) Following van den Driessche and Watmough [43], the reproduction number of the model system (24) is given byℛ_𝐵=⎷𝛽^𝑝𝑁_𝐸𝛾Λ_𝐻𝜇_𝑀𝜇_𝐻𝑃₀𝜇_𝐻+𝜔+𝑑_𝐵𝛽^𝑀𝜃Λ_𝑆=√𝜇_𝑃𝜇_𝑆𝑀₀𝜇_𝑆+𝑑_𝑆ℛ_𝐻ℛ_𝑆(27) where ℛ_𝐻=𝛽_𝑝𝑁_𝐸𝛾Λ_𝐻/𝜇_𝑀𝜇_𝐻𝑃₀(𝜇_𝐻+𝜔+𝑑_𝐵) represents the snail-man initial disease transmission and ℛ_𝑆=𝛽_𝑀𝜃Λ_𝑆/𝜇_𝑃𝜇_𝑆𝑀₀(𝜇_𝑆+𝑑_𝑆) is the man-snail initial disease transmission.

The threshold quantity ℛ_𝐵 measures the average number of new secondary cases generated by a single individual in a population where there is schistosomiasis treatment. An associated epidemiological threshold, ℛ_0𝐵, obtained using a similar technique of the next generation by considering model system (24) in the absence of schistosomiasis treatment is given byℛ_0𝐵=⎷𝛽^𝑝𝑁_𝐸𝛾Λ_𝐻𝜇_𝑀𝜇_𝐻𝑃₀𝜇_𝐻+𝑑_𝐵𝛽^𝑀𝜃Λ_𝑆=𝜇_𝑃𝜇_S𝑀₀𝜇_𝑆+𝑑_𝑆ℛ_0𝐻ℛ_0𝑆,(28) where ℛ_0𝐻=𝛽_𝑝𝑁_𝐸𝛾Λ_𝐻/𝜇_𝑀𝜇_𝐻𝑃₀(𝜇_𝐻+𝑑_𝐵) represents the snail-man initial disease transmission and ℛ_0𝑆=𝛽_𝑀𝜃Λ_𝑆/𝜇_𝑃𝜇_𝑆𝑀₀(𝜇_𝑆+𝑑_𝑆) is the man-snail initial disease transmission. Using Theorem 2 in [43], the following result is established.

Theorem 3. The disease-free equilibrium 𝒰_0𝐵 is locally asymptotically stable whenever ℛ_𝐵<1 and unstable otherwise.

Impact of Schistosomiasis Treatment in the Community
Here, the reproductive number ℛ_𝐵 is analyzed to determine whether or not treatment of schistosomiasis patients (modeled by the rate 𝜔) can lead to the effective control of schistosomiasis in the community. It follows from (27) that the elasticity [48] of ℛ_𝐵 with respect to 𝜔 can be computed using the approach in [49] as follows: 𝜔ℛ_𝐵𝜕ℛ_𝐵𝜕𝜔=−2𝜇_𝐻+𝜔𝜔+𝑑_𝐵<0.(29) The sensitivity index of the reproduction number is used to assess the impact on the relevant parameters to disease transmission. That is, the elasticity measures the effect a change in 𝜔, say, has as a proportional change in ℛ_𝐵, and from (29), we note that an increase in 𝜔 will lead to a decrease in ℛ_𝐵, thus (29) suggests that an increase in treatment of schistosomiasis patients does have a positive impact in controlling schistosomiasis in the community (assuming full compliance to the therapy, no treatment failure, and no development of resistance).

3.3. Global Stability of the Disease-Free Equilibrium

We shall use the following theorem of Castillo-Chavez et al. [50] in the sequel.

Theorem 4 (see [50]). If system (5) can be written in the form 𝑑𝑋𝑑𝑍𝑑𝑡=𝐹(𝐱,𝑍),𝑑𝑡=𝐺(𝑋,𝑍),𝐺(𝐱,0)=0,(30) where 𝑋∈ℝ^𝑚 denotes (its components) the number of uninfected individuals, 𝑍∈ℝ^𝑛 denotes (its components) the number of infected individuals including latent and infectious, and 𝐔₀=(𝐱^∗,0) denotes the disease-free equilibrium of the system. Assume that (i) for 𝑑𝑋/𝑑𝑡=𝐹(𝑋,0), 𝑋^∗ is globally asymptotically stable, (ii) 𝐺(𝑋,𝑍)=𝐴𝑍−𝐺(𝑋,𝑍), 𝐺(𝑋,𝑍)≥0 for (𝑋,𝑍)∈𝒟, where 𝐴=𝐷_𝑍𝐺(𝑋^∗,0) is an 𝑀-matrix (the off-diagonal elements of 𝐴 are nonnegative) and 𝒟 is the region where the model makes biological sense. Then the fixed point 𝐔₀=(𝐱^∗,0) is a globally asymptotic stable equilibrium of model system (5) provided ℛ_𝐵<1.

Applying Theorem 4 to model system (5) yields𝐺(𝑋,𝑌)=⎡⎢⎢⎢⎢⎢⎣𝐺¹𝐺(𝑋,𝑌)₂𝐺(𝑋,𝑌)₃𝐺(𝑋,𝑌)₄(𝑋,𝑌)⎤⎥⎥⎥⎥⎥⎦=⎡⎢⎢⎢⎢⎢⎢⎣𝛽^𝑃𝑃Λ^𝐻𝑃₀𝜇_𝐻−𝑆^𝐻𝛽𝑃₀+𝜖𝑃_𝑀𝑀Λ^𝑆𝑀₀𝜇_𝑆−𝑆^𝑆00𝑀₀+𝜖𝑀⎤⎥⎥⎥⎥⎥⎥⎦.(31) Since 𝑆^0𝐻(=Λ_𝐻/𝜇_𝐻)(1/𝑃₀)≥𝑆_𝐻/(𝑃₀+𝜖𝑃) and 𝑆_𝑆(=Λ_𝑆/𝜇_𝑆)(1/𝑀₀)≥𝑆_𝑆/(𝑀₀+𝜖𝑀), it follows that 𝐺(𝑋,𝑌)≥0. We summarise the result in Theorem 5.

Theorem 5. The disease-free equilibrium (𝒰_0𝐵) of model system (24) is globally asymptotically stable (GAS) if ℛ_𝐵<1 and unstable if ℛ_𝐵>1.

3.3.1. Schistosomiasis-Only Equilibrium

Model system (24) has an endemic equilibrium denoted by 𝒰^∗2, where𝒰^∗2⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩𝑆^𝐻∗∗=Λ^𝐻𝜇_𝐻+𝜔+𝜆^𝑃∗∗𝐼_𝐵∗∗,=Λ^𝐻𝜆^𝑃∗∗𝑀_∗∗𝜇_𝐻+𝜔+𝜆^𝑃∗∗𝜇_𝐻+𝜔+𝑑_𝐵,=𝑁^𝐸𝛾Λ_𝐻𝜆^𝑃∗∗𝑆_𝑆∗∗𝜇_𝑀𝜇_𝐻+𝜔+𝜆^𝑃∗∗𝜇_𝐻+𝜔+𝑑_𝐵,=Λ^𝑆𝜇_𝑆+𝜆^𝑆∗∗𝐼_𝑆∗∗,=Λ^𝑆𝜆^𝑀∗∗𝑃_∗∗𝜇_𝑆+𝜆^𝑀∗∗𝜇_𝑆+𝑑_𝑆,=𝜃Λ^𝑆𝜆^𝑆∗∗with𝜆𝜇_{𝑃𝑃∗∗}𝜇_𝑆+𝜆^𝑆∗∗𝜇_𝑆+𝑑_𝑆,=𝛽^𝑃𝑃^∗∗𝑃₀+𝜖𝑃^∗∗,𝜆^𝑀∗∗=𝛽^𝑀𝑀^∗∗𝑀₀+𝜖𝑀^∗∗.(32) The local asymptotic stability of the endemic equilibrium 𝒰^∗2 can also be analyzed using the centre manifold theory. In this case, the Jacobian matrix of the system at 𝒰_0𝐵 is given by𝐽=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝒰−_0𝐵𝜇_𝐻0000−𝛽^𝑃Λ_𝐻𝑃₀𝜇_𝐻0−𝜇_𝐻+𝜔+𝑑_𝐵000𝛽^𝑃Λ_𝐻𝑃₀𝜇_𝐻00𝑁0_𝐸𝛾−−𝛽_𝑀𝜇_𝑀000Λ_𝑆𝑀₀𝜇_𝑆00𝛽_𝑀−𝜇_𝑆00Λ_𝑆𝑀₀𝜇_𝑆00000−𝜇_𝑆𝜃+𝑑_𝑆−0𝜇_𝑃⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦.(33) If 𝛽_𝑃 is taken as a bifurcation parameter, and solving for 𝛽_𝑃 when ℛ_𝐵=1, we obtain𝛽_𝑃=𝛽^∗𝑃=𝜇^𝑀𝜇_𝐻𝑃₀𝜇_𝐻+𝜔+𝑑_𝐵𝑁_𝐸𝛾Λ_𝐻ℛ_𝑆.(34) The linearized system of the the model with 𝛽_𝑃=𝛽^∗𝑃 has a simple zero eigenvalue. Therefore, it can be shown that the above Jacobian has a right eigenvector given by 𝑤=[𝑤₁,𝑤₂,𝑤₃,𝑤₄,𝑤₅,𝑤₆]^𝑇, where𝑤₁=−𝛽^∗𝑃Λ_𝐻ℛ_0𝑆𝑤₃𝑃₀𝜇^2𝐻,𝑤₂=𝜇^𝑃𝛽^∗𝑃Λ_𝐻𝑤₃𝜃𝜇_𝐻+𝜔+𝑑_𝐵,𝑤³𝑤₄=𝑤₃,=−𝛽^𝑀Λ_𝑆𝑤₃𝑀₀𝜇^2𝑆,𝑤₅=𝛽^𝑀Λ_𝑆𝑤₃𝜇_𝑆+𝑑_𝑆𝑀₀𝜇_𝑆,𝑤₆=ℛ_0𝑆𝑤₃.(35) The left eigenvector of 𝐽(𝒰_0𝐵) associated with the zero eigenvalue at 𝛽_𝑃=𝛽^∗𝑃 is given by 𝑧=[𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝑧₅,𝑧₆]^𝑇, where𝑧₁=0=𝑧₄,𝑧₃>0,𝑧₂=𝑁^𝐸𝛾𝑧₃𝜇_𝐻+𝜔+𝑑_𝐵𝑧₅,=𝜇^𝑀𝑀₀𝜇_𝑆𝑧₃𝛽_𝑀Λ_𝑆,𝑧₆=𝛽^∗𝑃Λ_𝐻𝑁_𝐸𝛾𝑧₃𝑃₀𝜇_𝐻𝜇_𝑃𝜇_𝐻+𝜔+𝑑_𝐵.(36) Computation of the bifurcation coefficients 𝑎 and 𝑏 yields𝑎=−2𝑧₃𝑤²³𝑁^𝐸𝛾ℛ^2𝑆𝛽^∗𝑃Λ_𝐻𝛽^∗𝑃+𝜖𝜇_𝐻𝜇_𝐻+𝜔+𝑑_𝐵𝑃²⁰𝜇^2𝐻+𝑁^𝐸𝛾𝜃𝛽^∗𝑃Λ_𝐻𝛽_𝑀Λ_𝑆𝛽_𝑀+𝜖𝜇_𝑆𝑏=𝑁_𝐸𝜇_𝐻+𝜔+𝑑_𝐵𝜇_𝑆+𝑑_𝑆𝑃₀𝜇_𝐻𝑀²⁰𝜇^2𝑆<0,𝛾ℛ_𝑆Λ_𝑆𝑧₃𝑤₃𝜇_𝐻+𝜔+𝑑_𝐵𝑃₀𝜇_𝑆>0.(37) Thus, the following result is established.

Theorem 6. The unique endemic equilibrium 𝒰^∗2 is locally asymptotically stable for ℛ_𝐵>1.

Since 𝑎<0, local stability of 𝒰^∗2 implies its global stability.

4. HIV/AIDS and Schistosomiasis Model

Model system (5) has evident disease-free (DFE) given by𝒰₀==𝑆Λ^0𝐻_𝐻,𝐼^0𝐵,𝐼^0𝐻,𝐴^0𝐻,𝐴^0𝐻_𝑇𝐴,𝐼^0𝐻_𝐵,𝐴^0𝐻_𝐵,𝐴^0𝐻_𝑇𝐵,𝑀⁰,𝑆^0𝑆,𝐼^0𝑆,𝑃⁰𝜇_𝐻,0,0,0,0,0,0,0,0,Λ^𝑆𝜇_𝑆,0,0.(38) Following van den Driessche and Watmough [43], the reproduction number of the model isℛ_𝐻𝐵=maxℛ_𝐴,ℛ_𝐵(39) with ℛ_𝐴 and ℛ_𝐵 defined as earlier in Section 3 above. Using Theorem 2 in [43], the following result is established.

Theorem 7. The disease-free equilibrium 𝒰₀ is locally asymptotically stable whenever ℛ_𝐻𝐵<1 and unstable otherwise.

4.1. Sensitivity Analysis

In this section we investigate the effects of HIV/AIDS on schistosomiasis and vice versa, in the presence and absence of the aforementioned intervention strategies.

Impact of Schistosomiasis on HIV/AIDS in the Absence of Control Measures
To analyze the effects of schistosomiasis on HIV/AIDS and vice versa in the absence of control measures for either HIV/AIDS or schistosomiasis, we begin by introducing the following notation; in the absence of antiretroviral therapy (𝛼=0) the reproductive number is denoted by ℛ_0𝐴 and also in the absence of schistosomiasis treatment (𝜔=0), ℛ_𝐵 = ℛ_0𝐵. Thus, to express ℛ_0𝐵 in terms of ℛ_0𝐴, we solve for 𝜇_𝐻 and obtain 𝜇_𝐻=−𝜙¹ℛ_0𝐴+𝜙₂+𝜙₃ℛ²⁰_𝐴+𝜙₄ℛ_0𝐴+𝜙₅2ℛ_0𝐴,(40) where 𝜙𝜙₁₂=𝜌+𝑑_𝐴,𝜙=−𝛽_𝐻𝑐,₃𝜙₄=𝜌−𝑑_𝐴²,𝜙=2𝛽_𝐻𝑐𝑑_𝐴+𝜌(2𝜂−1)₅,=𝛽_𝐻𝑐².(41) Let 𝜙₃ℛ²⁰_𝐴+𝜙₄ℛ_0𝐴+𝜙₅=𝜙₆ℛ_0𝐴+𝜙₇, then, (40) becomes 𝜇_𝐻=𝜙⁶−𝜙₁ℛ_0𝐴+𝜙₇−𝜙₂2ℛ_0𝐴.(42) Substituting (42) into the expression for ℛ_0𝐵, we haveℛ²⁰_𝐵=4ℛ^0𝑆ℛ²⁰_𝐴𝛽_𝑃𝑁_𝐸𝛾Λ_𝐻𝜇_𝑀𝑃₀𝜙⁶−𝜙₁ℛ_0𝐴+𝜙₇−𝜙₂²+2𝑑_𝐴ℛ_0𝐴𝜙₆−𝜙₁ℛ_0𝐴+𝜙₇−𝜙₂.(43) Differentiating ℛ_0𝐵 partially with respect to ℛ_0𝐴 yields𝜕ℛ_0𝐵𝜕ℛ_0𝐴=4ℛ^0𝑆ℛ_0𝐴𝛽_𝑃𝑁_𝐸𝛾Λ_𝐻ℛ⁰_𝐴𝜙₇−𝜙₂𝜙₆−𝜙₁−𝑑_𝐴+𝜙₇−𝜙₂²𝜇_𝑀𝑃₀ℛ_0𝐵𝜙⁶−𝜙₁ℛ_0𝐴+𝜙₇−𝜙₂²+2𝑑_𝐴ℛ_0𝐴𝜙₆−𝜙₁ℛ_0𝐴+𝜙₇−𝜙₂².(44) Now, whenever (44) is greater than zero, an increase in HIV/AIDS cases results in an increase of schistosomiasis cases in the community. If (44) is equal to zero, this implies that HIV/AIDS cases have no effect on the transmission dynamics of schistosomiasis. Setting ℛ_0𝐵=1 and expressing 𝜇_𝐻 as the subject of formula, we have 𝜇_𝐻=−𝑑^𝐵𝜃₁ℛ_0𝐵+𝜃₁𝑑_𝐵ℛ_0𝐵²+4𝜃₂2𝜃₁ℛ_0𝐵,(45) where 𝜃₁=𝜇_𝑀𝑃₀and𝜃₂=𝜇_𝑀𝑃₀𝛽_𝑃𝑁_𝐸𝛾Λ_𝐻ℛ_0𝑆. Consider (𝜃₁𝑑_𝐵ℛ_0𝐵)²+4𝜃₂=𝜃₃ℛ_0𝐵+𝜃₄ such that (𝜃₃−𝑑_𝐵𝜃₁)ℛ_0𝐵+𝜃₄>0. Then, ℛ_0𝐴 expressed in terms of ℛ_0𝐵 reads ℛ_0𝐴=2𝜃¹𝛽_𝐻𝑐𝜅ℛ²⁰_𝐵+𝜃₄ℛ_0𝐵ℎ₁ℛ²⁰_𝐵+ℎ₂ℛ_0𝐵+ℎ₃,(46) where ℎℎ₁=𝜃₃−𝑑_𝐵𝜃₁²+4𝜌𝑑_𝐴𝜃²¹+2𝜃₁𝜃₃−𝑑_𝐵𝜃₁𝜌+𝑑_𝐴>0,₂ℎ=2𝜃₄𝜃₃−𝑑_𝐵𝜃₁+2𝜃₁𝜃₄𝜌+𝑑_𝐴>0,₃𝜅=𝜃²⁴>0,₁=𝜃₃+𝜃₁2𝜂𝜌+2𝑑_𝐴−𝑑_𝐵>0.(47) Partially differentiating ℛ_0𝐴 with respect to ℛ_0𝐵 yields 𝜕ℛ_0𝐴𝜕ℛ_0𝐵=𝜅¹ℎ₂−𝜃₄ℎ₁ℛ²⁰_𝐵+2𝜅₁ℎ₃ℛ_0𝐵+𝜃₄ℎ₃ℎ¹ℛ²⁰_𝐵+ℎ₂ℛ_0𝐵+ℎ₃².(48) Thus, whenever 𝜅₁ℎ₂≥𝜃₄ℎ₁, (48) is strictly positive meaning that schistosomiasis enhances HIV infection as a damaged urethra has increased chances of HIV entering the blood stream. The relationship between the HIV/AIDS basic reproduction number and the schistosomiasis basic reproduction number is illustrated graphically in Figure 2 using parameter values from Table 1.

Figure 2 

Relationship between the HIV/AIDS and the schistosomiasis basic reproduction numbers.

The graph in Figure 2 shows that an increase in the schistosomiasis-induced basic reproduction number results in an increase of the HIV/AIDS-induced basic reproduction number, suggesting that infection by schistosomiasis enhances the chances of HIV infection per sexual contact. This is as a result of the eggs of the parasites causing injury in the reproductive organs which enhance the transmission of sexually transmitted diseases such as HIV/AIDS and Gonorrhoea [51]. Thus, schistosomiasis control has a positive impact in controlling the transmission dynamics of HIV/AIDS.

Impact of Schistosomiasis Treatment on HIV/AIDS
Expressing ℛ_0𝐵 in terms of ℛ_𝐵, we obtain ℛ_0𝐵=𝜇^𝐻+𝜔+𝑑_𝐵ℛ_𝐵𝜇_𝐻+𝑑_𝐵.(49) Substituting (49) into (46) yieldsℛ_0𝐴=2𝛽^𝐻𝑐𝜃₁ℛ_𝐵𝜇_𝐻+𝜔+𝑑_𝐵𝜃₄𝜇_𝐻+𝑑_𝐵+𝜇_𝐻+𝜔+𝑑_𝐵𝜅₁ℛ_𝐵𝜇_𝐻+𝜔+𝑑_𝐵²ℎ₁ℛ^2𝐵+𝜇_𝐻+𝑑_𝐵𝜇_𝐻+𝜔+𝑑_𝐵ℎ₂ℛ_𝐵+𝜇_𝐻+𝑑_𝐵²ℎ₃.(50)Partially differentiating ℛ_0𝐴 with respect to 𝜔, we have 𝜕ℛ_0𝐴𝜕𝜔=−2𝛽^𝐻𝑐𝜃₁𝑘₃𝑘₄[Θ−1],(51) where Θ=𝑘₁𝑘₂/𝑘₃𝑘₄, with 𝑘𝑘₁₂=ℛ_𝐵𝜇_𝐻+𝑑_𝐵ℎ₂+2𝜁ℎ₁ℛ_𝐵,𝑘₃=𝜁ℛ_𝐵𝜇_𝐻+𝑑_𝐵𝜃₄+𝜁𝜅₁ℛ_𝐵,𝑘₄=ℛ_𝐵𝜇_𝐻+𝑑_𝐵𝜃₄+2𝜁𝜅₁ℛ_𝐵,𝜁==𝜁𝜇ℛ_𝐵𝐻𝜁ℎ₁ℛ_𝐵+𝜇_𝐻+𝑑_𝐵ℎ₂+ℎ₃,+𝜔+𝑑_𝐵.(52) Since ℛ_0𝐴 is a decreasing function of 𝜔, schistosomiasis treatment will have a positive impact on the dynamics of HIV/AIDS if Θ>1, no impact if Θ=1, and a negative impact if Θ<1. We summarize the result in lemma 5.

Lemma 5. Schistosomiasis (bilharzia) treatment for model system (5) only, will have (i)a positive impact on schistosomiasis and HIV/AIDS coinfection control if Θ>1,(ii)no impact on schistosomiasis and HIV/AIDS coinfection control if Θ=1,(iii)a negative impact on schistosomiasis and HIV/AIDS coinfection control if Θ<1.

The synergy between HIV and other diseases such as schistosomiasis provides more opportunities to combat HIV/AIDS by treating its coinfections with these other diseases.

4.2. Global Stability of the Disease-Free Equilibrium (𝒰₀)

We shall use the following theorem of Castillo-Chavez et al. [50] in the sequel.

Theorem 8 (see [50]). If system (5) can be written in the form 𝑑𝑋𝑑𝑍𝑑𝑡=𝐹(𝐱,𝑍),𝑑𝑡=𝐺(𝑋,𝑍),𝐺(𝐱,0)=0,(53) where 𝑋∈ℝ^𝑚 denotes (its components) the number of uninfected individuals, 𝑍∈ℝ^𝑛 denotes (its components) the number of infected individuals including latent, infectious, and so forth, 𝐔₀=(𝐱^∗,0) denotes the disease-free equilibrium of the system. Assume that (i) for 𝑑𝑋/𝑑𝑡=𝐹(𝑋,0),𝑋^∗ is globally asymptotically stable, (ii) 𝐺(𝑋,𝑍)=𝐴𝑍−𝐺(𝑋,𝑍), 𝐺(𝑋,𝑍)≥0 for (𝑋,𝑍)∈𝒟, where 𝐴=𝐷_𝑍𝐺(𝑋^∗,0) is an 𝑀-matrix (the off-diagonal elements of 𝐴 are nonnegative) and 𝒟 is the region where the model makes biological sense. Then the fixed point 𝐔₀=(𝐱^∗,0) is a globally asymptotic stable equilibrium of model system (5) provided ℛ_𝐻𝐵<1.

Applying Theorem 8 to model system (5) yields𝐺=(⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑋𝐺,𝑌)₁𝐺(𝑋,𝑌)₂𝐺(𝑋,𝑌)₃𝐺(𝑋,𝑌)₄𝐺(𝑋,𝑌)₅𝐺(𝑋,𝑌)₆𝐺(𝑋,𝑌)₇𝐺(𝑋,𝑌)₈𝐺(𝑋,𝑌)₉𝐺(𝑋,𝑌)₁₀(𝑋,𝑌)⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝛿𝜆^𝐻𝐼_𝐵+𝛽_𝑃𝑃𝑃(₀𝑡)𝜇Λ^𝐻_𝐻−𝑃(𝑡)𝑆^𝐻(𝑡)𝜆𝑃₀+𝜖𝑃(𝑡)_𝑃𝐼_𝐻+𝑁_𝐻𝜆_𝑃1−𝑁𝑆^𝐻_𝐻𝐴_𝐻𝜆_𝑃𝐴_𝐻𝑇𝐴−𝜆_𝑃𝐼_𝐻−𝛿𝜆_𝐻𝐼_𝐵−𝜆_𝑃𝐴_𝐻−𝜆_𝑃𝐴_𝐻𝑇𝐴𝛽_𝑀0𝑀𝑀₀Λ𝜇^𝑆_𝑆−𝑀𝑆^𝑆0𝑀₀+𝜖𝑀⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦.(54) The fact that 𝐺₄(𝑋,𝑌)<0, 𝐺₅(𝑋,𝑌)<0, and 𝐺₆(𝑋,𝑌)<0 implies that 𝐺(𝑋,𝑌) may not be greater or equal to zero. Consequently, 𝒰₀ may not be globally asymptotically stable for ℛ_𝐻𝐵<1. This suggests the possible existence of multiple equilibria.

4.3. Endemic Equilibria and Its Stability

For model system (5), there are three possible endemic equilibria: the case where there is HIV only, the case where there is schistosomiasis only (which have been discussed in Section 3), and the case when both schistosomiasis and HIV coexist.

4.3.1. Interior Endemic Equilibrium

This occurs when both infections coexist in the community. The interior equilibrium is given by𝒰^∗3=𝑆𝐴^{𝐻∗∗∗}_{∗∗∗𝐻𝑇𝐵},𝐼^{𝐵∗∗∗},𝐼^{𝐻∗∗∗},𝐴^{𝐻∗∗∗},𝐴^{∗∗∗𝐻𝑇}_𝐴,𝐼^{𝐻∗∗∗}_𝐵,𝐴^{𝐻∗∗∗}_𝐵,𝑀^∗∗∗,𝑆^{𝑆∗∗∗},𝐼^{𝑆∗∗∗},𝑃^∗∗∗.(55)The local asymptotic stability of this endemic equilibrium can be analyzed using the centre manifold theory similar to the analysis of 𝒰^∗1and𝒰^∗2, but it is not done here to avoid repetition. Thus, we claim the following result for the stability of 𝒰^∗1and𝒰^∗2.

Theorem 9. If ℛ_𝐻𝐵>1 with ℛ_𝐵>1 and ℛ_𝐴>1, then, the endemic equilibrium point 𝒰₃ is locally asymptotically stable whenever ℛ_𝐻𝐵>1.

5. Numerical Simulations

In order to illustrate the results of the foregoing analysis, numerical simulations of the full HIV-schistosomiasis model are carried out, using parameter values given in Table 1. The scarcity of data on HIV schistosomiasis codynamics limits our ability to calibrate, but, for the purpose of illustration, other parameter values are assumed. These parsimonious assumptions reflect the lack of information currently available on the coinfection of the two diseases.

Figure 3 depicts the effects of schistosomiasis on the dynamics of HIV in the community. The time series plots in Figure 3 suggest that the presence of schistosomiasis in the community might increase the prevalence of HIV/AIDS. These numerical results are in agreement with our analytical results. We note that 𝐼_𝐻 and 𝐴_𝐻 are not reflecting the disease-free equilibrium, and the convergence is simply due to scale.

(a)

(b)

Figure 3 

Numerical results of model system (5) showing time series plots of infectives either singly infected with HIV or dually infected with HIV and schistosomiasis for both cases (i.e., either displaying clinical symptoms of AIDS or not), using various initial conditions and parameter values from Table 1.

6. Summary and Conclusion

While schistosomiasis is the second most prevalent neglected tropical disease after hookworm infection (192 million cases worldwide) [5], HIV on the other hand which has killed more than 25 million people since first recognized in 1981 currently affects 33.4 million people, with deaths due to HIV/AIDS-related illnesses standing at about 2 million in 2008 [6]. A mathematical model for investigating the coinfection of schistosomiasis and HIV/AIDS is derived. Comprehensive and qualitative mathematical techniques were used to analyze steady states of the model. The disease-free equilibrium is shown to be locally asymptotically stable when the associated epidemic threshold known as the basic reproduction number for the model is less than unity. Center manifold theory is used to show that the schistosomiasis-only and HIV/AIDS-only endemic equilibria are locally asymptotically stable when the associated reproduction numbers are greater than unity. The impact of schistosomiasis and its treatment on the dynamics of HIV/AIDS is also investigated. Numerical results are provided to illustrate some of analytical results.

In this study, the impact of schistosomiasis and its treatment on the transmission dynamics of HIV/AIDS in the community is investigated by formulating a mathematical model that incorporates both key epidemiological parameters of both schistosomiasis and HIV/AIDS. Mathematical and numerical analysis of the model suggests that schistosomiasis may increase the prevalence of HIV/AIDS in the community. Analysis of the impact of schistosomiasis treatment has shown that the impact of this form of treatment depends on the sign of a certain threshold parameter Θ, and for Θ>1, schistosomiasis treatment will have a positive impact, for Θ=1, no impact, and for Θ<1, a negative impact on controlling the co-interaction of the two diseases. We, however, note that from schistosomiasis and HIV/AIDS epidemiology, realistic parameter values always yield 1<Θ. Consequently, schistosomiasis treatment will always have a positive impact on the control of both schistosomiasis and HIV/AIDS codynamics. Thus, schistosomiasis treatment can reduce the burden of schistosomiasis and HIV/AIDS coinfection in areas of extreme poverty, especially among the rural poor and some disadvantaged urban populations since it is less expensive and usually available in government clinics and hospitals. This outcome highlights the fact that global public health challenges require comprehensive and multipronged approaches to dealing with them. Current efforts that focus on a single infection at a time may be losing substantial rewards of dealing synergistically and concurrently with multiple infectious diseases [7].

Appendix

In order to establish the conditions for the existence of a bifurcation, we use Theorem 10 proven in [47].

Theorem 10. Consider the following general system of ordinary differential equations with a parameter 𝜙: 𝑑𝑥𝑑𝑡=𝑓(𝑥,𝜙),𝑓∶ℝ^𝑛×ℝ⟶ℝ,𝑓∈ℂ²(ℝ^𝑛×ℝ),(A.1) where 0 is an equilibrium of the system that is 𝑓(0,𝜙)=0 for all 𝜙, and assume that (𝐴1)𝐴=𝐷_𝑥𝑓(0,0)=((𝜕𝑓_𝑖/𝜕𝑥_𝑗)(0,0)) is linearization of system (A.1) around the equilibrium 0 with 𝜙 evaluated at 0. Zero is a simple eigenvalue of 𝐴, and other eigenvalues of 𝐴 have negative real parts, (𝐴2) matrix 𝐴 has a right eigenvector 𝑢 and a left eigenvector 𝑣 corresponding to the zero eigenvalue. Let 𝑓_𝑘 be the 𝐾th component of 𝑓 and 𝑎=^𝑛_𝑘_,𝑖,𝑗=1𝑣_𝑘𝑢_𝑖𝑢_𝑗𝜕²𝑓_𝑘𝜕𝑥_𝑖𝜕𝑥_𝑗𝑏=_𝑛(0,0),_𝑘_,𝑖=1𝑣_𝑘𝑢_𝑖𝜕²𝑓_𝑘𝜕𝑥_𝑖𝜕𝜙(0,0).(A.2) The local dynamics of (A.1) around 0 are totally governed by 𝑎 and 𝑏. (i)𝑎>0, 𝑏>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)𝑎<0, 𝑏<0. When 𝜙<0 with |𝜙|≪1, 0 is unstable and when 0<𝜙≪1, asymptotically stable, and there exists a positive unstable equilibrium. (iii)𝑎>0, 𝑏<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)𝑎<0, 𝑏>0. When 𝜙 changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly, a negative equilibrium becomes positive and locally asymptotically stable.

Computations of 𝑎 and 𝑏
For system (16), the associated nonzero partial derivatives of 𝐹 associated with 𝑎 at the disease-free equilibrium is given by 𝜕²𝑓₂𝜕𝑥₂𝜕𝑥₃=𝜕²𝑓₂𝜕𝑥₃𝜕𝑥₂=−𝛽^∗𝐻𝑐(1+𝜂)𝜇_𝐻Λ_𝐻𝜕₂,𝑓₂𝜕𝑥²²=−2𝛽^∗𝐻𝑐𝜇_𝐻Λ_𝐻𝜕₂,𝑓₂𝜕𝑥₂𝜕𝑥₄=𝜕²𝑓₂𝜕𝑥₄𝜕𝑥₂=−𝛽^∗𝐻𝑐(1+𝜅)𝜇_𝐻Λ_𝐻𝜕₂,𝑓₂𝜕𝑥²³=−2𝛽^∗𝐻𝑐𝜂𝜇_𝐻Λ_𝐻𝜕₂,𝑓₂𝜕𝑥₃𝜕𝑥₄=𝜕²𝑓₂𝜕𝑥₄𝜕𝑥₃=−𝛽^∗𝐻𝑐(𝜂+𝜅)𝜇_𝐻Λ_𝐻𝜕₂,𝑓₂𝜕𝑥²⁴=−2𝛽^∗𝐻𝑐𝜅𝜇_𝐻Λ_𝐻.(A.3) From (A.3), it follows that 𝑎=−2𝛽^∗𝐻𝑐𝜇_𝐻𝑣₂Λ_𝐻𝑢₂+𝑢₃+𝑢₄𝑢₂+𝜂𝑢₃+𝜅𝑢₄<0.(A.4) For the sign of 𝑏, it is associated with the following nonvanishing partial derivatives of 𝐹: 𝜕²𝑓₂𝜕𝑥₂𝜕𝛽^∗𝐻=𝑐,𝜕²𝑓₂𝜕𝑥₃𝜕𝛽^∗𝐻=𝑐𝜂,𝜕²𝑓₂𝜕𝑥₄𝜕𝛽^∗𝐻=𝑐𝜅,(A.5) from which it follows that 𝑏=𝑐𝑢₂+𝜂𝑢₃+𝜅𝑢₄𝑣₂>0.(A.6)Thus, 𝑎<0 and 𝑏>0 and from Theorem 10 item (iv), the result follows.

Acknowledgment

The authors thank the reviewers for comments and suggestions.