A Stochastic, Individual-Based Mathematical Model of Onchocerciasis Transmission in North-Eastern Nigeria

1. Introduction

Onchocerciasis, also known as river blindness is a tropical disease caused by the filarial nematode Onchocerca volvulus and transmitted through bites of Simulium blackflies,(maitumbi in Hausa, duduru in Kam language) . It is a disease of public importance and is highly associated with disability (Torsen et al., 2025). Onchocerciasis is one of the neglected torpical diseases and perhaps the most studied filarial diseases in Nigeria (Akogun & Onwuliri, 1991; Asha & Nyimvua, 2024). About seven million Nigerians are estimated to be infected while another forty two million are at the risk of infection. Each year, new foci of onchocerciasis transmission are discovered, revealing that, its distribution is far more extensive than has been assumed. The most notable foci of transmission in the central and northern Nigeria have been identified as Abuja FCT, Borno, Benue, Plateau, Nasarawa, Kogi, Adamawa and Taraba states (Akogun, 1992; Johanns et al., 2021).

Literatures showed that thirty seven (37) nations have the present of the disease. Of these nations, thirty (30) are African countries while six (6) and one (1) are in the Americas and Arabian Peninsula respectively. Onwubuya et al., (2023) confirmed that the estimated incidence revealed 500,000 for people with a visual impairment and 270,000 for people who are blind and that it is the second most common preventable cause of blindness after trachoma in sub-Saharan Africa. Onchocerca volvulus is carried by blackflies (Simulium blackflies), which spawn in fast-moving, highly oxygenated rivers and streams. An infected black fly releases one or more Onchocerca volvulus larvae into the human host when it feeds on the blood of an uninfected person. These larvae mature into fully fledged adult worms in less than a year. Usually over bony prominences, they often aggregate to create fibrous nodules beneath the skin. Onchocerciasis is currently mostly controlled in Africa through annual mass drug administration (MDA) with ivermectin. Periodic community wide mass drug administration (MDA) with ivermectin (Mectizan, Merck) can stop the spread of infection and prevent eye and skin illness. Significant progress has been made in the last three decades to control onchocerciasis in both Africa and the Americas. This is primarily due to global public private partnerships, ongoing funding from pertinent organizations and stakeholders, the creation of new tools and technology, and regional programs that have shifted from eradicating the disease to, whenever possible, eradicating the infection. However, given the current endemic situation of onchocerciasis in Africa, the MDA approach might not be sufficient to stop the disease's spread even in areas where vector control has been successful because the disease can be brought back into the area by migrant infected individuals (Onwubuya et al., 2023).

There have been endemic foci of onchocerciasis in Nigeria, with high prevalence of blindness and other complications associated with it compared with the other African nations (Akogun, 1992 in Idowu et al, 2013). Endemic communities of some of the states have been on annual ivermectin treatment since 1995. However, a survey by Idowu et al., (2013) reveals that the area still had incidences of the disease (though gradually becoming hypo endemic). Of concern also is the fact that about 27% of eligible persons were reported not to have been treated owing to various reasons including, fear of side effects of the drug; wrong perception by the people that the drug is not good; and absence of some people as at the time of drug administration. The appreciable proportion of those eligible for treatment that missed treatments overtime, erode the gains of ivermectin distribution (Onwujekwe et al., 2000; Awadzi et al., 2004). These people serve as reservoirs of the parasite and active transmission could be experienced in such communities where the potential vectors are available. Also, there are years that ivermectin distribution cannot be held in some areas as attested by researchers. The disease is adaptable in the riverine areas and has high rate of breeding in the fast flowing rivers. In most studies conducted, Nigeria has been the worst hit country in West and Africa at large while North-East Nigeria has the highest prevalence rate.

Johanns et al., (2022), reported that in the past 3 decades there has been a repeated annual treatments of onchocerciasis using ivermectin and that this has largely eliminated Onchocerca volvulus microfilariae (mf) from the skin and eyes of onchocerciasis patients. This treatment has halted the emergence of ocular pathologies which help to improving dermal health of those that were affected. The essence of all of these high coverage of ivermectin mass drug administration (MDA) to eligible populations, are aim to suppress the transmission of Onchocerca volvulus. Although the symptoms of the disease have largely disappeared, Ochocerca volvulus adult worms survive because single annual doses of ivermectin do not kill adult female Onchocerca volvulus and do not exert embryotoxic or embryostatic effects (Johanns et al., 2022). According to them, such medication does suppress the release of mf from gravid female Onchocerca volvulus for several months, and this measure must be applied repeatedly until fertile female Onchocreca volvulus either stop reproducing mf or die. In the area where this disease is high, it tend to affect the normal activities of the people living in that area, especially those living in the riverine area and if the spread of the disease is on the high rate of increase; their mode of livelihood suffers a great deal. The onchocerciasis which is a parasitic infectious disease, (the parasite Onchocerca volvulus), a filarial nematode lives in subcutaneous nodules, from where it produces and disseminate millions and millions of living embryos or microfilariae. Study has shown that the lifespan of the adult parasite is estimated at 10-11 years on average. The sessile female parasite can reach a length of 20-50 cm (Duke, 1990).

Onchocerciasis disease occurrences are more often in warm tropical settings of which Nigeria happen to be the tropical warm zone. The parasitic flies can survive in environments that arefavourable for their adaptation and growth all year round (Gbogbo et al., 2023). Onchocerciasis disease affects most of the people that live in rural areas. The rural populations suffer from disease and it has been a major cause of blindness, Onchocerciasis skin disease (OSD) and severe itching in the body of the patient that is affected with the disease. The endemic areas often suffer from the severe socioeconomic outcomes. The symptoms of the disease include severe itching, disfiguring skin conditions, and visual impairment globally. It is estimated that about 90 million people are at risk of infecting the disease (Gbogbo et al., 2023). Most of the study done on onchocerciasis employed mathematical modeling to predict and forecast the prevalence of the disease in an endemic areas.

A mathematical model is an abstract representation of a phenomenon created with the use of equations that generate perspectives of the general behaviour of an epidemic event, pointing out that the models can be improved mathematically to make them more similar to actual data. Mathematical model is a crude general behaviour of an epidemic as addressed by epidemic curves, which allows predictions about the duration of an epidemic, its magnitude in the population, and the evaluation of factors that influence the transmission dynamics, and consequently the number of cases (Costa et al., 2021). It is also a means of investigating the influence of determinate factors over disease spread. Mathematical modeling is a very versatile tool in the epidemiology of infectious diseases. It can be used to identify patterns in epidemics, extrapolate epidemic behaviours, and examine the effects of interventions like pharmaceutical treatment, vaccination, quarantine, social distancing, and hygiene measures in a dynamic setting. It is also inexpensive and allows for the simulation of experiments that are unethical in human subjects as well as experiments that have low economic viability in animal models. There are two broad categories of mathematical models used in infectious disease epidemiology: deterministic models that consider nonrandom rate flows in a population stratified in compartments; and stochastic models that take into account probabilities in the movements between the model's compartments, such as the likelihood that a susceptible person will become infected and the likelihood that the disease will spread throughout the population that a mathematical system is used to address (Costa et al., 2021).

Mathematical frameworks or model have been devised to investigate how these mechanisms affect parasite population dynamics and how control measures affect them. The assumptions about the adult female worm’s reproductive life duration and the distribution of its survival times at the longest living parasite stage are critical to all frameworks because they significantly impact post control forecasts. Creating a rational foundation for regulations intended to stop the spread of a disease is one of the goals of modeling epidemics (Onwubuya et al., 2023). The evaluation of the public health authorities’ intervention is made possible by including realistic optimal solutions in models. In an epidemic, decisions must be made on how much of the population should be treated to reduce both the spread of infection and the implementation costs of the treatment. Optimal control is a potent mathematical technique for making decisions that entail using the right tactics to eliminate infections from the populace (Onwubuya et al., 2023).

More often mathematical and computational models are ways used to understand the processes in complex systems that underlie in empirical observations and to generate possible future trajectories of these systems. In fact, mathematical models are the actual framework of composing a set of equations or theoretical approaches, while computational models refer to the numerical and computational approaches used to solve the mathematical framework. The overarching purpose of models is to allow a flexible framework in which different scenarios can be tested, different questions related to future behaviour can be posed and evaluated, and potential future behaviour can be predicted (Panovska-Griffiths et al., 2021).

Mathematical modeling seeks to create answers to unanswered questions from observations. Furthermore, the answers are then used to understand, manage and predict future behaviour ofcomplex systems and processes, for example, to inform public policy and future decision making. It is always important to bear in mind the assumptions underlying mathematical models and how these models are related to the real world scenarios. It is important to understand that modeling is like any other technology: it can be properly applied or not, it may produce output that admits a useful interpretation, or it may not (Panovska-Griffiths et al., 2021).

Rapid assessment of an outbreak of disease or natural disaster is made possible by mathematical modeling and simulation. When there are many experimental scenarios to test or when the expense of data collection is unaffordable, simulation is frequently utilized. Numerous methods have been put forth over the years to address the issue from various angles. These are divided into three main categories: mathematical models used in dynamical systems also known as state space models to predict the evolution of a hypothetical or ongoing epidemic spread; statistical methods for outbreak surveillance and the identification of spatial patterns in actual epidemics; and of recent the use of machine learning/expert methods for predicting the evolution of an ongoing epidemic. Again, there are various ways that weave a vast and varied literature for each of these three groups. This often depends in an attempt to map out these methods and explain the fundamental ideas that underlie these models (Constantinos & Lucia, 2013).

Mathematical modeling and optimization tools have been used rigorously for a assessing the epidemic of disease. However, there exists gaps that need to be addressesd. For example, despite the extensive application of mathematical models to study infectious diseases, including onchocerciasis, gaps remain in the literature. One critical challenge is that most of the worked done concerning the subject matter employed deterministic models but could not predict the prevalence of the disease. Secondly, there is lack of high-quality, region-specific data essential for calibrating and validating these models to ensure their accuracy and reliability. Additionally, most of the studies do not consider external factors such as climatic variables like temperature and rainfall being incorporated into the model. Also of important is that most of the work does not considered heterogeneity of the population but homogeneity. The available literatures reviewed do not consider a region such as this study.

2. Materials and Methods

The research design for this study is an explorative research design conducted across various communities that are residing in the riverine areas in the North-Eastern State of Nigeria viz: Adamawa, Bauchi, Borno, Gombe, Taraba and Yobe. The design of the Study sites include ,riverine areas, streams and vegetation cover as well as the climate and weather conditions of the study community settings. This reflects the diverse socio-economic, cultural, environmental, climate and healthcare factors influencing onchocerciasis ivermectin mass drug administration in the region. The design also considers onchocerciasis prevalence, accessibility, and data availability to ensure a thorough understanding of the impact of treatment and interventions.

1.4. Model Structure Assumptions

Most of the literatures reviewed in this work looked at the deterministic model with homogeneous populations; where the rate of transmission of the disease and the biting rate are considered uniform. This study looked at the heterogeneous nature of the populations of the study considering human, the blackfly and the environmental factors both having heterogeneous rate of transmission and biting rate. It is also assumed that the demographic characteristics and geographical locations are both heterogeneous. The assumptions of each are as seen below:

Table 1. Model Structure assumptions.

S/No	Assumptions	Meaning
1.	Stochastic model	The model output is randomly determined by the input
2.	Individual-based model	The model is individual-based, meaning that each human and blackfly is model separately
3.	Spatially-explicit model	The spatial location of each human and blackfly is taken into account
.4.	Time-dependent model	The model is time dependent, meaning that the output changes over time

Table 1. Model Structure assumptions.

S/No	Assumptions	Meaning
1.	Stochastic model	The model output is randomly determined by the input
2.	Individual-based model	The model is individual-based, meaning that each human and blackfly is model separately
3.	Spatially-explicit model	The spatial location of each human and blackfly is taken into account
.4.	Time-dependent model	The model is time dependent, meaning that the output changes over time

Table 2. Human Components assumptions.

S/No	Assumptions	Meaning
1.	Age structure	The human components have an age structure, with different age groups having different susceptibility to infection.
2.	Genetic variation	The human population is heterozygous, with different genotypes having different levels of susceptibility to infection.
3.	Heterogeneous mixing	Human mix non-randomly, with different age groups and genotypes having different contact rates.
4.	Spatial structure	The human population is spatially distributed, with different locations having different environmental and demographic characteristics.
5	Population varies	Human population is not constant i.e. birth rate is not equal to death rate.

Table 2. Human Components assumptions.

S/No	Assumptions	Meaning
1.	Age structure	The human components have an age structure, with different age groups having different susceptibility to infection.
2.	Genetic variation	The human population is heterozygous, with different genotypes having different levels of susceptibility to infection.
3.	Heterogeneous mixing	Human mix non-randomly, with different age groups and genotypes having different contact rates.
4.	Spatial structure	The human population is spatially distributed, with different locations having different environmental and demographic characteristics.
5	Population varies	Human population is not constant i.e. birth rate is not equal to death rate.

Table 3. Blackfly assumptions.

S/No	Assumptions	Meaning
1.	Age structure	The blackfly population have an age structure, with different age group having different feeding bahaviours and infections rate
2.	Genetic variation	The blackfly population is also heterozygous, with different genotypes having different levels of susceptibility to infection
3.	Heterogeneous feeding	Blackfly feed non-randomly, with different age groups and genotypes having different feeding preferences
4.	Spatial distribution	The blackfly population is spatially distributed, with different locations having different environmental and demographic characteristics

Table 3. Blackfly assumptions.

S/No	Assumptions	Meaning
1.	Age structure	The blackfly population have an age structure, with different age group having different feeding bahaviours and infections rate
2.	Genetic variation	The blackfly population is also heterozygous, with different genotypes having different levels of susceptibility to infection
3.	Heterogeneous feeding	Blackfly feed non-randomly, with different age groups and genotypes having different feeding preferences
4.	Spatial distribution	The blackfly population is spatially distributed, with different locations having different environmental and demographic characteristics

Table 4. Disease transmission assumptions.

S/No	Assumptions	Meaning
1.	Frequency transmission dependent	The rate of disease transmission depends on the frequency of contact between humans and blackflies.
2.	Genotic-specific transmission	The rate of disease transmission depends on the genotype of the human and blackfly.
3.	Age-specific transmission	The rate of disease transmission depends on the age of human and blackfly.
4.	Spatially explicit transmission	The rate of disease transmission depends on the spatial location of human and blackfly.

Table 4. Disease transmission assumptions.

S/No	Assumptions	Meaning
1.	Frequency transmission dependent	The rate of disease transmission depends on the frequency of contact between humans and blackflies.
2.	Genotic-specific transmission	The rate of disease transmission depends on the genotype of the human and blackfly.
3.	Age-specific transmission	The rate of disease transmission depends on the age of human and blackfly.
4.	Spatially explicit transmission	The rate of disease transmission depends on the spatial location of human and blackfly.

1.5. Human Population

The humans population are categorize into susceptible compliant ($S_c$),and non-compliant ($S_n$) respectively. The sensitization, education about onchocerciasis system; treatment and how to prevent themselves is given to non-compliant human beings. It is assume that the person who comply with ivermectin and protect him/herself will move to susceptible compliant human individuals at the rate ${\phi }_1$while exposed human individuals become educated at rate ${\phi }_2$and infected non-compliant individual become educated at rate ${\phi }_3$. The ${\phi }_2$and${\phi }_3$will later move to their respective compartment (see the definition below) (Asha & Nyimvua, 2020):

${\mu }_h{N}_h$: the rate of increase of the susceptible human individuals from new born

${\lambda }_h$: the rate at which susceptible compliant and non-compliant individuals become infected after being bitten by the blackfly.

${\theta }_b$: the biting rate of blackfly.

$L_b$: larvae of susceptible blackfly

$I_b$: infectious blackfly

$S_h$: Susceptible human individuals

$E_h$: Exposed human individuals

$I_h$: Infectious humans; these categories are carrying the Onchocercavolvulus and can transmit the disease to others.

$R_h$: Recovered humans, these groups have recovered from the blackfly disease and have developed some immunity to the disease; however, the immunity may wane over time. See the descriptions in the following Tables:

Table 5. Description of the model Variables.

Variable	Meaning
$S_n$	Population of a susceptible non-compliant humans
$S_c$	Population of a susceptible compliant humans
$E_n$	Population of exposed non-compliant human individuals
$E_c$	Population of exposed compliant human individuals
$I_n$	Non-compliant Infectious humans; these categories are carrying the Onchocercavolvulus and can transmit the disease to others
$I_c$	Compliant Infectious humans; these categories are carrying the Onchocercavolvulus and can transmit the disease to others
$R_h$	Recovered compliant and non-compliant humans, these groups have recovered from the blackfly disease
$E_g$	Eggs laid by the blackfly
$L_b$	larvae of susceptible blackfly
$I_b$	infectious blackfly

Table 5. Description of the model Variables.

Variable	Meaning
$S_n$	Population of a susceptible non-compliant humans
$S_c$	Population of a susceptible compliant humans
$E_n$	Population of exposed non-compliant human individuals
$E_c$	Population of exposed compliant human individuals
$I_n$	Non-compliant Infectious humans; these categories are carrying the Onchocercavolvulus and can transmit the disease to others
$I_c$	Compliant Infectious humans; these categories are carrying the Onchocercavolvulus and can transmit the disease to others
$R_h$	Recovered compliant and non-compliant humans, these groups have recovered from the blackfly disease
$E_g$	Eggs laid by the blackfly
$L_b$	larvae of susceptible blackfly
$I_b$	infectious blackfly

Table 6. Description of the model Parameters in the human population.

parameters	Meaning
${\phi }_1$	The rate of the movement of a compliant individuals
${\phi }_2$	The rate at which exposed individuals become educated
${\phi }_3$	The rate at which infected non-compliant humans become educated
${\mu }_h{N}_h$	The rate of increase of a susceptible individuals from a new born (recruitment rate)
${\beta }_h$	Probability that contact between the infected blackfly and susceptible human causes infections to human
$1-\delta$	Reduction rate of transmission of a compliant individual
$0<\delta <1$	The efficacy of public health education
$\delta$	The rate of movement of exposed compliant and non-compliant after developing a symptoms of onchocerceciaisis
$\gamma$	Rate of treatment of infected humans
$\rho$	Rate of recovery of the infectious humans
${\nu }_1$	Modification parameters that reduce transmission on infectious humans
${\nu }_2$	Modification parameters that reduce transmission on infectious treated humans
$\varphi$	The rate at which recovered individuals lose their immunity and become susceptible again

Table 6. Description of the model Parameters in the human population.

parameters	Meaning
${\phi }_1$	The rate of the movement of a compliant individuals
${\phi }_2$	The rate at which exposed individuals become educated
${\phi }_3$	The rate at which infected non-compliant humans become educated
${\mu }_h{N}_h$	The rate of increase of a susceptible individuals from a new born (recruitment rate)
${\beta }_h$	Probability that contact between the infected blackfly and susceptible human causes infections to human
$1-\delta$	Reduction rate of transmission of a compliant individual
$0<\delta <1$	The efficacy of public health education
$\delta$	The rate of movement of exposed compliant and non-compliant after developing a symptoms of onchocerceciaisis
$\gamma$	Rate of treatment of infected humans
$\rho$	Rate of recovery of the infectious humans
${\nu }_1$	Modification parameters that reduce transmission on infectious humans
${\nu }_2$	Modification parameters that reduce transmission on infectious treated humans
$\varphi$	The rate at which recovered individuals lose their immunity and become susceptible again

Table 7. Description of the model Parameters in the blackfly population.

parameters	Meaning
$K$	Carrying capacity
${\mu }_l$	the rate at which larvae population decreases
$\eta$	Death of larvae due to larviciding
$\psi$	the rate of surviving mature larvae that will move to susceptible blackfly
$1-\epsilon$	The rate of the reduction or decrease of surviving blackfly
${\lambda }_b$	The rate at which susceptible blackfly get infectious when taking blood meal
${\beta }_b$	Probability that contact between susceptible blackflies and infectious human cause infections to blackfly
$\alpha$	The rate at which susceptible blackfly die
${\mu }_b$	Rate at which susceptible and infected blackfly die naturally
$\alpha$	The rate which susceptible and infected blackfly are trap

Table 7. Description of the model Parameters in the blackfly population.

parameters	Meaning
$K$	Carrying capacity
${\mu }_l$	the rate at which larvae population decreases
$\eta$	Death of larvae due to larviciding
$\psi$	the rate of surviving mature larvae that will move to susceptible blackfly
$1-\epsilon$	The rate of the reduction or decrease of surviving blackfly
${\lambda }_b$	The rate at which susceptible blackfly get infectious when taking blood meal
${\beta }_b$	Probability that contact between susceptible blackflies and infectious human cause infections to blackfly
$\alpha$	The rate at which susceptible blackfly die
${\mu }_b$	Rate at which susceptible and infected blackfly die naturally
$\alpha$	The rate which susceptible and infected blackfly are trap

1.6. Model Formulation

Blackflies and human hosts are the populations taken into account in the model. The human host population is classified as either compliant (those who are aware of the disease and take ivermectin) or non-compliant (those who have a negative attitude toward the mass therapy, ivermectin, and those who are ineligible to the drug) (Asha & Nyimvua, 2020; Torsen et al., 2025).The model will have twelve compartments at any given time t from these two groups of human hosts: susceptible human individuals who are at risk of contracting infections, represented by S_i; exposed human individuals who are infected but not infectious, represented by E_i; infected human individuals who are infectious and capable of spreading the disease, represented by I_i; and human individuals undergoing microfilariae treatment, represented by T; R_i (R, c,n), where c and n are constant, R represents a recovered human individual who has been treated from the microfilariae and is still with the adult worms; in this case, c stands for compliant individuals and n for non-compliant individuals (Asha & Nyimvua, 2020; Torsen et al., 2025). Also, the adult female blackfly population is separated into susceptible and infected black flies, represented by S_b and I_b, respectively. The larvae of the black fly (vector) population begin with L_b. Since the sickness does not cause mortality, it is assumed that the human population is constant (birthrate = deathrate), while the population of blackflies varies.

Furthermore, it is believed that the number of susceptible persons in the human population increases at a rate of μ_hN_h from newborns. Ivermectin is administered to all people age five years and older, regardless of their disease status, with the exception of those who are ineligible, such as children under five, pregnant women, and seriously ill humans, even though it kills microfilariae and reduces fertility to adult female worms. Therefore, vulnerable non-compliant people are educated about onchocerciasis (symptoms, treatment, and how to avoid blackflies), and those who take ivermectin and take precautions against blackflies spread to susceptible compliant people at a rate of ${\phi }_1$. Furthermore, non-compliant individuals who are exposed or infected receive education at rates of ${\phi }_2$ and ${\phi }_3$, respectively, and relocate to their corresponding compartments. However, the biting rate of an infectious female blackflies is given by λ_h,. After this, susceptible non-compliant and compliant individuals may become infected and these group of people will move to exposed classes. The human force of infection, represented by the parameter λ_h, is (Asha & Nyimvua, 2020; Torsen et al., 2025):

$${\lambda }_h={\displaystyle \frac{\theta {\beta }_h{I}_b}{N_h}}$$

(1)

where β_h is the likelihood that interaction between susceptible humans and infected blackflies will result in human infection, and θ is the biting rate of blackflies. Additionally, onchocerciasis can be prevented from spreading within the most affected communities through public health education. By preventing blackfly bites, the disease can be avoided. This may involve educating individuals on the use of insect repellent and appropriate clothes. Hence, the rate of transmission is lowered by (1−σ), where 0<σ<1 represents the effectiveness of public health education, and it is believed that vulnerable complying individuals are aware of the disease and know how to protect themselves from female blackflies. After a given period of time, exposed compliant and non- compliant individuals develop symptoms of onchocerciasis and move to infectious compartments at the rate δ. Since ivermectin kills microfilariae and not the macrofilariae, the infectious compliant human will be treated for microfilariae at a rate γ and move to the treated compartment. The treated human will recover from onchocerciasis after being cleared from the microfilariae parasites due to their natural death at a rate ρ and move to susceptible compliant and non-compliant human individuals.

On the other hand, it is assumed that the growth of larvae depends on hatched eggs laid by susceptible and infected blackflies taking into account the carrying capacity of the breeding site given by $\pi \left(S_b+L_b\right)\left(1-{\displaystyle \frac{L_b}{K}}\right)$, where π is the number of laid eggs and K is the carrying capacity depending on the space available, food and fresh air. The larvae population decreases by natural death at a rate μ_l and death due to larviciding at a rate η. The surviving larvae will mature at a rate ψ and move to susceptible adult female black flies.

But the maturation rate is reduced by a rate $\left(1-\epsilon \right)$, where $\epsilon$ is the efficacy of larvicide. Male blackflies are not considered as susceptible because they are not playing any role in the transmission of onchocerciasis. Susceptible female blackflies get infections when taking blood meal from infectious human individuals (compliant, non-compliant or treated and recovered human) at a rate λ_b. The parameter λ_b is the force of infection for blackflies given by:

$${\lambda }_b={\displaystyle \frac{\theta {\beta }_b\left(v_1{I}_c+I_n+v_{2}R\right)}{N_h}}$$

(2)

where β_b is the probability that contact between susceptible black flies and infectious human cause infections to blackflies,ν₁ and ν₂ are modification parameters which reduce transmission on infectious compliant and recovered humans respectively. Using trapping as a strategy of reducing the number of black flies, susceptible and infected blackflies are trapped at a rate α. The model assumes that the infected blackflies do not recover. The susceptible and infected female blackflies (vector) die naturally at a rate μ_b.

Furthermore, literature showed that a number of climatic variables such as temperature and rainfall play significant role in the transmission dynamic of vector born disease such as Onchocerca Volvulus. The significant role temperature plays on the dynamic of transmission of onchocerciasis are in different forms. For example, conducive temperatures generally cause the vector to feed and breed more frequently. On the other hand, higher temperatures make it at times difficult for the vector to breed. The high temperatures also destroy the Onchocerca eggs.

The effect of rainfall on onchocerciasis transmission dynamics is directly tied to availability of Onchocerca breeding sites. The use of precipitation as an empirical of incidence, as well as its direct impact on Onchocerca abundance has been established (Kamaldeen and Abba, 2016). Although, total rainfall increases the Onchocerca availability and productivity of breeding sites; increasing rainfall increases the vector abundance. Yet, excessive rainfall can lead to wash out of its breeding sites. For this model, these two (2) climatic variables are incorporated in the model (see Figure 1).

The variables $T=T\left(t\right)$ and $R=R\left(t\right)$ denote temperature and rainfall at time t. We assumed each of these functions to be continuous, bounded and positive and a function of time. The temperature dependent eggs deposition rate is denoted by ${\alpha }_3\left(T\right)$ and the population of immature is decreased by maturation to adult at a temperature and rainfall dependent rate of ${\xi }_2\left(T,R\right)$, a natural death at a temperature dependent rate of${\mu }_2\left(T\right)$. The population of the immature Onchocerca Vulvolus population growth is limited by the carrying capacity $K_v={\alpha }_3\left(T\right)\left(1-{\displaystyle \frac{V_e}{K_v}}\right)$. The force of interaction is given by ${\lambda }_k$which is govern by the equation:

$${\lambda }_k={\displaystyle \frac{{\alpha }_3\left(T\right)+{\mu }_2\left(T\right)}{{\xi }_2\left(T,R\right)}}$$

(3)

The interaction between host, vector populations, temperature and rainfall is illustrated in Figure 1.

The equation governing the above equation is given below:

$${\displaystyle \frac{d{S}_n}{dt}}={\mu }_h{N}_h-\left({\lambda }_n+{\mu }_h+{\phi }_1\right)S_n$$

(4)

$${\displaystyle \frac{d{S}_c}{dt}}={\phi }_1{S}_n+\rho R-((1-\delta ){\lambda }_h+{\mu }_h)S_c$$

(5)

$${\displaystyle \frac{d{E}_n}{dt}}={\lambda }_h{S}_n-(\left(\delta +{\phi }_2+{\mu }_n\right)E_n$$

(6)

$${\displaystyle \frac{d{E}_c}{dt}}=\left(1-\delta \right){\lambda }_n{S}_c+{\phi }_3{E}_n-\left(\delta +{\mu }_h\right)E_c$$

(7)

$${\displaystyle \frac{d{I}_n}{dt}}=\delta E_n-\left({\phi }_3+{\mu }_h\right)I_n$$

(8)

$${\displaystyle \frac{d{I}_c}{dt}}=\delta E_c+{\phi }_3{I}_n-\left(\gamma +{\mu }_h\right)I_c$$

(9)

$${\displaystyle \frac{dT}{dt}}=\gamma I_c-\left(\omega +{\mu }_h\right)T$$

(10)

$${\displaystyle \frac{dR}{dt}}=\varpi T-\left(\rho +{\mu }_h\right)R$$

(11)

$$N_h=S_n+S_c+E_n+E_c+I_n+I_c+T+R$$

(12)

$$S_h=S_n+S_c$$

(13)

The equation governing the blackfly interaction with its host and the two climatic variables is given below:

The total population for the blackfly is given by:

$$N_b=L_b+S_b+I_b$$

(14)

$${\displaystyle \frac{d{E}_e}{dt}}=\prod -[\left(\psi +{\nu }_e\right)+\left(1-\epsilon \right)\psi +{\alpha }_1\left(T_m\right)]E_e$$

(15)

$${\displaystyle \frac{d{L}_b}{dt}}=\pi \left(L_b+S_b\right)\left(1-{\displaystyle \frac{L_b}{K}}\right)-\left(1-\epsilon \right)\psi -\left(\eta +{\mu }_l\right)-{\alpha }_2\left(R_f\right)$$

(16)

$${\displaystyle \frac{d{S}_b}{dt}}=\left(1-\epsilon \right)\psi L_b-{\lambda }_b{S}_b-\left(\alpha +{\mu }_b\right)I_b$$

(17)

$${\displaystyle \frac{d{I}_b}{dt}}={\lambda }_b{S}_b-\left(\alpha +{\mu }_b\right)I_b$$

(18)

The model has basic properties as described below:

The positivity of the solution, invariant region that is its boundedness as well as the existence and the stability of the equilibrium points of the model system are subjectively analyzed.

1.7. Positivity and Boundedness of the Solution

To show that the model is biologically meaningful and very well posed, we have to prove beyond reasonable doubt that the solution to all the state variables are positives at all level and all the time (Nyimvua, 2024; Huo & Qui, 2014; Huo et al., 2014). Consider the following Lemma:

1.8. Lemma

The solution for all state variables of the system equation 3.4 to 3.18 with positive initial values are all positive for all t >0.

Proof: Nyimvua (2024); Huo & Qiu; here we prove the positivity by contradiction as seen below:

Given the non-negative initial value conditions S_n(0)>0, S_n(t₀)=0, S`n(t₀) ≤ 0 S_c(t) > 0 E_n(t) >0, E_c(t) >0, I_n(t) >0, I_c(t) > 0, L_b(t) > 0, S_b(t) > 0, 0≤ t < t₀. However, from equation (3.4), we have:

$${\displaystyle \frac{d{S}_n\left(t_0\right)}{dt}}={\mu }_h{N}_h-\left({\lambda }_n+{\mu }_h+{\phi }_1\right)S_n\left(t_0\right)$$

(19)

This implies that

$${\displaystyle \frac{d{S}_n\left(t_0\right)}{dt}}={\mu }_h{N}_h>0,$$

(20)

This contradict our assumption and this means that

$S_n\left(t\right)>0$for t >0. Also, we assumed that there exist a slightly time t₁ such that conditions S_n(0)>0, S_n(t₁)=0, S`n(t₁) ≤ 0 S_c(t₁) > 0 E_n(t) >0, E_c(t) >0, I_n(t) >0, I_c(t) > 0, L_b(t) > 0, S_b(t) > 0, 0≤ t < t_1.

Also from equation (3.5), we apply the same scenario as follows:

$${\displaystyle \frac{d{S}_c\left(t_1\right)}{dt}}={\phi }_1{S}_n+\rho R-((1-\delta ){\lambda }_h+{\mu }_h)S_c\left(t_1\right)$$

(21)

which gives:

$${\displaystyle \frac{d{S}_c\left(t_1\right)}{dt}}={\phi }_1{S}_n+\rho R>0$$

(22)

This also contradicts with the assumption set forth; therefore it means that

$$S_n\left(t\right)>0$$

We applied the same process and shows that:

$$E_n\left(t\right)\ge 0,E_c\left(t\right)\ge 0,I_n\left(t\right)\ge 0,I_c\left(t\right)\ge 0,L_b\left(t\right)\ge 0,S_b\left(t\right)\ge 0,I_b\left(t\right)\ge 0$$

for all t >0. Hence, the solutions for all state variables of the model system are non-negative.

1.9. Invariant Region

Here the region in which the solution of the model has biological meaning is obtained is seen below. There are two populations hence, let the feasible region be given as $\mathsf{\Phi }={\mathsf{\Phi }}_h+{\mathsf{\Phi }}_b\in {\Re }_{+}^{8}X{\Re }_{+}^{3}where,{\mathsf{\Phi }}_h=\left\{S_n,S_c,E_n,E_c,I_n,I_c,T,R\in {\Re }_{+}{}^8:S_n+S_c+E_n+E_c+I_n+I_c+T+R=N_h\right\}$and

$${\mathsf{\Phi }}_b=\left\{L_b,S_b,I_b\in {\Re }_{+}^3{L}_b\le K,S_b+I_b\le N_b\le {\displaystyle \frac{\left(1-\epsilon \right)K}{\alpha +{\mu }_b}}\right\}.$$

(23)

To prove the feasible region we add the equations to give us the following:

$${\displaystyle \frac{d{\sout{N}}_h}{dt}}={\mu }_h{N}_h-{\mu }_h\left(S_n+S_c+E_n+E_c+I_n+I_c+T+R\right)$$

(24)

Which gives

$${\displaystyle \frac{d{\sout{N}}_h}{dt}}={\mu }_h{N}_h-{\mu }_h{N}_h=0$$

(25)

Equation (25) shows a constant population. If we consider the blackfly population, we have :

$${\displaystyle \frac{d{\sout{N}}_b}{dt}}=\left(1-\epsilon \right)K-\left(\alpha +{\mu }_b\right)\left(S_b+I_b\right)$$

(26)

Implying that

$${\displaystyle \frac{d{\sout{N}}_b}{dt}}=\left(1-\epsilon \right)K-\left(\alpha +{\mu }_b\right)N_b$$

(27)

Since the human population is assumed to be constant from equation (25) , this implies that:

$$S_n\le N_h,S_c\le N_h,E_n\le N_h,E_c\le N_h,I_n\le N_h,I_c\le N_h,T\le N_h,R\le N_h$$

Solving equation (3.27) yields

$$N_b\left(t\right)={\displaystyle \frac{\left(1-\epsilon \right)K}{\alpha +{\mu }_b}}+\left(N_b\left(0\right)-{\displaystyle \frac{\left(1-\epsilon \right)K}{\alpha +{\mu }_b}}\right)e^{-\left(\alpha +{\mu }_b\right)t}$$

(28)

We applied the differential inequally theorem in (28) and obtained :

$0\le N_b\left(t\right)\le {\displaystyle \frac{\left(1-\epsilon \right)K}{\alpha +{\mu }_B}}$, as t tends to infinity.

Hence, the solution on the two populations (human and blackfly) enters the invariant region

$$\begin{array}{l}{\mathsf{\Phi }}_h=\left\{\left(S_n,S_c,E_n,E_c,I_n,I_c,T,R\right)\in {\Re }^8{}_{+}+S_n+S_c+E_n+E_c+I_n+I_c+T+R=N_h\right\}and\\ {\mathsf{\Phi }}_b=\left\{\left(L_b,S_b,I_b\right)\in {\Re }^3{}_{+}:L_b\le K,N_b\le {\displaystyle \frac{\left(1-\epsilon \right)K}{\alpha +{\mu }_b}}\right\}.\end{array}$$

(29)

Hence, the region is bounded and it attracts all solutions in $\mathsf{\Phi }={\mathsf{\Phi }}_h\cup {\mathsf{\Phi }}_b$. We can therefore say that the model is biologically meaningful.

1.10. Disease Free Equilibrium

The disease free equilibrium point is a steady state solution of the model when there is no disease in both human host and blackfly vector populations. This is to say that at disease free equilibrium, $E_n=E_c=I_n=I_c=T=R=I_b=0$. Putting all infected classes equal to zero in the model system equations and solve, we have the following equilibrium points. We infer that trivial disease free equilibrium point is when blackflies population are completely extinct which is not realistic in a real situation as it’s not possible to eliminate all the blackflies population. This is given by the following:

$$E^0=\left(S^{*}{}_n,S^{*}{}_c,0,0,0,0,0,0,0,0\right)=\left({\displaystyle \frac{{\mu }_h{N}_h}{{\phi }_1+{\mu }_h}},{\displaystyle \frac{{\phi }_1{N}_h}{{\phi }_1+{\mu }_h}},0,0,0,0,0,0,0,0,0\right)$$

(30)

Equation (30) is not biologically meaningful. The biologically meaningful steady state is given below:

$$\begin{array}{l}{E}^{*}=\left(S^{*}{}_n,S^{*}{}_c,0,0,0,0,L_b,S_b^{*},0,0\right),where\\ S_n^{\ast }={\displaystyle \frac{{\mu }_{h1}{N}_h}{{\phi }_1+{\mu }_h}},S_c^{\ast }={\displaystyle \frac{{\phi }_1{N}_h}{{\phi }_1+{\mu }_h}},L_b^{\ast }={\displaystyle \frac{K\left(\left(\alpha +{\mu }_b\right)\left(1-\epsilon \right)\psi +\eta +{\mu }_l\right)\left(\frac{\left(1-\epsilon \right)\pi \psi }{\left(\alpha +{\mu }_b\right)\left(\left(1-\epsilon \right)\psi +\eta +{\mu }_l\right)}\right)}{\pi \left(1-\epsilon \right)\psi }},\\ S_b^{\ast }={\displaystyle \frac{K\left(\left(1-\epsilon \right)\psi +\eta +{\mu }_l\right)\left(\alpha +{\mu }_b\right)}{\pi }}\end{array}$$

(31)

1.11. Reproduction Number $\left(R_0\right)$

The reproduction number is very significant in modeling epidemiological phenomenon. A very important threshold quantity is the basic reproduction number, sometimes called the basic reproductive number or basic reproductive ratio (Heffernan et al., 2005), which is usually denoted by R₀. The epidemiological definition of R₀ is the average number of secondary cases produced by one infected individual introduced into a population of susceptible individuals, where an infected individual has acquired the disease, and susceptible individuals are healthy but can acquire the disease. In reality, the value of R₀ for a specific disease depends on many variables, such as location and density of population. The threshold (R₀) has three interpretations: for instance, if R₀ < 1, then the number of infectious individuals decreases monotonically to 0. Meaning that the disease is going to die out or be eradicated. If R₀ > 1, the disease will increase or the incidence rate will increase or there will be a spread of the disease and if R₀ = 1, the disease will steadily increase but at a slower rate . In this work, we used the next generation matrix (NGM) to compute the R₀ as seen below:First, the Jacobian matrix is given by the system of equation represented in equation 32:

$$J=\left[\begin{array}{cccccccc}{\displaystyle \frac{\partial S_n^{\prime }}{\partial S_n}}& {\displaystyle \frac{\partial S_n^{\prime }}{\partial S_c}}& {\displaystyle \frac{\partial S_n^{\prime }}{\partial E_n}}& {\displaystyle \frac{\partial S_n^{\prime }}{\partial E_c}}& {\displaystyle \frac{\partial S_n^{\prime }}{\partial I_n}}& {\displaystyle \frac{\partial S_n^{\prime }}{\partial I_c}}& {\displaystyle \frac{\partial S_n^{\prime }}{\partial T}}& {\displaystyle \frac{\partial S_n^{\prime }}{\partial R}}\\ {\displaystyle \frac{\partial S_c^{\prime }}{\partial S_n}}& {\displaystyle \frac{\partial S_c^{\prime }}{\partial S_c}}& {\displaystyle \frac{\partial S_c^{\prime }}{\partial E_n}}& {\displaystyle \frac{\partial S_c^{\prime }}{\partial E_c}}& {\displaystyle \frac{\partial S_c^{\prime }}{\partial I_n}}& {\displaystyle \frac{\partial S_c^{\prime }}{\partial I_c}}& {\displaystyle \frac{\partial S_c^{\prime }}{\partial T}}& {\displaystyle \frac{\partial S_c^{\prime }}{\partial R}}\\ {\displaystyle \frac{\partial E_n^{\prime }}{\partial S_n}}& {\displaystyle \frac{\partial E_n^{\prime }}{\partial S_c}}& {\displaystyle \frac{\partial E_n^{\prime }}{\partial E_n}}& {\displaystyle \frac{\partial E_n^{\prime }}{\partial E_c}}& {\displaystyle \frac{\partial E_n^{\prime }}{\partial I_n}}& {\displaystyle \frac{\partial E_n^{\prime }}{\partial I_c}}& {\displaystyle \frac{\partial E_n^{\prime }}{\partial T}}& {\displaystyle \frac{\partial E_n^{\prime }}{\partial R}}\\ {\displaystyle \frac{\partial E_c^{\prime }}{\partial S_n}}& {\displaystyle \frac{\partial E_c^{\prime }}{\partial S_c}}& {\displaystyle \frac{\partial E_c^{\prime }}{\partial E_n}}& {\displaystyle \frac{\partial E_c^{\prime }}{\partial E_c}}& {\displaystyle \frac{\partial E_c^{\prime }}{\partial I_n}}& {\displaystyle \frac{\partial E_c^{\prime }}{\partial I_c}}& {\displaystyle \frac{\partial E_c^{\prime }}{\partial T}}& {\displaystyle \frac{\partial E_c^{\prime }}{\partial R}}\\ {\displaystyle \frac{\partial I_n^{\prime }}{\partial S_n}}& {\displaystyle \frac{\partial I_n^{\prime }}{\partial S_c}}& {\displaystyle \frac{\partial I_n^{\prime }}{\partial E_n}}& {\displaystyle \frac{\partial I_n^{\prime }}{\partial E_c}}& {\displaystyle \frac{\partial I_n^{\prime }}{\partial I_n}}& {\displaystyle \frac{\partial I_n^{\prime }}{\partial I_c}}& {\displaystyle \frac{\partial I_n^{\prime }}{\partial T}}& {\displaystyle \frac{\partial I_n^{\prime }}{\partial R}}\\ {\displaystyle \frac{\partial I_c^{\prime }}{\partial S_c}}& {\displaystyle \frac{\partial I_c^{\prime }}{\partial S_c}}& {\displaystyle \frac{\partial I_c^{\prime }}{\partial E_n}}& {\displaystyle \frac{\partial I_c^{\prime }}{\partial E_c}}& {\displaystyle \frac{\partial I_c^{\prime }}{\partial I_n}}& {\displaystyle \frac{\partial I_c^{\prime }}{\partial I_c}}& {\displaystyle \frac{\partial I_c^{\prime }}{\partial T}}& {\displaystyle \frac{\partial I_c^{\prime }}{\partial R}}\\ {\displaystyle \frac{\partial T^{\prime }}{\partial S_n}}& {\displaystyle \frac{\partial T^{\prime }}{\partial S_c}}& {\displaystyle \frac{\partial T^{\prime }}{\partial E_n}}& {\displaystyle \frac{\partial T^{\prime }}{\partial E_c}}& {\displaystyle \frac{\partial T^{\prime }}{\partial I_n}}& {\displaystyle \frac{\partial T^{\prime }}{\partial I_c}}& {\displaystyle \frac{\partial T^{\prime }}{\partial T}}& {\displaystyle \frac{\partial T^{\prime }}{\partial R}}\\ {\displaystyle \frac{\partial R^{\prime }}{\partial S_n}}& {\displaystyle \frac{\partial R^{\prime }}{\partial S_c}}& {\displaystyle \frac{\partial R^{\prime }}{\partial E_n}}& {\displaystyle \frac{\partial R^{\prime }}{\partial E_c}}& {\displaystyle \frac{\partial R^{\prime }}{\partial I_n}}& {\displaystyle \frac{\partial R^{\prime }}{\partial I_c}}& {\displaystyle \frac{\partial R^{\prime }}{\partial T}}& {\displaystyle \frac{\partial R^{\prime }}{\partial R}}\end{array}\right]$$

(32)

$$J=\left[\begin{array}{cccccccc}-\left({\lambda }_n+{\mu }_n+{\phi }_1\right)& 0& -{\lambda }_n{S}_n^{\ast }& 0& 0& 0& 0& 0\\ \phi & -\left(\left(1-\delta \right){\lambda }_n+{\mu }_n\right)& 0& -\left(1-\delta \right){\lambda }_n{S}_c^{\ast }& 0& 0& 0& \rho \\ {\lambda }_n{S}_n^{\ast }& 0& -\left(\delta +{\phi }_2+{\mu }_n\right)& 0& 0& 0& 0& 0\\ 0& \left(1-\delta \right){\lambda }_n{S}_c^{\ast }& {\phi }_3& -\left(\delta +{\mu }_n\right)& 0& 0& 0& 0\\ 0& 0& \delta & 0& -\left({\phi }_3+{\mu }_n\right)& 0& 0& 0\\ 0& 0& 0& \delta & {\phi }_3& -\left(\gamma +{\mu }_n\right)& 0& 0\\ 0& 0& 0& 0& 0& \gamma & -\left(\omega +{\mu }_b\right)& 0\\ 0& 0& 0& 0& 0& 0& \omega & -\left(\psi +{\mu }_n\right)\end{array}\right]$$

(33)

The Jacobian matrix above was obtained by differentiating the system of the above equations. The F matrix is given in equation (34):

$$F=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& {\displaystyle \frac{\theta S^{*}{}_n\beta _h}{N_h}}\\ 0& 0& 0& 0& 0& {\displaystyle \frac{\left(1-\sigma \right)\theta S^{*}{}_c\beta _h}{N_h}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{S^{*}{}_b\beta _b}{N_h}}& {\displaystyle \frac{\theta S^{*}{}_b\beta _b{\nu }_1}{N_h}}& {\displaystyle \frac{\theta S^{*}{}_b\beta _b{\nu }_2}{N_h}}& 0\end{array}\right]$$

(34)

The matrix V of the system of the disease if given in equation (35) below

$$V=\left[\begin{array}{cccccc}\delta +{\phi }_2+{\mu }_h& 0& 0& 0& 0& 0\\ -{\phi }_2& \delta +{\mu }_h& 0& 0& 0& 0\\ -\delta & 0& {\phi }_3+{\mu }_h& 0& 0& 0\\ 0& -\delta & -{\phi }_3& \gamma +\mu & 0& 0\\ 0& 0& 0& -\gamma & \rho +{\mu }_h& 0\\ 0& 0& 0& 0& 0& \alpha +{\mu }_b\end{array}\right]$$

(35)

The inverse of the matrix V is given in equation ( 36) below

$$V^{-1}=\left[\begin{array}{cccccc}{\displaystyle \frac{1}{\delta +{\phi }_2+{\mu }_i}}& 0& 0& 0& 0& 0\\ {\displaystyle \frac{{\phi }_2}{\left(\delta +{\phi }_2+{\mu }_h\right)\left(\delta +{\mu }_h\right)}}& {\displaystyle \frac{1}{\left(\delta +{\mu }_h\right)}}& 0& 0& 0& 0\\ {\displaystyle \frac{\delta }{\left(\delta +{\phi }_1+{\mu }_h\right)\left({\phi }_3+{\mu }_h\right)}}& 0& {\displaystyle \frac{1}{\left({\phi }_3+{\mu }_h\right)}}& 0& 0& 0\\ {\displaystyle \frac{\delta \left(\delta +{\phi }_3\right)}{\left(\delta +{\phi }_2+{\mu }_h\right)\left({\phi }_3+{\mu }_h\right)\left(\gamma +{\mu }_h\right)}}& {\displaystyle \frac{\delta }{\left(\sigma +{\mu }_h\right)\left({}^{\gamma }+{\mu }_h\right)}}& {\displaystyle \frac{{\phi }_3}{\left({\phi }_3+{\mu }_h\right)\left({}^{\gamma }+{\mu }_h\right)}}& {\displaystyle \frac{1}{\gamma +{\mu }_b}}& 0& 0\\ {\displaystyle \frac{\delta \gamma \left(\delta +{\phi }_3\right)}{\left(\delta +{\phi }_2+{\mu }_h\right)\left({\phi }_3+{\mu }_h\right)\left(\rho +{\mu }_h\right)}}& {\displaystyle \frac{\delta \gamma }{\left(\delta +{\mu }_h\right)\left({}^{\gamma }+{\mu }_h\right)\left(\rho +{\mu }_h\right)}}& {\displaystyle \frac{\gamma {\phi }_3}{\left({\phi }_3+{\mu }_h\right)\left({}^{\gamma }+{\mu }_h\right)\left({}^{\rho }+{\mu }_h\right)}}& {\displaystyle \frac{\gamma }{\left({}^{\gamma }+{\mu }_b\right)\left({}^{\rho }+{\mu }_h\right)}}& {\displaystyle \frac{1}{\left({}^{\rho }+{\mu }_h\right)}}& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{1}{\left(\rho +{\mu }_b\right)}}\end{array}\right]$$

(36)

The product V^-1 and V is given by equation (37) below

$$F{V}^{-1}=\left[\begin{array}{l}\begin{array}{cccccc}0& 0& 0& 0& 0& {\displaystyle \frac{\theta S_n^{\ast }{\beta }_h}{N_h\left(\alpha +{\mu }_b\right)}}\\ 0& 0& 0& 0& 0& {\displaystyle \frac{\left(1-\sigma \right)\theta S_c^{\ast }{\beta }_h}{N_h\left(\alpha +{\mu }_b\right)}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ {\displaystyle \frac{S_b^{\ast }\delta {\beta }_b\left({\phi }_3+{\mu }_h\right)\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_h\right)+\theta {\nu }_1\left(\delta +{\phi }_3\right)\left(\rho +{\mu }_h\right)+\theta {\nu }_2\gamma \left(\delta +{\phi }_3\right)}{N_h\left(\delta +{\phi }_2+{\mu }_h\right)\left({\phi }_3+{\mu }_h\right)\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_h\right)}}& \begin{array}{l}{\displaystyle \frac{S_b^{\ast }{\beta }_b\theta \sigma \left({\nu }_1\left(\rho +{\mu }_b\right)+{\nu }_2\gamma \right)}{N_h\left(\delta +{\mu }_h\right)\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_h\right)}},\end{array}& & & & 0\end{array}\\ {\displaystyle \frac{S_b^{\ast }{\beta }_b\left(\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_b\right)+\theta {\nu }_1{\phi }_3\left(\rho +\mu \right)+\theta {\nu }_1{\phi }_3\right),}{N_h\left({\phi }_3+{\mu }_h\right)\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_h\right)}}{\displaystyle \frac{S_b^{\ast }{\beta }_b\theta \left({\nu }_1\left(\rho +{\mu }_b\right)+{\nu }_2\gamma \right)}{N_h\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_h\right)}},\\ {\displaystyle \frac{S_b^{\ast }{\beta }_b\theta {\nu }_2}{N_h\left(\rho +{\mu }_h\right)}}\end{array}\right]$$

(37)

From the analysis, let A be the Onchocerca to human entry denoted by the given by (38)

$$A={\displaystyle \frac{\theta S_n^{\ast }{\beta }_h}{N_h\left(\alpha +{\mu }_b\right)}}$$

(38)

and let B, be the human entry of FV^-1 given by equation (39)

$$B={\displaystyle \frac{S_b^{\ast }{\beta }_b\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_h\right)+\theta {\nu }_1{\phi }_2\left(\gamma +{\mu }_b\right)+\theta {\nu }_2\gamma {\phi }_2}{N_h\left({\phi }_3+{\mu }_h\right)\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_h\right)}}$$

(39)

The reproduction number R₀ is obtained as seen below:

$$R_0=\sqrt{{\displaystyle \frac{\theta S_n^{\ast }{\beta }_h{S}_b^{\ast }\left(\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_h\right)+\theta {\nu }_1{\phi }_3\left(\rho +{\mu }_h\right)+\theta {\nu }_2\gamma {\phi }_3\right)}{N_h^2\left({\phi }_3+{\mu }_h\right)\left(\gamma +{\mu }_b\right)\left(\rho +{\mu }_h\right)\left(\alpha +{\mu }_b\right)}}}$$

(40)

The parameters values are seen in Table 8

3. Results

The data in Table 8 above were applied to equation (4) and analyzed using R software. The reproduction number (R₀) obtained is 0.020. The annual prediction was done based on the value of R₀ obtained using the same R software and the result is depicted in Table 9.

The result of Table 9 above were set randomly plot the graph for five (6) years that from 2026 to 2030 inclusive. The year 2030 is the next target year for the elimination of the disease as set by WHO and the government of Nigeria.

The basic reproduction number ($R_0$) was computed using the Next Generation Matrix and was found to be 0.0202. Since $R_0$= 0.0202 < 1, the model predicts that the disease will not persist but die out in the long run as new incidence cases are perceivably weak. The predicted values shown in Table 9 revealed the same results as the incidence and prevalence cases are also not high. However, if interventions are not met and ivermectin distributions are slowed down, the disease is going to rise again due to lack of control measures. The prevalence prediction rate is plotted Figure 1. The infected non-compliant, the infected compliant and the treated values were plotted and depicted in Figure 2.

The figure revealed the same result as the reproduction number also revealed. This indicates that the model fit well and by this results the disease will eventually die out if all the necessary measures to eliminate the disease are put in place.

Figure 2 represents the human infections dynamics: the red curve is for the infected non-compliant, the blue is for the infected compliant while the green is for the treated human individuals. The results revealed that the infected non-compliant has more infections at the onset on their non-compliant to the drug or any other means of protecting themselves. However, as they become aware through the awareness campaign and resort to taking drugs and protecting themselves. By so doing, they drastically showed great reduction in infections and are likely going to be free from the disease. Infected and treated compliant has shown drastic health recovery and it is an indication that the disease might be eliminated at the set target year of elimination. The human prevalence and the vector infections rate graph also revealed exponential decay same as the results as in figure 1 and 2 (see Figure 3)

Figure 3 above shows the human and vector infections rate. The blue dash line is the infections rate of vector and it was scaled up to be an independent variable. The human infections rate is shown in red. This also revealed the exponential decay indicating the drastic decay that shows dying rate of the disease in the long run. The impact interventions scenario is depicted in Figure 4.

The intervention base line is indicated by black while education and treated are shown in blue and green colours respectively. The combined intervention is shown in red. The education given to people on the need to take drug, protect themselves from the bites of the blackfly plays vital role as well as the treatments. These will greatly aid in eradicating or eliminating the disease. However, if these variables are combined and apply simultaneously at the same time, the elimination will be faster than doing it separately.

Discussion of Results

One of the most important and essential threshold numbers in mathematical model is the basic reproduction number (R₀). In the absence of any control measures, it is defined as the average number of secondary infections caused by a single infectious individual introduced into a community that is completely susceptible.

The calculated value of R₀ = 0.0202 obtained for Onchocerciasis transmission dynamics in this study has a profound significant. Since R₀ < 1, the interpretation is clear: it shows that the chain of transmission is self-limiting and mathematically assured to decay over time since each infectious individual typically produces fewer than one new infection. This is referred to as the disease free equilibrium (DFE), and it is locally asymptotically stable when R₀ < 1, suggesting that minor disturbances or the entrance of the disease into the population won't result in long lasting outbreaks. This indicates that, if present intervention efforts, particularly ivermectin Mass Drug Administration (MDA) public health campaigns and treatement are strictly maintained. It also shows that the onchocerciasis burden in the study zone is headed toward elimination.

It also revealed that the disease cannot sustain itself in the population since R₀ is less than the critical threshold of unity. With each new generation of transmission, the force of infection decreases and the epidemic curve sharply declines. This stands in contrast to a situation where R₀ > 1, which would indicate a stable endemic equilibrium that actively opposes eradication efforts and exponential development in the case load. The result is in tandem with the study of Torsen et al., (2025) whose study revealed that the disease can be eliminated with concerted efforts in public education since the R₀ < 1 . Hethcote ( 2000) results also revealed the same thing.

Figure 1 shows the five-year forecast of the onchocerciasis prevalence rate (%) plotted against time in years. From an initial prevalence of roughly 10% at time zero, the graph shows a fast monotonic fall that asymptotically converges to near-zero prevalence by year five. This decay is typical of a sub critical epidemic system, in which the degree of deviation from the threshold R₀ = 1 controls the pace of decline. The illness is markedly downward and the epidemic force is weak since R₀ = 0.0202 is significantly below unity. For instance, the nonlinear nature of epidemic extinction processes is seen in the prevalence curve initial gradient and subsequent gradual flattening. When active infections are treated and the transmission chain is broken, the disease rapidly diminishes in the early stages.

Figure 2 revealed that at the beginning of the simulation period, the infected non-compliant compartment (red line) shows the biggest peak, reaching over 2.5 million individuals before gradually declining over time. This substantial initial burden among non-compliant individuals highlights how important community adherence to ivermectin MDA is in reducing transmission. Even when overall coverage rates seem adequate, non-compliance with treatment caused by cultural beliefs, fear of negative medication effects, ignorance, or practical difficulties constitutes a significant epidemiological reservoir that perpetuates transmission. The model-predicted prevalence curve is shown in Figure 4.3. This is to determine how best the mathematical model fit the data. The model structure and parameter estimates are validated by the close correspondence between the two progressions, which both follow a quick initial fall followed by an asymptotic approach to near zero.

The current trajectory (progression) of disease decrease under current intervention levels is represented by the baseline scenario (blackline line). It acts as a benchmark for evaluating the added value of improved interventions. It slower pace of reduction in comparison to all other scenarios highlights the fact that, even though the current situation is improving, it will not be enough to accomplish elimination within the WHO's 2030 target timeframe without augmentation (see Figure 3). We see that the combined intervention scenario (blue line) that is simultaneous scale-up of ivermectin MDA, health education and treatment attains the fastest and most complete elimination, reaching near-zero prevalence well before year five. This finding is consistent with the synergistic principle in public health: no single intervention is maximally effective in isolation and the biggest gains come from integrated, multi-pronged approaches that address the biological (treatment) and the behavioral (education and compliance) dimensions of disease transmission. The results is also in line with the finding of Nyimvua et al.,(2020) who reported that effective elimination and control of onchocerciasis can best be done when all the four strategies they applied in their studies can be implemented simultaneously.

Conclusion

This work created a stochastic, individual-based mathematical model that incorporates human compliance behavior, blackfly vector stages, and climate influences to capture the diverse transmission dynamics of onchocerciasis in northeastern Nigeria. The model produced a basic reproduction number of R₀ = 0.0202, showing that the disease is on a trajectory toward elimination under present intervention efforts. It was demonstrated to be well-posed, with confirmed positivity and boundedness of solutions. Based on the simulations, infection trends and prevalence would decline through 2030. When compared to single-strategy methods, the combination of ivermectin distribution and public health education produced the fastest elimination outcomes. To achieve Nigeria's 2030 elimination goal, we recommend that health authorities should maintain and step up integrated interventions simultaneously increasing ivermectin mass drug administration, community education, and vector control instead of depending on any one tactic.