Local Stability and Bifurcation Analysis of Onchocerciasis Transmission Dynamics With Nonlinear Incidence Functions in Two Interacting Populations

1. Introduction

Onchocerciasis is one of the neglected tropical diseases caused by the parasite Onchocerca Volvulus, a filarial nematode [2] . The disease is transmitted from one person to another by repeated bites of black flies. The disease is endemic in Sub-saharan Africa. Many researchers have worked on many ways to reduce the spread of the disease. For instance, Remme et al. [11] used skin snip survey in West Africa to investigate the impact of controlling black flies by larviciding. Plaisier et al. [10] used micro simulation model to determine the period required for combining annual ivermectin treatment and vector control in the onchocerciasis Control Programme in West Africa. Alley et al. [1] used a computer simulation model to study prevention of onchocerciasis by using macrofilaricide which kills the adult worms. Asha Hassan & Nyimvua Shaban [3] investigated the effects of four control strategies on the spread of the disease.

In this paper, we consider onchocerciasis transmission dynamics with nonlinear incidence functions. The human population is sub-divided into four compartments and the vector population is sub-divided into three compartments. We show local asymptotic behaviour in disease-free and endemic equilibria and also the bifurcation analysis examined. The rest of the paper is organized as follows: the description of the model and theorems on positivity of solutions are given in section 2 while section 3 is devoted to the proof local stability theorems and the bifurcation analysis of the model is done in section 4.

2. Model Description

Two interacting populations are considered; the humans and the black-flies populations. The human population is partitioned into four compartments: the susceptible human compartment; $S_h$,, the exposed compartment; $E_h$, the infectious human compartment; $I_h$ and the recovered human compartment; $R_h$. The black-fly population is partitioned into three compartments: susceptible vector; $S_v$, the exposed vector compartment; $E_v$ and the infective vector compartment. The total human and vector populations at any given time, t, are respectively given by; $N=S_h\left(t\right)+E_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)$ and $V=S_v\left(t\right)+E_v\left(t\right)+I_v\left(t\right)$. We assume that the transmission of onchocerciaisis in susceptible hosts is only through contact with infectious vector. We also assume that susceptible vector becomes infectious as a result of contact with infectious hosts during blood meal. The population under study is assumed to be large enough to be modelled deterministically. The following system of non-linear ordinary differential equations, with non-negative initial conditions, describes the dynamics of onchocerciaisis epidemics.

(2.1)

subject to the following initial conditions:

$$\begin{array}{c}{S}_h(0,x_i)=S_{0h}\left(x_i\right),E_h(0,x_i)=E_{0h}\left(x_i\right),\\ I_h(0,x_i)=I_{0h}\left(x_i\right),R_h(0,x_i)=R_{0h}\left(x_i\right)S_m\left(0\right)=S_{0m},E_m\left(0\right)=E_{0m},I_m\left(0\right)=I_{0m}\end{array}$$

(2.2)

Symbols	Definitionss
$S_h(t,x_i)$	Number of susceptible humans at time t and discrete age $x_i$
$E_h(t,x_i)$	Number of exposed humans at time t and discrete age $x_i$
$I_h(t,x_i)$	Number of infectious humans at time t and discrete age $x_i$
$R_h(t,ai)$	Number of recovered humans at time t and discrete age $x_i$
$S_v\left(t\right)$	Number of susceptible black-flies at time t
$E_v\left(t\right)$	Number of exposed black-flies at time t
$I_v\left(t\right)$	Number of infectious black-flies at time t
${\mathsf{\Psi }}_h\left(x_i\right)$	Recruitment term of the susceptible humans at discrete age $x_i$
${\mathsf{\Psi }}_v$	Recruitment term of the susceptible vectors
$\delta$	Biting rate of the vector
${\lambda }_h\left(x_i\right)$	Probability that a bite by an infectious vector results in transmission of disease to human at discrete age $x_i$
${\lambda }_v$	Probability that a bite results in transmission of parasite to a susceptible vector
${\mu }_h\left(x_i\right)$	Per capita death rate of humans at discrete age $x_i$
${\mu }_v$	Per capita death rate of vector
${\gamma }_h\left(x_i\right)$	Disease-induced death rate of humans at discrete age $x_i$
${\gamma }_v$	Disease-induced death rate of vectors
${\alpha }_h\left(x_i\right)$	Per capita rate of progression of humans from the exposed state to the infectious state at discrete age $x_i$
${\alpha }_v$	Per capita rate of progression of vectors from the exposed state to the infectious state
$r\left(x_i\right)$	Per capita recovery rate for humans from the infectious state to the recovered state due to treatment at discrete age $x_i$
$\omega \left(x_i\right)$	Per capita transition rate of recovered humans to the susceptible state at discrete age $x_i$
${\nu }_h\left(x_i\right)$	Humans disease-inhibiting factor at discrete age $x_i$
${\nu }_v$	Vectors disease-inhibiting factor

Model assumptions

The formulation of the compartmental model is based on the following assumptions:

1.

That all humans are born susceptible. That is, humans are liable to contract the disease.

2.

That the susceptible humans, when infected, becomes exposed humans who are not yet infectious.

3.

That the exposed humans progress to become infectious only.

4.

That the infectious humans may either die naturally or as a result of the disease, and if not, they become recovered humans due to treatment.

5.

That the recovered humans become susceptible again.

6.

All black-flies are born susceptible.

7.

That the susceptible black-flies, when infected, becomes exposed black-flies who are not yet infectious.

8.

That the exposed black-flies progress to become infectious only.

9.

That the infectious black-flies remain infectious for life. That is, there is no recovered class for black-fly population.

2.1. Existence and Positivity of Solutions

In this section, we analyse the general properties of the system (2.1) with positive initial conditions. It describes the population dynamics both in human and black-fly populations. The system is biologically relevant in the set given by

$$\mathsf{\Omega }=(S_h(t,x_i),E_h(t,x_i),I_h(t,x_i),R_h(t,x_i))\in {\mathbb{R}}_{+}^4:N_h\le \sum _{i=0}^L\frac{{\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)},(S_v\left(t\right),E_v\left(t\right),I_v\left(t\right))\in {\mathbb{R}}_{+}^3:N_v\le \frac{{\mathsf{\Psi }}_v}{{\mu }_v}$$

Here, the following results are provided which guarantee that the model governed by system (2.1) is mathematically well-posed in a feasible region $\mathsf{\Omega }$ defined by:

$$\mathsf{\Omega }={\mathsf{\Omega }}_h\times {\mathsf{\Omega }}_v\subset {\mathbb{R}}^4\times {\mathbb{R}}^3$$

Theorem 1: There exists a domain $\mathsf{\Omega }$ in which the solution set $S_h(t,x_i),E_h(t,x_i),I_h(t,x_i),R_h(t,x_i),S_v\left(t\right),E_v\left(t\right),I_v\left(t\right)$ is contained and bounded.

Proof If the total human population size is given by $N_h=S_h(t,x_i)+E_h(t,x_i)+I_h(t,x_i)+R_h(t,x_i)$, and the total size of black-fly population is $N_v=S_v\left(t\right)+E_v\left(t\right)+I_v\left(t\right)$. From model (2.1), we have that

$$\frac{d{\mathit{N}}_h(t,x_i)}{dt}\le {\mathsf{\Psi }}_h\left(x_i\right)-\sum _{i=0}^L{\mu }_h\left(x_i\right)N_h(t,x_i)$$

(2.3)

and

$$\frac{d{\mathit{N}}_v}{dt}\le {\mathsf{\Psi }}_v-{\mu }_v{N}_v$$

(2.4)

It follows from (2.3) and (2.4) that

$$N_h(t,x_i)\le \frac{{\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}[1-e^{1-{\mu }_h\left(x_i\right)t]+N_h(0,x_i)e^{-{\mu }_h\left(x_i\right)t}]}$$

and

$$N_v\le \frac{{\mathsf{\Psi }}_v}{{\mu }_v}[1-e^{-{\mu }_{v}t}]+N_v\left(0\right)e^{-{\mu }_{v}t}$$

Taking the $limsup$ as $t\to \infty$ gives $N_h\le \frac{{\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$ and $N_v\le \frac{{\mathsf{\Psi }}_v}{{\mu }_v}$. This shows that all solutions of the humans population only are confined in the solution set ${\mathsf{\Omega }}_h$ and all solutions of the black-fly population are confined in ${\mathsf{\Omega }}_v$. It also suffices to say that $\mathsf{\Omega }$ is positively invariant as $N_h(t,x_i)\le {\sum }_{i=0}^L\frac{{\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$ whenever $N_h(0,x_i)\le \frac{{\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$ and $N_v\left(t\right)\le \frac{{\mathsf{\Psi }}_v}{{\mu }_v}$ if $N_v\left(0\right)\le \frac{{\mathsf{\Psi }}_v}{{\mu }_v}$, Therefore the solution set for the model (2.1) exists and is given by $\mathsf{\Omega }={\mathsf{\Omega }}_h\times {\mathsf{\Omega }}_v\subset {\mathbb{R}}_{+}^4\times {\mathbb{R}}_{+}^3$ □

It remains to show that the solutions of system (2.1) are nonnegative in $\mathsf{\Omega }$ for any time $t>0$ since the variables represent human and black-fly populations.

Theorem 2: The solutions, $S_h(t,x_i)$, $E_h(t,x_i)$, $I_h(t,x_i)$, $R_h(t,x_i)$, $S_v\left(t\right)$, $E_v\left(t\right)$, $I_v\left(t\right)$, of model (2.1) with nonnegative initial conditions in $\mathsf{\Omega }$, remain nonnegative in $\mathsf{\Omega }$ for all $t>0$.

Proof: Given that the initial conditions, $S_{0h}\left(x_i\right)$, $E_{0h}\left(x_i\right)$, $I_{0h}\left(x_i\right)$, $R_{0h}\left(x_i\right)$, $S_{0v}$, $E_{0v}$, $I_{0v}$, are non-negative and from (2.1),

$$\frac{d{\mathit{S}}_h(t,x_i)}{dt}+\sum _{i=0}^L\left[\frac{b{\lambda }_h\left(x_i\right)I_v\left(t\right)}{1+{\nu }_h\left(x_i\right)I_v\left(t\right)}+{\mu }_h\left(x_i\right)\right]S_h(t,x_i)\ge 0$$

so that

$$\frac{d}{dt}\left[\sum _{i=0}^L{S}_h(t,x_i)exp\left({\int }_0^t\frac{b{\lambda }_h\left(x_i\right)I_v\left(\eta \right)}{1+{\nu }_h\left(x_i\right)I_v\left(\eta \right)}d\eta +{\mu }_h\left(x_i\right)t\right)\right]\ge 0,$$

(2.5)

Integrating (2.5), we have

$$\sum _{i=0}^L{S}_h(t,x_i)\ge \sum _{i=0}^L{S}_{0h}\left(x_i\right)exp\left[-\left({\int }_0^t\frac{b{\lambda }_h\left(x_i\right)I_v\left(\eta \right)}{1+{\nu }_h\left(x_i\right)I_v\left(\eta \right)}d\eta +{\mu }_h\left(x_i\right)t\right)\right]\ge 0,$$

which implies that for all $t>0$ and for all $a\in {\mathbb{R}}_{+}$, we have

$$S_h(t,x_i)\ge \sum _{i=0}^L{S}_{0h}\left(x_i\right)exp\left[-\left({\int }_0^t\frac{b{\lambda }_h\left(x_i\right)I_v\left(\eta \right)}{1+{\nu }_h\left(x_i\right)I_v\left(\eta \right)}d\eta +{\mu }_h\left(x_i\right)t\right)\right]\ge 0.$$

Hence, $S_h(t,x_i)>0$ for any arbitrary $x_i$. Also, we have

$$\frac{d{\mathit{E}}_h(t,x_i)}{dt}+\sum _{i=0}^L\left(({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))\right)E_h(t,x_i)\ge 0$$

so that

$$\frac{d}{dt}\left[\sum _{i=0}^L{E}_h(t,\left(x_i\right))exp({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right)t)\right]\ge 0$$

(2.6)

Integrating (2.6), we have for all $t>0$ and for all $a\in mathbb{R}_{+}$, that

$$E_h(t,a_{)}\ge \sum _{i=0}^L{E}_{0h}\left(x_i\right)exp\left[-({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))t\right]$$

Hence, $E_h(t,x_i)>0$ for any arbitrary $x_i$ Also we have

$$\frac{d{\mathit{I}}_h(t,x_i)}{dt}\ge -\sum _{i=0}^L(r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))I_h\left(t\right)$$

so that

$$\frac{d}{dt}\left[I_h\left(t\right)exp(r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))t\right]\ge 0$$

(2.7)

Similarly, (2.7) becomes

$$I_h(t,a)\ge \sum _{i=0}^L{I}_{0h}exp\left[-(r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))t\right]>0\mathrm{for}\mathrm{all}t>0\mathrm{for}\mathrm{all}a\in {\mathbb{R}}_{+}$$

Hence, $I_h(t,x_i)>0$ for any arbitrary $x_i$. Also from (2.1), we have

$$\frac{d{\mathit{R}}_h(t,x_i)}{dt}+\sum _{i=0}^L({\mu }_h\left(x_i\right)+w\left(x_i\right))R_h(t,x_i)\ge 0$$

and we have

$$\frac{d}{dt}\left[\sum _{i=0}^L{R}_h(t,x_i)exp(({\mu }_h\left(x_i\right)+w\left(x_i\right))t\right]\ge 0$$

(2.8)

Integrating (2.8), we have, for all $t>0$ and $a\in \mathbb{R}$, that

$$R_h(t,a)\ge \sum _{i=0}^L{R}_{0h}\left(x_i\right)exp(-({\mu }_h\left(x_i\right)+w\left(x_i\right))t)>0$$

Hence, $R_h(t,x_i)>0$ for any arbitrary $x_i$. In a similar manner, we have

$$\frac{d{\mathit{S}}_v}{dt}+\left[\sum _{i=0}^L\frac{b{\lambda }_v{I}_h\left(t\right))}{1+{\nu }_v{I}_h\left(t\right)}+{\mu }_v\right]S_v\left(t\right)\ge 0$$

so that

$$\frac{d}{dt}\left[S_v\left(t\right)exp\left({\int }_0^t\frac{b{\lambda }_v{I}_h\left(\eta \right))}{1+{\nu }_v{I}_h\left(\eta \right)}d\left(\eta \right)+{\mu }_{v}t\right)\right]\ge 0$$

(2.9)

Integrating (2.9), we have

$$S_v\left(t\right)\ge S_{0v}exp\left[-\left({\int }_0^t\frac{b{\lambda }_v{I}_h\left(\eta \right))}{1+{\nu }_v{I}_h\left(\eta \right)}d\left(\eta \right)+{\mu }_{v}t\right)\right]>0\forall t>0$$

Also we have

$$\frac{d{\mathit{E}}_v}{dt}\ge -({\alpha }_v+{\mu }_v)E_v\left(t\right)$$

which on integration gives

$$E_v\left(t\right)\ge E_v\left(0\right)exp\left[-({\alpha }_v+{\mu }_v)t\right]>0\forall t>0$$

(2.10)

And finally, we have

$$\frac{d{\mathit{I}}_v}{dt}+({\mu }_v+{\gamma }_v)I_v\left(t\right)$$

so that

$$\frac{d}{dt}\left[I_v\left(t\right)exp({\mu }_v+{\gamma }_v)t\right]\ge 0$$

(2.11)

And we have

$$I_v\left(t\right)\ge I_v\left(0\right)exp\left[-({\mu }_v+{\gamma }_v)t\right]>0,\forall t>0$$

This completes the proof □

3. Existence and stability of the equilibrium points

3.1. Disease-free equilibrium

The disease-free equilibrium (DFE) points are steady state solutions that depict the absence of infection in both the human host and black-fly vector populations, i.e, onchocerciasis does not exist in the population. Thus, the disease-free equilibrium point, $E_0$, for the model (2.1) implies that $S^{*}{\left(x_i\right)}_h\ne 0$, $E_h^{*}\left(x_i\right)=I_h^{*}=0\left(x_i\right)=R_h^{*}\left(x_i\right)=0$, $S_v^{*}\ne 0$, $E_v=I_v=0$ and putting these into (2.1), we have $S^{*}{\left(x_i\right)}_h=\frac{{\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$ and $S_v^{*}=\frac{{\mathsf{\Psi }}_v}{{\mu }_v}$. Consequently we obtain $E_0$ as

$$E_0=\left(\frac{{\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)},0,0,0,\frac{{\mathsf{\Psi }}_v}{{\mu }_v},0,0\right)$$

(3.1)

A key notion in the analysis of infectious disease models is the basic reproduction number ${\mathcal{R}}_0$ , an epidemiological threshold that determines whether disease dies out or persists in the population.The basic reproduction number ${\mathcal{R}}_0$ of the system (2.1) is computed using the next generation matrix method and is given by

$${\mathcal{R}}_0=\sqrt{{\mathcal{R}}_h{\mathcal{R}}_v}$$

where ${\mathcal{R}}_h={\sum }_{i=0}^L\frac{\delta {\alpha }_h{\lambda }_h\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))(r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))}$ and ${\mathcal{R}}_v=\frac{\delta {\alpha }_v{\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_v({\alpha }_v+{\mu }_v)({\gamma }_v+{\mu }_v)}$. The basic reproduction number ${\mathcal{R}}_0$, determines whether onchocerciasis dies out or persists in the population. Therefore, ${\mathcal{R}}_h$ describes the number of humans that one infectious black-fly infects over its expected infectious period in a completely susceptible humans population, while ${\mathcal{R}}_v$ is the number of blac-flies infected by one infectious human during the period of infectiousness in a completely susceptible black-fly population.

3.2. Local Stability of the Disease-free Equilibrium Point $E_0$

Using the basic reproduction number obtained for the model (2.1), we analyse the stability of the equilibrium point in the following result.

Theorem 3:The disease-free equilibrium point, $E_0$, is locally asymptotically stable if ${\mathcal{R}}_0<1$, and unstable if ${\mathcal{R}}_0>1$.

Proof: The Jacobian matrix of the system (2.1) evaluated at the disease-free equilibrium point $E_0$, is obtained as

$$M\left(E_0\right)=\left(\begin{array}{ccccccc}{M}_{11}& 0& 0& M_{14}& 0& 0& M_{17}\\ 0& M_{22}& 0& 0& 0& 0& M_{27}\\ 0& M_{32}& M_{33}& 0& 0& 0& 0\\ 0& 0& M_{43}& M_{44}& 0& 0& 0\\ 0& 0& M_{53}& 0& M_{55}& 0& 0\\ 0& 0& M_{63}& 0& 0& M_{66}& 0\\ 0& 0& 0& 0& 0& M_{76}& M_{77}\end{array}\right)$$

where $M_{11}=-{\mu }_h\left(x_i\right)$, $M_{14}=w\left(a_1\right)$, $M_{17}=-{\sum }_{i=0}^L\frac{\delta {\lambda }_h\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$, $M_{22}=-({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))$, $M_{27}={\sum }_{i=0}^L\frac{\delta {\lambda }_h\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$, $M_{32}={\alpha }_h\left(x_i\right)$, $M_{33}=-(r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))$, $M_{43}=r\left(x_i\right)$, $M_{44}=-({\mu }_h\left(x_i\right)+w\left(x_i\right))$, $M_{53}=-\frac{\delta {\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_v}$, $M_{55}=-{\mu }_v$, $M_{63}=\frac{\delta {\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_v}$, $M_{66}=-({\alpha }_v+{\mu }_v)$, $M_{76}={\alpha }_v$, $M_{77}=-({\mu }_v+{\gamma }_v)$ We need to show that all the eigenvalues of $M\left(E_0\right)$ are negative. As the first and fifth columns form the two negative eigenvalues, $-{\mu }_h\left(x_i\right)$ and $-{\mu }_v$, the other five eigenvalues can be obtained from the sub-matrix, $M^1\left(E_0\right)$, formed by excluding the first and fifth rows and columns of $M\left(E_0\right)$. Hence

$$M^1\left(E_0\right)=\left(\begin{array}{ccccc}{M}_{11}^{\prime }& 0& 0& 0& M_{15}^{\prime }\\ {\alpha }_h\left(x_i\right)& M_{22}^{\prime }& 0& 0& 0\\ 0& r\left(x_i\right)& M_{33}^{\prime }& 0& 0& 0\\ 0& \frac{\delta {\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_v}& 0& -({\alpha }_v+{\mu }_v)& 0\\ 0& 0& 0& {\alpha }_v& -({\mu }_v+{\lambda }_v)\end{array}\right)$$

In the same way, the third column of $M^1\left(E_0\right)$ contains only the diagonal term which forms a negative eigenvalue, $-({\mu }_h\left(x_i\right)+w\left(x_i\right))$. The remaining four eigenvalues are obtained from the sub-matrix $M^2\left(E_0\right)$ given by

$$M^2\left(E_0\right)=\left(\begin{array}{cccc}{M}_{11}^{″}& 0& 0& M_{14}^{″}\\ {\alpha }_h\left(x_i\right)& M_{22}^{″}& 0& 0\\ 0& \frac{\delta {\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_v}& -({\alpha }_v+{\mu }_v)& 0\\ 0& 0& {\alpha }_v& -({\mu }_v+{\lambda }_v)\end{array}\right)$$

Thus, the eigenvalues of the matrix $M^2\left(E_0\right)$ are the roots of the characteristic equation of the form

$$(\xi +{\alpha }_h\left(x_i\right))(\xi +r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))(\xi +{\mu }_v+\gamma )-\sum _{i=0}^L\frac{{\delta }^2{\alpha }_h\left(x_i\right){\lambda }_h\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right){\alpha }_v{\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_h\left(x_i\right){\mu }_v}=0$$

(3.2)

If we let $Y_1={\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right)$, $Y_2=r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right)$, $Y_3={\alpha }_v+{\mu }_v$, and $Y_4={\mu }_v+{\gamma }_v$, then (3.2) becomes

$$X_4{\xi }^4+X_3{\xi }^3+X_2{\xi }^2+X_1\xi +X_0=0,$$

(3.3)

where

$$\left.\begin{array}{c}{X}_4=1\hfill \\ X_3=Y_1+Y_2+Y_3+Y_4\hfill \\ X_2=(Y_1+Y_2)(Y_2+Y_4)+Y_1{Y}_2+Y_3{Y}_4\hfill \\ X_1=(Y_1+Y_2)Y_3{Y}_4+(Y_3+Y_4)Y_1{Y}_2\hfill \\ X_0=Y_1{Y}_2{Y}_3{Y}_4-\sum _{i=0}^L\frac{{\delta }^2{\alpha }_h\left(x_i\right){\lambda }_h\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right){\alpha }_v{\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_h\left(x_i\right){\mu }_v}\hfill \end{array}\right\}$$

(3.4)

Expressing $X_0$ in terms of reproduction number ${\mathcal{R}}_0$, we have

$$X_0=Y_1{Y}_2{Y}_3{Y}_4(1-{\mathcal{R}}_0^2)$$

(3.5)

We can see from (3.4) that $X_1>0$, $X_2>0$, $X_3>0$, $X_4>0$, since all $Y_{i}s$ are positive. Moreover, if ${\mathcal{R}}_0<1$, it follows from (3.5) that $X_0>0$. Thus, using the Routh-Hurwitz criterion, we have $H_1=X_3>0$

$$H_2=\left|\begin{array}{cc}{X}_3& X_4\\ X_1& X_2\end{array}\right|=Y_1(Y_2+Y_3+Y_4)(Y_1+Y_2+Y_3+Y_4)+(Y_2+Y_3)(Y_2+Y_4)(Y_3+Y_4)>0$$

Similarly we have $H_3>0$ and $H_4>0$ where

$$H_3=\left|\begin{array}{ccc}{X}_3& X_4& 0\\ X_1& X_2& X_3\\ 0& X_0& X_1\end{array}\right|and{H}_4=\left|\begin{array}{cccc}{X}_3& X_4& 0& 0\\ X_1& X_2& X_3& X_4\\ 0& X_0& X_1& X_2\\ 0& 0& 0& X_0\end{array}\right|$$

Therefore, all the eigenvalues of the Jacobian matrix $M\left(E_0\right)$ have negative real parts when ${\mathcal{R}}_0<1$ and the disease-free equilibrium point is locally asymptotically stable. However, when ${\mathcal{R}}_0>1$, we see that $X_0<0$ and there is one eigenvalue with positive real part and therefore the disease-free equilibrium point is unstable □

3.3. Endemic Equilibrium Point $E_e$

We shall show that the formulated model (2.1) has an endemic equilibrium point, $E_e$. The endemic equilibrium point is a positive steady state solution where the disease persists in the population.

Theorem 4: The model (2.1) has a unique endemic equilibrium $E_e$ whenever ${\mathcal{R}}_0>1$.

Proof: Let $E_e=(S_h^{″}\left(x_i\right),E_h^{″}\left(x_i\right),I_h^{″}\left(x_i\right),R_h^{″}\left(x_i\right),S_v^{″},E_v^{″},I_v^{″})$ be a nontrivial equilibrium of the model (2.1). That is, all components of $E_e$ are positive. Then the onchocerciasis model (2.1) at steady-state becomes

$${\mathsf{\Psi }}_h\left(x_i\right)-\sum _{i=0}^L\left(\frac{\delta {\lambda }_h\left(x_i\right)S_h^{″}\left(x_i\right)I_v}{1+{\nu }_h\left(x_i\right)I_v^{″}}-{\mu }_h\left(x_i\right)S_h^{″}\left(x_i\right)+\omega \left(x_i\right)R_h^{″}\left(x_i\right)\right)=0$$

(3.6)

$$\sum _{i=0}^L\left(\frac{\delta {\lambda }_h\left(x_i\right)S_h^{″}\left(x_i\right)I_v}{1+{\nu }_h\left(x_i\right)I_v^{″}}-({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))E_h^{″}\left(x_i\right)\right)=0$$

(3.7)

$$\sum _{i=0}^L({\alpha }_h\left(x_i\right)E_h^{″}\left(x_i\right)-(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h\left(x_i\right))I_h^{″}\left(x_i\right))=0$$

(3.8)

$$\sum _{i=0}^{L}r\left(x_i\right)I_h^{″}\left(x_i\right)-({\mu }_h\left(x_i\right)+\omega \left(x_i\right))R_h^{″}\left(x_i\right)=0$$

(3.9)

$${\mathsf{\Psi }}_v-\frac{\delta {\lambda }_v{S}_v^{″}{I}_h\left(x_i\right)}{1+{\nu }_v\left(x_i\right)I_h^{″}\left(x_i\right)}-{\mu }_v{S}_v^{″}=0$$

(3.10)

$$\frac{\delta {\lambda }_v{S}_v^{″}{I}_h\left(x_i\right)}{1+{\nu }_v\left(x_i\right)I_h^{″}\left(x_i\right)}-({\alpha }_v+{\mu }_v)E_v^{″}=0$$

(3.11)

$${\alpha }_v{E}_v^{″}-({\mu }_v+{\gamma }_v)I_v^{″}=0$$

(3.12)

From the last three equations, we have

$$I_v^{″}=\frac{{\alpha }_v{E}_v^{″}}{{\mu }_v+{\gamma }_v}$$

(3.13)

$$E_v^{″}=\frac{\delta {\lambda }_v{S}_v^{″}{I}_h\left(x_i\right)}{1+{\nu }_v\left(x_i\right)I_h^{″}\left(x_i\right)({\alpha }_v+{\mu }_v)}$$

(3.14)

and

$$S_v^{″}=\frac{{\mathsf{\Psi }}_v}{\frac{\delta {\lambda }_v{S}_v^{″}{I}_h\left(x_i\right)}{1+{\nu }_v\left(x_i\right)I_h^{″}\left(x_i\right)}+{\mu }_v}$$

(3.15)

Substituting (3.14) and (3.15) into (3.13) yields

$$I_v^{″}=\frac{{\mathcal{R}}_v{\mu }_v{I}_h^{″}\left(x_i\right)}{{\mu }_v+(\delta {\lambda }_v+{\mu }_v{\nu }_v)I_h^{″}\left(x_i\right)}$$

(3.16)

From (3.8) and (3.9), we have

$$E_h^{″}\left(x_i\right)=\sum _{i=0}^L\frac{(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h\left(x_i\right))I_h\left(x_i\right)}{{\alpha }_h\left(x_i\right)}$$

(3.17)

and

$$R_h^{″}\left(x_i\right)=\sum _{i=0}^L\frac{r\left(x_i\right)I^{″}\left(x_i\right)}{{\mu }_h\left(x_i\right)+\omega \left(x_i\right)}$$

(3.18)

If we put (3.16) and (3,17) in (3.7) in terms of ${\mathcal{R}}_0$, we have

$$S_h^{″}\left(x_i\right)=\frac{{\sum }_{i=0}^L{\mathsf{\Psi }}_h\left(x_i\right)[{\mu }_v+(\delta {\lambda }_v+{\mu }_v{\nu }_v+{\nu }_h\left(x_i\right){\mu }_v{\mathcal{R}}_v)I_h^{″}\left(x_i\right)]}{{\mu }_h\left(x_i\right){\mu }_v{\mathcal{R}}_0^2}$$

(3.19)

Finally, using (3.16), (3.18) and (3.19) in (3.7), we have

$$I_h^{″}\left(x_i\right)=\sum _{i=0}^L{\mu }_h\left(x_i\right){\mu }_v{\mathsf{\Psi }}_h\left(x_i\right)({\mu }_h\left(x_i\right)+\omega \left(x_i\right))\frac{({\mathcal{R}}_0^2-1)}{\rho }$$

(3.20)

where

$$\rho =\sum _{i=0}^L({\mu }_h\left(x_i\right)+\omega \left(x_i\right))[\delta {\lambda }_h\left(x_i\right){\mu }_v{\mathcal{R}}_v+{\mathsf{\Psi }}_h\left(x_i{\mu }_h\left(x_i\right)(\delta {\lambda }_v+{\mu }_v{\nu }_v+{\nu }_h\left(x_i\right){\mu }_v){\mathcal{R}}_m\right)]-\sum _{i=0}^L{\mu }_h\left(x_i\right){\mu }_v\omega \left(x_i\right)r\left(x_i\right){\mathcal{R}}_0^2.$$

(3.21)

If in (3.20), $\omega \left(x_i\right)=0$ then $\rho >0$. From this, one sees that model (2.1) has no positive solution when ${\mathcal{R}}_0<1$. However, with $\omega \left(x_i\right)=0$, a unique endemic equilibrium exists when ${\mathcal{R}}_0>1$. This completes the proof. □

Remark 1:It is important to have a remark that positive solution exists for the model (2.1) in a case where $\rho <0$ and ${\mathcal{R}}_0<1$. This implies that the disease-free equilibrium co-exists with the endemic equilibrium state when ${\mathcal{R}}_0$ is slightly less than unity resulting into a phenomenon of subcritical (backward) bifurcation.

4. Bifurcation Analysis

The mathematical examination of changes in the qualitative behaviour of a dynamical system as its parameter passes through a critical value called a bifurcation point is referred to as bifurcation analysis.

Theorem 5(Castillo-Chavez and Song (2004):

Consider the following general system of ordinary differential equations with a parameter $\phi$:

$$\frac{dy}{dt}=f(y,\phi ),f:{\mathbb{R}}^n\times \mathbb{R}⟶\mathbb{R}\mathrm{and}f\in \mathbb{C}({\mathbb{R}}^n\times \mathbb{R})$$

(4.1)

where 0 is an equilibrium point of the system (that is, $f(0,\phi )\equiv 0$ for all $\phi$) and we have the following assumptions:

1.

$A=D_{y}f(0,0)=\left(D_{y}f(0,0),\right)$ is the linearization matrix of the system (4.1) around the equilibrium 0 with $\phi$ evaluated at 0. Zero is a simple eigenvalue of A and other eigenvalues of A have negative real parts;

2.

Matrix A has a nonnegative right eigenvector $\kappa$ and a left eigenvector $\eta$ corresponding to the zero eigenvalue.

Let $f_k$ be the $kth$ component of f and

$$m=\sum _{k,i,j=1}^n{\eta }_k{\kappa }_i{\kappa }_j\frac{{\partial }^2{f}_k}{\partial y_i\partial y_j}(0,0),$$

$$n=\sum _{i,j=1}^n{\eta }_k{\kappa }_i\frac{{\partial }^2{f}_k}{\partial y_i\partial \phi }(0,0)$$

The local dynamics of the system around 0 is totally determined by the signs of m and n.

(i)

$m>0$, $n>0$. When $\phi <0$ with $|\phi <<1$, 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when $0<\phi <<1$, 0 is unstable and there exists a negative, locally asymptotically stable equilibrium;

(ii)

$m<0$, $n<0$. When $\phi <0$ with $\left|\phi \right|1$, 0 is unstable; when $0<\phi <<$1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium;

(iii)

$m>0$, $n<0$. When $\phi <0$ with $\left|\phi \right|<<1$, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when $0<\phi <<1$, 0 is stable, and a positive unstable equilibrium appears;

(iv)

$m<0$, $n>0$. When $\phi$ 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 particular, if $m>0$ and $n>0$, then there exists a backward bifurcation.

To demonstrate the possibility of the co-existence of the equilibria of model (2.1) at ${\mathcal{R}}_0<1$ but near ${\mathcal{R}}_0=0$, the Center Manifold Theory described by Castillo-Chavez and Song (2004). Let the onchocerciasis model (2.1) be written in the vector form

$$\frac{dY}{dt}=F\left(Y\right)$$

where $Y={(y_1,y_2,y_3,y_4,y_5,y_6,y_7)}^T$ and $F={(f_1,f_2,f_3,f_4,f_5,f_6,f_7)}^T$ so that so that $Sh(t;x_i)=y_1$; $E_h(t;x_i)=y_2$; $Ih(t;x_i)=y_3$; $R_h(t;x_i)=y_4$; $S_v\left(t\right)=y_5$; $E_v\left(t\right)=y_6$; and $I_v\left(t\right)=y_7$: Then model (2.1) becomes

(4.2)

Let ${\lambda }_h\left(x_i\right)$ be a bifurcation parameter such that ${\lambda }_h\left(x_i\right)={\lambda }_h^{*}\left(x_i\right)$ when ${\mathcal{R}}_0=1$. Then

$${\lambda }_h^{*}\left(x_i\right)=\sum _{i=0}^L\frac{{\mu }_h\left(x_i\right){\mu }_v({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))(r\left(x_i\right)+\mu {\left(x_i\right)}_h+{\gamma }_h\left(x_i\right))({\mu }_v+\gamma )({\alpha }_v+{\mu }_v)}{{\delta }^2{\alpha }_h\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right){\alpha }_v{\lambda }_v{\mathsf{\Psi }}_v}$$

The linearization matrix $M(E_0;{\lambda }_h^{*}\left(x_i\right))$ of the model (2.1), evaluated at the disease-free equilibrium $E_0$ and the bifurcation parameter ${\lambda }_h^{*}\left(x_i\right)$, is given

$$M(E_0;{\lambda }_h^{*}\left(x_i\right))=\left(\begin{array}{ccccccc}-{\mu }_h\left(x_i\right)& 0& 0& \omega \left(x_i\right)& 0& 0& A\\ 0& B& 0& 0& 0& 0& C\\ 0& {\alpha }_h\left(x_i\right)& D& 0& 0& 0& 0\\ 0& 0& r\left(x_i\right)& E& 0& 0& 0\\ 0& 0& G& 0& -{\mu }_v& 0& 0\\ 0& 0& K& 0& 0& L^{\prime }& 0\\ 0& 0& 0& 0& 0& {\alpha }_v& M\end{array}\right)$$

where $A=-{\sum }_{i-0}^L\frac{\delta {\lambda }_h^{*}\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$, $B=-({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))$, $C={\sum }_{i-0}^L\frac{\delta {\lambda }_h^{*}\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$, $D=-(r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))$, $E=-({\mu }_h\left(x_i\right)+\omega \left(x_i\right))$, $G=-\frac{\delta {\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_v}$, $K=\frac{\delta {\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_v}$, $L^{\prime }=-({\alpha }_v+{\mu }_v)$, $M=-({\mu }_v+{\gamma }_v)$

The eigenvalues of $M(E_0,{\lambda }_h\left(x_i\right))$ can be obtained by solving the characteristic equation

$$(\xi +{\mu }_h\left(x_i\right))(\xi +{\mu }_h)(\xi +{\mu }_h\left(x_i\right)+\omega \left(x_i\right))P_1\left(\xi \right)=0,$$

where $P_1\left(\xi \right)$ is a polynomial of degree four. Consequently, $M(E_0,{\lambda }^{*}\left(x_i\right))$ has a zero eigenvalue which is simple and other eigenvalues are real and negative. Let the right eigenvector corresponding to this zero eigenvalue be denoted by $\kappa ={\left({\kappa }_1,{\kappa }_2,{\kappa }_3,{\kappa }_4,{\kappa }_5,{\kappa }_6,{\kappa }_7\right)}^7$ so that $M(E_0,{\lambda }_h^{*}\left(x_i\right))\kappa =0$. Then we have

$$\left.\begin{array}{c}-{\mu }_h\left(x_i\right){\kappa }_1+\omega \left(x_i\right){\kappa }_4+A{\kappa }_7=0\hfill \\ B{\kappa }_2-A{\kappa }_7=0\hfill \\ {\alpha }_h\left(x_i\right){\kappa }_2+D{\kappa }_3=0\hfill \\ r\left(x_i\right){\kappa }_3+E{\kappa }_4=0\hfill \\ G{\kappa }_3-{\mu }_v{\kappa }_5=0\hfill \\ -G{\kappa }_3+L^{\prime }{\kappa }_6=0\hfill \\ {\alpha }_v{\kappa }_6+M{\kappa }_7=0\hfill \end{array}\right\}$$

(4.3)

From (4.2), the components of the right eigenvector $\kappa$ are given by

$$\left.\begin{array}{c}{\kappa }_1=U{\kappa }_6\hfill \\ {\kappa }_2=\frac{(r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))({\alpha }_v+{\mu }_v){\mu }_v{\kappa }_6}{\delta {\lambda }_v{\mathsf{\Psi }}_v\alpha \left(x_i\right)}\hfill \\ {\kappa }_3=\frac{({\alpha }_v+{\mu }_v){\mu }_v{\kappa }_6}{\delta {\lambda }_v{\mathsf{\Psi }}_v}\hfill \\ {\kappa }_4={\frac{r\left(x_i\right)({\alpha }_v+{\mu }_v){\mu }_v}{\kappa }}_6\delta {\lambda }_v{\mathsf{\Psi }}_v({\mu }_h\left(x_i\right)+\omega \left(x_i\right))\hfill \\ {\kappa }_5=-\frac{({\alpha }_v+{\mu }_v){\kappa }_6}{{\mu }_v}\hfill \\ {\kappa }_6=-\frac{1}{T_1}\hfill \\ {\kappa }_7=\frac{{\alpha }_v{\kappa }_6}{{\mu }_v+{\gamma }_v}\hfill \end{array}\right\}$$

(4.4)

$$\begin{array}{c}U=-\sum _{i=0}^L\frac{({\alpha }_v+{\mu }_v){\mu }_v}{\delta {\lambda }_v{\mathsf{\Psi }}_v{\alpha }_h\left(x_i\right){\mu }_h\left(x_i\right)({\mu }_h\left(x_i\right)+\omega \left(x_i\right))}(r\left(x_i\right){\alpha }_h\left(x_i\right)({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right)+\omega \left(x_i\right))\hfill \\ \hfill +({\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))(r\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_h\left(x_i\right))({\mu }_h\left(x_i\right)+\omega \left(x_i\right))\end{array}$$

(4.5)

and

$$\begin{array}{c}{T}_1={\alpha }_v[({\alpha }_v+{\mu }_v)({\mu }_v+{\gamma }_v)(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h\left(x_i\right)+{\mu }_v({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right)))\hfill \\ \hfill +\sum _{i=0}^L({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h\left(x_i\right))(2{\mu }_v+{\gamma }_v+{\alpha }_v)]\end{array}$$

(4.6)

Further, $M(E_0,{\lambda }^{*}\left(x_i\right))$ has a left eigenvector, $\eta =({\eta }_1,{\eta }_2,{\eta }_3,{\eta }_4,{\eta }_5,{\eta }_6,{\eta }_7)$, associated with the zero eigenvalue, which can be obtained from $\eta M(E_0,{\lambda }^{*}\left(x_i\right))=0$ as follows: $=-\left(\alpha h\right(ai)+\mu h(ai\left)\right)\left(r\right(ai)+\mu h(ai)+\delta h(ai\left)\right)(\alpha m+\mu m)(\mu m+\delta m)$

$$\left.\begin{array}{c}{\eta }_1=0\hfill \\ {\eta }_2=\frac{\delta {\lambda }_v{\alpha }_v{\alpha }_h\left(x_i\right){\mathsf{\Psi }}_v{\eta }_7}{({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h\left(x_i\right))({\alpha }_v+{\mu }_v)}\hfill \\ {\eta }_3=\frac{\delta {\lambda }_v{\alpha }_v{\mathsf{\Psi }}_v{\eta }_7}{(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h\left(x_i\right))({\alpha }_v+{\mu }_v)}\hfill \\ {\eta }_4=0\hfill \\ {\eta }_5=0\hfill \\ {\eta }_6=\frac{{\alpha }_v{\eta }_7}{{\alpha }_v}+{\mu }_v\hfill \\ {\eta }_7=-({\alpha }_h\left(x_i\right)+{\mu }_h\left(a_i\right))(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h\left(x_i\right))({\alpha }_v+{\gamma }_v)\hfill \end{array}\right\}$$

(4.7)

We see from (4.1) that all the second-order partial derivatives at $E_0$ and $\lambda \left(x_i\right)$ are zero except the following:

$$\frac{{\partial }^2{f}_1}{\partial y_1\partial y_7}=\frac{{\partial }^2{f}_1}{\partial y_7\partial y_1}=-\delta {\lambda }_h^{*}\left(x_i\right),$$

$$\frac{{\partial }^2{f}_2}{\partial y_1\partial y_7}=\frac{{\partial }^2{f}_2}{\partial y_7\partial y_1}=\delta {\lambda }_h^{*}\left(x_i\right)$$

,

$$\frac{{\partial }^2{f}_5}{\partial y_3\partial y_5}=-\frac{{\partial }^2{f}_5}{\partial y_5\partial y_3}=-\delta {\lambda }_v,\frac{{\partial }^2{f}_6}{\partial y_3\partial y_5}=\frac{{\partial }^2{f}_6}{\partial y_5\partial y_3}=\delta {\lambda }_v$$

$$\frac{{\partial }^2{f}_1}{d{y}_7^2}=\sum _{i=0}^L\frac{2\delta {\lambda }_h\left(x_i\right){\nu }_h\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)},\frac{{\partial }^2{f}_2}{d{y}_7^2}=-\sum _{i=0}^L\frac{2\delta {\lambda }_h\left(x_i\right){\nu }_h\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$$

$$\frac{{\partial }^2{f}_5}{\partial y_3^2}=\frac{2\delta {\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_v},\frac{{\partial }^2{f}_6}{\partial y_3^2}=\frac{2\delta {\lambda }_v{\mathsf{\Psi }}_v}{{\mu }_v}$$

tegether with

$$\frac{{\partial }^2{f}_1}{{\partial }_7{\partial }_h\left(x_i\right)}=-\delta \frac{{\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$$

and

$$\frac{{\partial }^2{f}_2}{{\partial }_7{\partial }_h\left(x_i\right)}=b\frac{{\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}$$

If we let the bifurcation coefficients m and n of Theorem be $m_1$ and $n_1$ respectively, it then follows from the above partial derivatives that

$$\begin{array}{c}{m}_1=\sum _{i,j,k=1}^7{\eta }_k{\kappa }_i{\kappa }_j\frac{{\partial }^2{f}_k}{\partial y_i\partial y_j}(E_0,{\lambda }_h^{*}\left(x_i\right))\hfill \\ =2{\eta }_1{\kappa }_1{\kappa }_7\frac{{\partial }^2{f}_1}{\partial y_1\partial y_7}+{\eta }_1{\kappa }_7^2\frac{{\partial }^2{f}_1}{{\partial }^2{y}_j}+2{\eta }_2{\kappa }_1{\kappa }_7\frac{{\partial }^2{f}_2}{\partial y_1\partial y_7}\\ +{\eta }_2{\kappa }_7^2\frac{\partial f_2}{\partial y_7^2}+2{\eta }_5{\kappa }_3{\kappa }_5\frac{\partial f_5}{\partial y_3\partial y_5}+{\eta }_5{\kappa }_3^2\frac{{\partial }^2{f}_5}{\partial y_3^2}\\ \hfill +2{\eta }_6{\kappa }_3{\kappa }_5\frac{\partial f_6}{\partial y_3\partial y_5}+{\eta }_6{\kappa }_3^2\frac{{\partial }^2{f}_6}{{\partial }^2{y}_3}\end{array}$$

(4.8)

Substituting the appropriate right eigenvector (4.3), the left eigenvector (4.6) and the second-order partial derivatives obtained earlier into $m_1$ yields

$$\begin{array}{c}{m}_1=\sum _{i=0}^L\frac{{\delta }^2{\lambda }_h^{*}\left(x_i\right){\lambda }_v{\alpha }_v^2{\alpha }_h\left(x_i\right){\mathsf{\Psi }}_v{\eta }_7{\kappa }_6^{2}U}{({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h)({\alpha }_v+{\mu }_v)({\mu }_v+{\gamma }_v)}\hfill \\ +\sum _{i=0}^L\frac{\delta {\lambda }_v{\alpha }_v^2{\alpha }_h\left(x_i\right){\mathsf{\Psi }}_v{\eta }_7{\kappa }_6^2}{({\alpha }_h\left(x_i\right)+{\mu }_h\left(x_i\right))(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h\left(x_i\right))({\alpha }_v+{\mu }_v){({\mu }_v+{\gamma }_v)}^2}\\ \hfill \times \left[\frac{-2\delta {\lambda }_h^{*}\left(x_i\right){\nu }_h\left(x_i\right){\mathsf{\Psi }}_h\left(x_i\right)}{{\mu }_h\left(x_i\right)}\right]-\frac{2{\alpha }_v({\alpha }_v+{\mu }_v){\eta }_7{\kappa }_6^2}{{\mathsf{\Psi }}_v}-\frac{{\alpha }_v{\mu }_v({\alpha }_v+{\mu }_v){\nu }_v{\eta }_7{\kappa }_6^2}{\delta {\lambda }_v{\mathsf{\Psi }}_v}\end{array}$$

(4.9)

Similarly, we obtain

$$\begin{array}{c}{n}_1=\sum _{k,i=0}^7{\eta }_k{\kappa }_i\frac{{\partial }^2{f}_k}{\partial y_i\partial {\lambda }_h}(E_0,{\lambda }_h^{*}\left(x_i\right))=2{\eta }_1{\kappa }_7\frac{{\partial }^2{f}_1}{\partial y_7\partial {\lambda }_h}+2{\eta }_2{\kappa }_7\frac{{\partial }^2{f}_2}{\partial y_7\partial {\lambda }_h\left(x_i\right)}\hfill \\ \hfill =\sum _{i=0}^L\frac{2{\delta }^2{\alpha }_h\left(x_i\right){\alpha }_v^2{\lambda }_v{\mathsf{\Psi }}_v{\mathsf{\Psi }}_h\left(x_i\right){\eta }_7{\kappa }_6}{({\mu }_v+{\gamma }_v)\left({\alpha }_h\left(x_i\right){\mu }_h\left(x_i\right)\right)(r\left(x_i\right)+{\mu }_h\left(x_i\right)+{\gamma }_h\left(x_i\right)({\alpha }_v+{\mu }_v){\mu }_h\left(x_i\right))}\end{array}$$

(4.10)

The signs of $m_1$ and $n1$ are important in determining the direction of bifurcation at $mathcalR=1$ (Castillo-Chavez and Song, 2004). One sees that $m_1>0$ and $n_1>0$ and by item(iv) of Theorem 5, the model (2.1) is capable of exhibiting a backward bifurcation phenomenon. Thus, we have the following result:

Theorem 6: The onchocerciasis model given by (2.1) undergoes a phenomenon of backward bifurcation at $E_0$ and ${\mathcal{R}}_0=1$.