38 International Journal of Development Mathematics Vol 1 Issue 3, Page No.

038 - 059

INTERNATIONAL
INTERNATIONALJOURNAL
JOURNALOF
OFDEVELOPMENT
DEVELOPMENTMATHEMATICS
MATHEMATICS

ISSN: 3026-8656 (Print) | 3026-8699 (Online)

journal homepage: [Link]

Modeling the Impacts of Diagnosis and Treatment of Snakebite Envenoming Victims with
Anti-Snake Venom in a Community
Joseph S. Ayubaa*, Musa Abdullahib and Patience B. Iliyaa
aDepartment of Mathematics and Statistics, Federal University of Kashere, Gombe State, Nigeria.
bDepartment of Mathematics, Modibbo Adama University, Yola, Adamawa State, Nigeria.

ABSTRACT
ARTICLE INFO Snakebite envenoming (SBE) is one of the major causes of deaths and disabilities in tropical
and subtropical countries especially in the continents of Asia and Africa. Report revealed that
Article history: on average, 5-8 and about 3,000 cases of SBE are recorded daily and yearly respectively in
Received 17 March 2024 northeast Nigeria which population are largely agriculturists. Thus, snakebite (SB) led to
Received in revised form 05 June 2024 economy setback and slows down productivity among the workforce as a result of disabilities
Accepted 29 July 2024 induced by SBE. Consequently, the revenue generation even including income facing
declined. This paper proposes population-based mathematical model for studying the dynamic
impacts of diagnosis and treatment of snakebite envenoming victims with anti-snake venom
Keywords:
in a community. The model is rigorously analyzed using relevant techniques and theorems in
Anti-snake venom, Mathematical
model, Nneglected tropical diseases dynamical systems. The model was validated using monthly reported data on snakebite
Snakebite envenoming, Stability envenoming from 2020–2022, collected from snakebite treatment and research hospital
analysis (STRH) Kaltungo in northeast Nigeria. The qualitative study revealed that the model has
snakebite-envenoming-free and snakebite-envenoming-endemic equilibrium points. It is
shown via the Lyapunov function in conjunction with the Lassale invariance principle that
MSC 2020 Subject classification: these two equilibrium points are globally asymptotically stable. Furthermore, numerical
92-10 92D30 34D20 34D23 65L07 simulations of the model were carried out to investigate the impact of the diagnosis of
snakebite envenoming victims and the treatment of the same on the burden of snakebite
envenoming. The results of the simulations revealed that appropriate diagnosis and treatment
of snakebite envenoming victims with anti-snake venom could significantly reduce the
number of snakebite-induced deaths and disabilities. This work can serve as a framework for
making policy and decision in addressing the issues of SBE in the public healthcares in the
study region and beyond.

1. Introduction
Snakebite (SB) is an injury caused by a bite from a snake, resulting in puncture wounds inflicted by the fangs,
while snakebite envenoming (SBE) is a potentially life-threatening disease which occurs as a result of an injection of
a complex mixture of different toxins following the bite of a venomous snake (Warrell, 2010 and Williams, et al,.
2019). The severity of snakebite depends on some factors, such as the kind of snake, the part of the body bitten, the
quantity of venom injected, and the general health of the person bitten (Abdullahi, Habib & Hussaini, 2021). Generally,
there are more than 3000 species of snakes in the world, and they live in terrestrial and aquatic ecosystems. They are
predatory carnivores with a wide range of prey species. Approximately 15% of the 3000 species of snakes worldwide
are considered dangerous to humans and responsible for significant morbidity and mortality, (WHO, 2017; Gold, Dart
& Barish, 2002). Snake venom is the most complex of all natural venoms and poisons and it varies in composition
from species to species, (Menez, 2003). There is a considerable variation also in the venom constituents within a single
species within a geographical distribution at different season (Warrel, 1996). The four key toxins of snake venom
includes serine protein (SVSP), three finger toxins (3FTx), phospholipase A 2 (PLA2) and snake venom
metalloproteinease (SVMP), (Casewell & Clevers, 2020).

SBE is a significant public health concern in many parts of the world, particularly in rural areas of tropical and
subtropical regions. It is estimated that 1.7 billion lives globally are threatened annually by neglected tropical diseases,
of which SBE is among them. An estimated 5.4 million cases of SB occur every year, and it is responsible for causing
around 2.7 million cases of envenoming, with about 138,000 cases of deaths and more than 400,000 people suffering
from disabilities around the world, (Warrell, 2010; Williams, et al,. 2019 and WHO, 2019). SBE is most prevalent in
Africa, Asia, and Latin America. According to research, individuals living in rural areas of these regions are at a higher

∗ Corresponding author. Tel.: +23407038262468

E-mail address: solojayuba@[Link] (Joseph S. A.)
[Link] © 2024 Department of Mathematics, Modibbo Adama University.
39 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

risk of SBE due to their occupational exposure to snakes. SBE causes morbidity and mortality among farmers,
pastoralists, hunters, and children, (WHO, 2019; Longbottom, et al., and Gampin, Nassouri, Chippaux & Semde, 2016).
Those who suffered the most from SB are the young agriculturists, the most economically productive group, followed
by children, (Warrell, 2010 & WHO, 2019).

SBE can lead to a wide range of clinical manifestations, including local effects (such as pain, swelling, and tissue
damage), systemic effects (such as bleeding disorders, cardiovascular collapse, and renal failure), and neurotoxic
effects (such as paralysis and respiratory failure). Limited access to healthcare facilities, especially in rural areas, is a
major challenge in the management of SBE. This often results in delayed treatment and increased morbidity and
mortality, (Gutierrez et al., 2017). SBE can be treated when adequately diagnosed. Usually, the clinical outcomes of
the diagnosis of SBE cases for management can be classified as mild, moderate, or severe. The toxicity effects induced
by SBE are currently treated with intravenous administration of anti-snake venom (ASV) in combination with analgesic
fluid therapy, hemodialysis, and/or antibiotics, (Harrison & Williams, 2019; Gutierrez, et al., 2007 and Gutierrez et
al., 2017). The availability and accessibility of effective ASV are crucial for treating SBE. However, some studies have
shown significant challenges related to the production, distribution, and affordability of ASV in many affected regions
(Habib & Brown 2018 and Gutiérrez, 2018).

Despite the dangers posed by snakes, they benefit humans by killing unwanted insects and rodents in food stores
and farms, (Katie & Kennymac, 2013 and Abdullahi, Musa, He & Bello 2021). Snake-skins are used to make shoes,
handbags, and other articles. And snake park serve as tourism attractions. The venom is used for producing life-saving
ASV, biomedical research and other medicinal products, and to control infections such as protozoan infections, cancers,
and so on, (Abubakar et al., 2010; Agom, Moses & Priti,.2022; Alams & Gomes, 2003; Zainal Abidin, Lee, Lekhsan
& Rakesh, 2019 and Hussaini, Okuneye & Gumel, 2017).

It is noteworthy to mention that several mathematical models have been developed to understand the
transmission dynamics and control of some neglected tropical diseases, such as Leishmaniasis, Rabies, Chikungunya,
Dengue, Chagas disease, Onchocerciasis, Lymphatic Filariasis, and Schistosomiasis See, for instance, Adamu &
Hussaini, 2019; Shigui, 2017; Zhang, Jin, Sun, Zhou & Ruan, 2011; Yakob & Clements 2013; Sabiu & Hussaini,
2017; Bhunu & Mushayabasa, 2012; Dumont & Chiroleu 2010; Omondi, Nyabadza, Bonyah, & Badu, 2017; Sabiu,
Hussaini, Zhao & Daihai, 2020 and Woolhouse, 1992) and references therein. Unlike some of those mentioned
neglected tropical diseases, it has been recognized to the best of the authors’ knowledge that not many mathematical
models have been constructed to analyze the dynamics and control of SBE. However, some few researches include
prediction of snake distributions in relation to its transmission using ecological niche modeling, estimation of incidence
of snakebite by application of the law of mass action, Goldstein, Erinjery, Martin, Kasturiratne, Ediriweera & de Silva
(2021) predicting snakebite risk using species distribution model Casewell & Clevers (2020), mathematical model
based on social demography, economic and health infrastructure to examine the impact of behavior on the overall
mortality outcomes of snakebite victims Goldstein et al., (2021) and evaluation of intervention strategies with effective
cost in management of snakebite using mechanistic model (Yañez-Arenas, Peterson, Mokondoko, Rojas-Soto &
Martínez-Meyer, 2014). Some of other few ones can be found in the following studies (Abdullahi, Habib & Hussaini,
2021; Bravo-Vega, Santos-Vega, Juan & Cordovez, 2021; Martin, Erinjery, Ediriweera, de Silva, Lalloo, Iwamura &
Murray, 2022; Kumar, Vedat Suat Erturk, Govindaraj, & Baleanu, 2023 and Abdullahi, Habib. & Hussaini, 2024) and
references therein.

In management of snakebite, effective treatment is associated to proper diagnosis (mild, moderate and severe)
of victim of snakebite envenoming which were not captured in previous SBE models. Therefore, this study aims to
develop a population-based mathematical model to investigate the impacts of the diagnosis of SBE victims and their
treatment with ASV on the dynamics of SBE. A notable feature of the model is the incorporation of the rate of diagnosis
of SBE victims with three different conditions, mild, moderate, and severe. The model also includes the death rate of
snakes due to human activities. The paper is organized as follows: The model is formulated in Section 2, and the
analysis of the mathematical properties of the model and its validation via model fitting are presented in Section 3.
The numerical simulations are explored in Section 4, and the conclusion is presented in Section 5.

© 2024 Department of Mathematics, Modibbo Adama University.

40 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

2. Methods

2.1 Model formulation

In this section, the procedures for formulating the model are presented.

2.1.1 Assumptions of the model

The following are some of the major assumptions made in the formulation of the model:
i. all species of venomous snakes are grouped into a single class;
ii. effective treatment is associated with a proper diagnosis;
iii. recovered individuals are not immune to snakebite thus they transit to the susceptible class after
recovery;
iv. snakebite-induced death is negligible in all treatment classes.

2.1.2 Description of the model

The human population at time t , denoted by ( N ( t ) ) , is divided into nine (9) mutually exclusive
h

compartments, namely: population of susceptible persons ( S ( t ) ) , population of undiagnosed SBE victims ( I ( t ) ) ,

population of SBE victims diagnosed with mild conditions ( I ( t ) ) , population of SBE victims diagnosed with
ml

moderate conditions ( I ( t ) ) , population of SBE victims diagnosed with severe conditions ( I ( t ) ) , population of
mo s

SBE victims receiving treatment with ASV ( T ( t ) ) , population of persons suffering from early adverse reactions

during treatment with ASV ( E ( t ) ) , population of persons who recovered with disability ( R ( t ) ) and population
ar d

of persons who recovered without disability ( R ( t ) ) . Thus, the total human population is given by
w

N h ( t ) = S ( t ) + I ( t ) + I ml ( t ) + I mo ( t ) + I s ( t ) + T ( t ) + Ear ( t ) + Rd (t ) + Rw (t ) . (1)

The total population of venomous snakes is represented by (V ( t ) ) . In this work, susceptible persons refer to
s
those at risk of being bitten by a venomous snake. The population of susceptible persons is generated by a constant
recruitment rate given by  h , recovered persons with disabilities and without disabilities who become susceptible
again at the rates1 and  2 respectively. It diminishes by susceptible persons who become victims of SBE at a rate 
and natural death rate h . Thus, the equation governing the dynamics of the susceptible population is given by
dS
=  h + 1Rd +  2 Rw − (  + h ) S , (2)
dt
where
Vs
= and  is the effective contact rate between venomous snakes and humans.
Nh
The population of undiagnosed SBE victims is generated by susceptible persons who are bitten by snakes and
become victims of SBE as given by  S . The population decreased due to the diagnosis of SBE victims with mild,
moderate and severe conditions at the rates, 1 , 2 and 3 respectively. It also diminishes due to SBE-induced
mortality at a rate  1 and natural mortality at the rate h .. Hence, we have the following equation describing the
dynamics of undiagnosed SBE victims:
dI
=  S − (1 + 2 + 3 + 1 + h ) I . (3)
dt

© 2024 Department of Mathematics, Modibbo Adama University.

41 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

The population of SBE victims diagnosed with mild conditions is generated by those SBE victims diagnosed
with mild conditions at the rate 1 and diminishes by those receiving treatment at the rate  1 and natural mortality at
the rate h .. Thus, we obtain the following equation for SBE victims diagnosed with mild conditions:
dI ml
= 1I − (1 + h ) I ml . (4)
dt
The population of SBE victims diagnosed with moderate conditions is generated by those SBE victims
diagnosed with mild conditions at the rate  2 and the population diminishes by those receiving treatment at a rate  2 ,
and natural mortality at the rate h . Thus, we obtain the following equation:
dI mo
= 2 I − ( 2 + h ) I mo . (5)
dt
The population of SBE victims diagnosed with severe conditions is generated by those SBE victims diagnosed
with severe conditions at the rate 3 , it diminishes by those receiving treatment at a rate  3 , SBE-induced mortality
at the rate  2 , and natural mortality at the rate h . Hence, we have the following equation:
dI s
= 3 I − ( 3 +  2 + h ) I s . (6)
dt
The population of SBE victims receiving treatment with ASV is generated by SBE victims diagnosed with
mild, moderate, and severe conditions at the rates 1 , 2 and  3 , respectively. This population is reduced by
proportions of persons who recovered with or without disability at a rate  , those who suffer from early adverse
reactions at a rate  and natural mortality at a rate h . Hence, we have the following equation:
dT
= 1I mi +  2 I mo +  3 I s − ( +  + h ) T . (7)
dt
The population of persons suffering from early adverse reactions due to ASV treatment is generated by SBE
victims on treatment who suffer from early adverse reactions due to ASV at the rate  . This population diminishes
due to those persons who recovered with or without disability at a rate  and natural mortality at a rate h . Thus, we
have the following equation:
dEar
=  T − ( + h ) Ear . (8)
dt
The population of persons recovered with disability is generated by the proportions of persons who recovered
with disabilities in the treatment class at the rate  and in the early adverse reaction class at the rate  . It
diminishes because of movement to a susceptible class at the rate 1 and natural mortality at the rate h . Hence, we
have the following equation:
dRd
= T + Ear − (1 + h ) Rd . (9)
dt
The population of persons who recovered without disabilities is generated by the proportions of persons in the
treatment class who recovered without disabilities at the rate  (1 −  ) and in the early adverse reaction classat the
rate  (1 −  ) . It diminishes because of transition to susceptible class at the rate  2 and natural mortality at the rate
h . Thus, we have the following equation:

© 2024 Department of Mathematics, Modibbo Adama University.

42 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

dRw
=  (1 −  ) T + (1 −  ) Ear − ( 2 + h ) Rw . (10)
dt
The venomous snake population increases at a rate  s ,which is reduced by natural mortality rate  s and human
activities induce mortality at a rate  . Thus, we obtain the equation governing the dynamics of the venomous snakes
as
dVs
=  s − ( + s )Vs . (11)
dt
2.1.3 The model equation
Therefore, the model describing the dynamics of SBE is given by the following system of first-order nonlinear
ordinary differential equations (the model schematic diagram is depicted in Fig. 1, and the associated state variables
and parameters are presented in Tables 1 and 2, respectively).

dS
=  h + 1 Rd +  2 Rw − (  + P1 ) S ,
dt
dI
=  S − P2 I ,
dt
dI ml
= 1 I − P3 I ml ,
dt
dI mo
=  2 I − P4 I mo ,
dt
dI s
= 3 I − P5 I s ,
dt
dT
=  1 I ml +  2 I mo +  3 I s − P6T ,
dt
dEar
=  T − P7 Ear ,
dt
dRd
= T +  Ear − P8 Rd ,
dt
dRw
= 1T +  2 Ear − P9 Rw , (12)
dt
dVs
=  s − P10Vs ,
dt

with the following initial conditions

S ( 0 )  0, I ( 0 )  0, I ml ( 0 )  0, I mo ( 0 )  0, I s ( 0 )  0, T ( 0 )  0, Ear ( 0 )  0, Rd ( 0 )  0, Rw ( 0 )  0, Vs ( 0 )  0,
where
P1 = h , P2 = (1 + 2 + 3 + 1 + h ) , P3 = (1 + h ) , P4 = ( 2 + h ), P5 = ( 3 +  2 + h ), P6 = ( +  + h ) ,
Vs
P7 = ( + h ) , P8 = (1 + h ) , P9 = ( 2 + h ) , P10 = ( + s ) , 1 = (1 −  ) ,  2 = (1 −  ) ,  = .
Nh
Table 1: Definition of State Variables of the Model

Variable(s) Description

S(t) Population of susceptible persons at time t

© 2024 Department of Mathematics, Modibbo Adama University.

43 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

I(t) Population of undiagnosed SBE victims at time t

Iml(t) Population of SBE victims diagnosed with mild conditions at time t
Imo(t) Population of SBE victims diagnosed with moderate conditions at time t
Is(t) Population of SBE victims diagnosed with severe conditions at time t
T(t) Population of SBE victims receiving treatment with ASV at time t
Ear(t) Population of persons who suffer from the early adverse reaction at time t
Rd(t) Population of persons who recovered with disabilities at time t
Rw(t) Population of persons who recovered without disabilities at time t
Vs(t) Population of venomous snakes at time t

Table 2: Definition of Parameters of the Model

Parameter(s) Description

h ( s ) Recruitment rate into susceptible class(Venomous snakes)

Rate of susceptible persons who become victims of SBE

i ( i = 1, 2,3) Rates of diagnosis of SBE victims that produce mild, moderate and severe
conditions respectively
 i ( i = 1, 2,3) Rate of treatment of SBE victims with mild, moderate and severe conditions
respectively
 ( i = 1, 2 ) SBE-induced death rates in I and Is classes respectively
i
 Proportion of persons who recovered with disabilities in T class
 Proportion of persons who recovered with disabilities in Ear class
 Recovery rate of persons in T class
 Recovery rate of persons in Ear class
 Death rate of venomous snakes due to human activities

h ( s ) Humans (venomous snakes) natural mortality rates

Transition rates to susceptible classof recovered persons with disability (without
1 (  2 ) disability)
Effective contact rate between venomous snakes and humans
 Rate of persons who suffered from early adverse reactions


Figure 1: Schematic diagram of the model

© 2024 Department of Mathematics, Modibbo Adama University.

44 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

3. Results

3.1 Model Analysis

For model analysis, we assume that all the model parameters are positive. Further, the assumption that the
recovered person transit to the susceptible class is relaxed (i.e., 1 = 2 = 0 ) for simplicity and mathematical
convenience.

3.1.1 Basic properties of the model

The basic properties of the model in terms of positivity and boundedness of solution are explored in the following
theorems:
Theorem 1: For any non-negative initial conditions of the model equations given by system (12), the solution
( S , I , I ml , I mo , I s , T , Ear , Rw , Rd ,Vs ) of the model equation will be positive for all t  0 .
Proof. Let
t1 = sup t  0 : S  0, I  0, I ml  0, I mo  0, I s  0, T  0, Ear  0, Rd  0, Rw  0, Vs  0 , therefore, t>
0. Consider the first equation of system (12) as follows:
dS dS
=  h − (  + P1 ) S, since  h  0 we have  − (  + P1 ) S.
dt dt
Separating the variables and integrating both sides from t = 0 to t = t1 we obtain
S ( t1 )  S ( 0 ) exp  −  Pt +   ( y )dy    0, thus, S ( t )  0 , for every t  0 . Using a similar
t1

   
1 1 0

approach it can be shown that the remaining state variables, i.e.,

I ( t )  0, I ml ( t )  0, I mo ( t )  0, I s ( t )  0, T ( t )  0,
Ear ( t )  0, Rd ( t )  0, Rw ( t )  0,Vs ( t )  0 for all t  0 , hence the proof is complete.
Lemma 1. The solution of the model equations given by system (12) is bounded in the closed set
 h 
Ω = Ω h Ω s  9
+  + ,where Ω h = ( S , I , I ml , I mo , I s , T , Ear , Rd , Rw )  9
+ : 0  Nh   and
 h 
 s 
Ω s = Vs  + : 0  Vs  .
 s 
Furthermore, the set Ω is positively invariant with respect to the model equations.
Proof. The total human population is given by
N h = S + I + I ml + I mo + I s + T + Ear + Rd + Rw . (13)
Thus, by differentiating equation (13), we have
dNh dS dI dI ml dI mo dI s dT dEar dRd dRw
= + + + + + + + + , (14)
dt dt dt dt dt dt dt dt dt dt
substituting the right-hand sides of the system (12) into equation (14) and simplifying we obtain
dNh
=  h − Nh h − 1I −  2 I s .
dt
Neglecting SBE-induced mortality in I and Is classes, we have
dNh dNh
  h − Nh h and  − N h h
dt dt

© 2024 Department of Mathematics, Modibbo Adama University.

45 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

Separating the variables and integrating both sides we have

  −  h N h ( 0 )  − h t
h
N h ( t )  N h ( 0 ) e − ht and N h ( t ) 
− h  e . This can be rewritten as
h  h 
   −  h N h ( 0 )  − h t
N h ( 0 ) e − h t  N h ( t )  h −  h  e . It follows that
h  h 

0  lim inf N h ( t )  lim sup N h ( t )  h .Similarly, considering the equation for the venomous snake population
t → t → h
dVs
from system (12), we have =  s − sVs −  Vs . Neglecting the venomous snakes’ death rate due to human
dt
dVs dVs
activities, we have   s − sVs and  − sVs . Separating the variables and integrating both sides, we
dt dt
   s −  sVs ( 0 )  − s t
obtain Vs ( 0 ) e s  Vs ( t )  s − 
− t
 e . It follows that
s  s 
s h s
0  lim inf Vs ( t )  lim sup Vs ( t )  .Now, if Nh ( 0)  and Vs ( 0 )  , then we have
t → t → s h s
h s
Nh (t )  and Vs ( t )  respectively. Thus the region Ω is positively Invariant meaning that the solution
h s
h s
starting in Ω stays in Ω for t>0. In addition, if N h ( 0 )  and Vs ( 0 )  , then either the solution enters the
h s
h s
region Ω in finite time or N h ( t ) approaches asymptotically, and N s ( t ) approaches asymptotically.
h s
10
Thus, the region Ω attracts all solutions in + , so the proof is complete.
Since the region Ω is epidemiologically and mathematically realistic, it suffices to study the dynamics of the
model equations given by system (12) in the region, where the usual existence, uniqueness and continuation of
solution hold.

3.2 Model Fitting and Parameter Estimation

Monthly data on snakebite envenoming for the period 2020–2022 were collected from the snakebite treatment
and research hospital (STRH) in Kaltungo, Gombe State. The data is used to estimate the unknown values of parameters
as well as to validate the model. This is done by using the nonlinear least squares method implemented in Matlab. The
graphs depicted in Figures 2–10 demonstrate how well the model fits the reported data with respect to the cumulative
number of monthly cases of the following: undiagnosed SBE victims; SBE victims diagnosed with mild conditions;
SBE victims diagnosed with moderate conditions; SBE victims diagnosed with severe conditions; SBE victims
receiving treatment with ASV; persons suffering from early adverse reactions; persons recovered with disabilities;
persons recovered without disabilities; and SBE-induced deaths. In addition, the estimated values of the initial
conditions and the parameters in the model are presented in Table 3 and Table 4, respectively. Since the estimated
outcomes of the model are in good agreement with the actual reported data, the model can be used to predict the cases
of SBE and study its dynamics in a given population.

© 2024 Department of Mathematics, Modibbo Adama University.

46 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

Table 3: Reported data and the estimated values of the state variables
Variable(s) Value(s) Reference

S(0) 1.9315×106 Fitted

I(0) 69 Data(STRH)
Iml(0) 9 Data(STRH)
Imo(0) 8 Data(STRH)
Is(0) 29 Data(STRH)
T(0) 46 Data(STRH)
Ear(0) 1 Data(STRH)
Rd(0) 0 Data(STRH)
Rw(0) 46 Data(STRH)
Vs(0) 1.2431×104 Fitted

Table 4: Estimated values of the state model parameters

Parameter(s) Value(s) Unit Reference

Λh 1336 Day-1 (Abdullahi, Habib & Hussaini, 2021)
Λs 0.2206 Day-1 Fitted
𝛳 i(i = 1, 2, 3) 0.2775, 0.1830, 1.0000 Day-1 Fitted
𝜏i(i = 1, 2, 3) 1.0000, 1.0000, 0.4888 Day-1 Fitted
γ 0.0166 Day-1 Fitted
δ1, δ2 0.0014, 1.9066×10−8 Day-1 Fitted
ρ 1.8805×10−4 Nill Fitted
𝜎 0.8751 Day-1 Fitted
ψ 0.1713 Nill Fitted
 0.0537 Day-1 Fitted
5.04×10−6, Day-1 (Abdullahi, Habib & Hussaini, 2021)
µh
2.283×10−4 Day-1 Fitted
µs
0.9973, 0.8290 Day-1 Fitted
1 ,  2 0.9995 Day-1 Fitted
ϕ 0.08367 Nill Fitted
𝛽 4.1040×108 Day-1 Fitted
Ks

© 2024 Department of Mathematics, Modibbo Adama University.

47 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

Figure 2: Graph showing the results of the model fitting with the cumulative number of reported data on
undiagnosed SBE victims.

Figure 3: Graph showing the model results fitting with the cumulative number of reported data on SBE
victims diagnosed with mild conditions.

Figure 4: Graph showing the results of the model fitting with the cumulative number of reported data on SBE
victims diagnosed with moderate conditions.

© 2024 Department of Mathematics, Modibbo Adama University.

48 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

Figure 5: Graph showing the results of the model fitting with the cumulative number of reported data on SBE
victims diagnosed with severe conditions.

Figure 6: Graph showing the results of the model fitting with the cumulative number of reported data on SBE
victims receiving treatment.

Figure 7: Graph showing the results of the model fitting with the cumulative number of reported data on
individuals suffering from early adverse reactions uring treatment.

© 2024 Department of Mathematics, Modibbo Adama University.

49 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

Figure 8: Graph showing the results of the model fitting with the cumulative number of reported data on
persons who recovered with disabilities.

Figure 9: Graph showing the results of the model fitting with the cumulative number of reported data on
persons who recovered without disabilities.

Figure 10: Graph showing the results of the model fitting with the cumulative number of reported data on
SBE-induced deaths.

3.3 Equilibrium Points of the Model

In this section, we establish the conditions for the existence of equilibrium points in the model equations
given by the system (12).
3.3.1 Existence of snakebite envenoming free equilibrium point
The Snakebite-envenoming free equilibrium point describes a circumstance where a community will be free
from SBE. Biologically speaking, this equilibrium point can only exist if there is no contact between venomous snakes

© 2024 Department of Mathematics, Modibbo Adama University.

50 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

and humans, i.e.,  = 0 which also implies that the force of envenoming,  = 0 . Therefore, the existence of snakebite-
*
envenoming free equilibrium point denoted by E is established and is presented as follows:
  
E * = ( S * , I * , I ml
* *
, I mo , I s* , T * , Ear* , Rd* , Rw* ,Vs* ) =  h , 0, 0, 0, 0, 0, 0, 0, 0, s  .
 P1 P10 
3.3.2 Existence of snakebite envenoming endemic equilibrium point
The Snakebite-envenoming endemic equilibrium point demonstrates a situation where SBE will persist in a
given community. Biologically, this can due to interaction between venomous snakes and humans, i.e.,   0 , which
also implies the force of envenoming,   0 . Now, to obtain the conditions for the existence of the snakebite-
envenoming endemic equilibrium point,
Let E
**
= ( S ** , I ** , I ml
** **
, I mo , I s** , T ** , Ear** , Rd** , Rw** ,Vs** ) be an arbitrary endemic equilibrium point of the
model equation given by system (12). In addition, let
V **
 = ** s
, (15)
N h**
be the associated force of envenoming of the model equation at steady state. Now, the equations of the model given by
system (12) are solved in terms of the force of envenoming at a steady state given by equation (15). Solving for the
state variables in the model equation given by system (12) at steady-state gives
h  **  h 1 h  **  2  h  ** 3 h  **
S ** = , I ** = **
, I ml = *
, I mo = , I s** = ,
 + P1 P2 (  + P1 ) P2 P3 (  + P1 ) P2 P4 (  + P1 ) P2 P5 (  ** + P1 )
** ** ** **

  h ( P4 P51 + P3 P52 + P3 P43 )  ** **   h ( P4 P51 + P3 P52 + P3 P43 )  **

T ** = , Ear = ,
P2 P3 P4 P5 P6 (  ** + P1 ) P2 P3 P4 P5 P6 P7 (  ** + P1 )
 h (  +   P7 )( P3 P43 + P3 P52 + P4 P51 )  **
R = **
,
P2 P3 P4 P5 P6 P7 P8 (  ** + P1 )
d

 h (  2 +  1P7 )( P4 P51 + P3 P52 + P3 P43 )  ** **  s

R = **
,Vs =
P2 P3 P4 P5 P6 P7 P9 (  ** + P1 )
w .
P10
Substituting the values of these state variables at steady state into the corresponding force of envenoming given by
equation (15), we obtain the following quadratic equation in terms of  ** , A2 **2 + A1 ** − A0 = 0 ,
(16)
where
A2 =     h P9 P10 ( P4 P51 + P3 P5 2 + P3 P43 ) +     2  h P3 P8 P10 ( P4 P51 + P5 2 + P43 )
+     h P7 P9 P10 ( P4 P51 + P3 P5 2 + P3 P43 ) +   1 h P7 P8 P10 ( P4 P51 + P3 P5 2 + P3 P43 )
+   h P8 P9 P10 ( P4 P51 + P3 P5 2 + P3 P43 ) +   h ( P4 P51 + P3 P5 2 + P3 P43 )
+ P3 P4 P5 P6 P7 P8 P9 (  h P10 +  s P2 ) +  h P6 P7 P8 P9 P10 ( P4 P51 + P3 P52 + P3 P43 ) ,
A1 = P2 P3 P4 P5 P6 P7 P8 P9 (  s P1 +  h P10 −  s  ) ,
A0 =  s  PP
1 2 P3 P4 P5 P6 P7 P8 P9 .

Now, for  0 which also implies 

**
 0 , it is clear that Ai  0, ( i = 0, 2 ) , since all the model parameters are
positive. Also, A1  0 or A1  0 whenever  s P1 +  h P10   s  or  s P1 +  h P10   s  , respectively.

© 2024 Department of Mathematics, Modibbo Adama University.

51 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

Applying the Descartes Rule of Signs in each case showsthat only one sign changes, establishing the existence of one
**
positive real root, i.e., a unique snakebite-envenoming endemic equilibrium point of the model denoted by E is
established. Thus, the following proposition is presented:
Proposition [Link]   0 with either  s P1 +  h P10   s  or  s P1 +  h P10   s  , then the model equation
given by system (12) has a unique snakebite-envenoming endemic equilibrium point denoted by E ( ).
**

3.4 Local Stability of the Equilibrium Point of the Model

The local asymptotic stability of the equilibrium points of the model is established using the method of
linearization, and the results are reported in the following theorems:
Theorem [Link] snakebite-envenoming free equilibrium point, E ( ) *
and the snakebite-envenoming endemic

equilibrium point, E ( ) of the model equation given by system (12) are locally asymptotically stable (LAS) in the
**

region  .
Proof. The Jacobian corresponding to the model equation given by system (12) is given as
  Vs  S
 − N * − P1 0 0 0 0 0 0 0 0 −
N h* 
 h 
  Vs S 
 −P2 0 0 0 0 0 0 0 
 N h* N h* 
 0 1 −P3 0 0 0 0 0 0 0 
,
 
 0 2 0 −P4 0 0 0 0 0 0 
J (X ) = 
 0 3 0 0 −P5 0 0 0 0 0 
 0 0 1 2 3 −P6 0 0 0 0 
  (17)
 0 0 0 0 0  −P7 0 0 0 
 0 0 0 0 0   −P8 0 0 
 
 0 0 0 0 0 1  2 0 −P9 0 
 −P10 
 0 0 0 0 0 0 0 0 0

(
where X = S , I , I ml , I mo , I s , T , Ear , Rd , Rw , Vs ) and Nh = Nh* are taken as a constant for simplicity. Let
 = 0 , so that the snakebite-envenoming free equilibrium, E * exists. Now, the Jacobian of the system (12)
*
evaluated at E is given by
 − P1 0 0 0 0 0 0 0 0 0 
 0 −P 0 0 0 0 0 0 0 0 
 2

 0 1 −P3 0 0 0 0 0 0 0 
 
 0 2 0 −P4 0 0 0 0 0 0 
. (18)
 3 −P5 0 
J ( E* ) = 
0 0 0 0 0 0 0

 0 0 1 2 3 −P6 0 0 0 0 
 0 0 0 0 0  −P7 0 0 0 
 
 0 0 0 0 0   −P8 0 0 
 0 0 0 0 0 1  2 0 −P9 0 
 
 0 0 0 0 0 0 0 0 0 −P10 
The characteristic equation of the Jacobian matrix evaluated at the snakebite-envenoming free equilibrium point given
by equation(18)is presented as

© 2024 Department of Mathematics, Modibbo Adama University.

52 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

det ( J ( E * ) − I10 ) = 0 . (19)

The eigenvalues of equation (19) are presented as follows:
i = − Pi , (i = 1, 2...,10) . (20)
It follows from equation (20) that all the eigen values,  i , for i = 1, 2,...,10 ,have satisfied the negativity
requirement for local asymptotic stability of E * for the reason that Pi  0, for i = 1,2,...,10 , since, all the model
parameters are positive. Thus, the snakebite-envenoming free equilibrium point E * for the system (12) is locally
asymptotically stable. Similarly, suppose that the result in Proposition 1 holds so that a unique snakebite-envenoming
endemic equilibrium point E ** exists. The Jacobian matrix evaluated at the snakebite-envenoming endemic equilibrium
point E ** is presented as
  Vs**  S ** 
 − * − P1 0 0 0 0 0 0 0 0 − 
 Nh N h* 
  V ** S ** 
 s
−P2 0 0 0 0 0 0 0 
 N h* N h*  . (21)
 0 1 −P3 0 0 0 0 0 0 0 

 2 −P4 0 
J ( E ** ) = 
0 0 0 0 0 0 0

 0 3 0 0 −P5 0 0 0 0 0 
 0 0 1 2 3 −P6 0 0 0 0 
 
 0 0 0 0 0  −P7 0 0 0 
 0 0 0 0 0   −P8 0 0 
 
 0 0 0 0 0 1  2 0 −P9 0 
 −P10 
 0 0 0 0 0 0 0 0 0
The characteristic equation of the Jacobian matrix evaluated at the snakebite-envenoming endemic equilibrium point
is given by
det ( J ( E ** ) − I10 ) = 0 . (22)
The eigenvalues of equation (22) are presented as follows:

 j = − Pi , (i = 2,3...,10),( j = 1, 2,3...,9) and 10 = − Vs +*N h P1 .

** *
(23)
Nh
From equation (23), it is obvious that 10  0 . Also, observe that the first nine eigenvalues of equation (22) and
equation (19) are the same. Hence, following similar arguments that established the LAS of the snakebite-envenoming
free equilibrium point, we conclude that the snakebite-envenoming endemic equilibrium point of the model is also
locally asymptotically stable (LAS), and hence the proof is complete.

3.5 Global Stability of Equilibrium Point of the Model

To ensure that the stabilities of snakebite-envenoming free and snakebite-envenoming endemic equilibrium
points are independent of the initial conditions, we shall carry out global stability analysis using Lyapunov function
theory in conjunction with LaSalle’s invariance principle. These will be achieved through the following theorems:
Theorem 3. The snakebite-envenoming free equilibrium point of the model is globally asymptotically stable in the
regionΩ.
Proof. Consider the following candidate Lyapunov function given by:
  V 
( t ) = f1I + f2 I ml + f3 I mo + f4 I s + f5T + f6 Ear + f7 Vs − Vs* − Vs* ln  s*   , (24)
  Vs  
where fi ,(i = 1, 2...,7) are positive constants to be determined. Now, the time derivative of the Lyapunov function
along the solution of system (12) is given by

© 2024 Department of Mathematics, Modibbo Adama University.

53 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

 V* 
(t ) = f1 I + f 2 I ml + f 3 I mo + f 4 I s + f 5T + f 6 Ear + f 7 1 − s  Vs . (25)
 Vs 
Using the right-hand sides of system (12) in equation (24) we have
( t ) = f1 (  S − P2 I ) + f 2 (1I − P3 I ml ) + f3 (2 I − P4 I mo ) + f 4 (3 I − P5 I s )
 V* 
+ f5 ( 1 I ml +  2 I mo +  3 I s − P6T ) + f 6 (  T − P7 Ear ) + f 7 1 − s  (  s − P10Vs ) . (26)
 Vs 
Expanding and collecting like terms in equation (26) and keeping into consideration that S  N h the following
equation is obtained:
( t )  ( f 21 + f32 + f 43 − f1 (1 + 2 + 3 ) ) I − f1 (1 + h ) I + ( f51 − f 21 ) I ml − f 2 h I ml
+ ( f5 2 − f3 2 ) I mo − f3 h I mo + ( f5 3 − f 4 3 ) I s − f 4 ( 2 + h ) I s + ( f 6 − f5 P6 ) T − f 6 P7 Ear
 sVs*
+ ( f1 − f 7 P10 )Vs − f 7 + f 7  s + f 7Vs* P10 . (27)
Vs
From equation (27) the following Lyaponuv coefficients were derived
P6 
f1 = f 2 = f3 = f 4 = f5 = 1, f 6 = and f 7 = .
 P10
(28)
Thus, equation (27) becomes
P6 P7    sVs* 
( t )  − (1 + h ) I − h ( I ml + I mo ) − ( 2 + h ) I s − Ear −  −  s − Vs* P10  (29)
 P10  Vs 
Since  = 0 , it follows that ( t )  0, with ( t ) = 0, if I = I ml = I mo = I s = T = Ear = 0 . Since,  is

positively invariant and attracting, The largest positive invariant set in  is the singleton E .Thus, according to
*

*
LaSalle’s invariance principle [38],the snakebite-envenoming free equilibrium point E is globally attractive.
Therefore, the snakebite-envenoming free equilibrium point is globally asymptotically stable.
Theorem 4. The snakebite-envenoming endemic equilibrium point of the model is globally asymptotically stable in
the region Ω.
Proof. Suppose the result in Proposition 1 holds so that the existence of unique snakebite-envenoming endemic
equilibrium point is assured. Consider the following quadratic Lyapunov function with the following candidate:
G( X ) = 1 ( S − S ** ) + 2 ( I − I ** ) + 3 ( I ml − I ml ) + g24 ( I mo − I mo** ) + g25 ( I s − I s** ) + g26 (T − T ** )
g 2 g 2 g ** 2 2

2 2 2
+ 7 ( Ear − Ear** ) + 8 ( Rd − Rd** ) + 9 ( Rw − Rw** ) + 10 (Vs − Vs** ) ,
g 2 g 2 g 2 g 2
(30)
2 2 2 2
where X = S , I , I ml , I mo , I s , T , Ear , Rd , Rw ,Vs .
Now, differentiating equation (30) along the solution of the model equation given by system (12) and choosing
gi = 1, ( i = 1, 2,3,...10 ) gives
= ( S − S ** ) + ( I − I ** ) + ( I ml − I ml** ) + ( I mo − I mo ) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** )
dG **

dt 
d dV
 S + I + I ml + I mo + I s + T + Ear + Rd + Rw  + Vs − Vs**  s . (31)
dt dt
From system (12), we have
d
 S + I + Iml + Imo + I s + T + Ear + Rd + Rw  = h − h S + I + Iml + Imo + I s + T + Ear + Rd + Rw  − 1I −  2 I s (32)
dt

© 2024 Department of Mathematics, Modibbo Adama University.

54 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

and
dVs
=  s − P10Vs . (33)
dt
Using equations (32) and (33) in equation (31), we obtain
= ( S − S ** ) + ( I − I ** ) + ( I ml − I ml** ) + ( I mo − I mo ) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** )
dG **

dt 
  h − h  S + I + I ml + I mo + I s + T + Ear + Rd + Rw  − 1I −  2 I s  + Vs − Vs**    s − P10Vs  . (34)
At steady state, we have
+ I s** + T ** + Ear** + Rd** + Rw*  + 1I ** +  2 I s** and  s = P10Vs .
**
 h = h  S ** + I ** + I ml
**
+ I mo
**
(35)
Substituting equation (35) in equation (34) gives
= ( S − S ** ) + ( I − I ** ) + ( I ml − I ml** ) + ( I mo − I mo ) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** )
dG **

dt 
 h  S ** + I ** + I ml
 
**
+ I mo
**
+ I s** + T ** + Ear** + Rd** + Rw*  + 1I ** +  2 I s** − h  S + I + I ml + I mo + I s + T + Ear + Rd + Rw 

 −1I −  2 I s  + Vs − Vs**   P10Vs** − P10Vs  . (36)

Collecting like terms gives
= ( S − S ** ) + ( I − I ** ) + ( I ml − I ml** ) + ( I mo − I mo ) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** )
dG **

dt 
−h ( S − S ** ) − h ( I − I ** ) − h ( I ml − I ml** ) − h ( I mo − I mo
**
) − h ( Is − Is** ) − h (T − T ** ) − h ( Ear − Ear** ) − h ( Rd − Rd** ) − h ( Rw − Rw** )

−1 ( I − I ** ) −  2 ( I s − I s** ) − P10 (Vs − Vs** )(Vs − Vs** ) . (37)
 
After expansion and some series of simplifications, the following equation is obtained:
= − ( S − S ** )  h ( S − S ** ) + 1  +  − ( I − I ** )  h ( I − I ** ) +  2  +  
dG
dt    
( ) (
 ) **
( 
 ) ( ** 
 )
− I ml − I ml h  I ml − I ml + 3  +  − I mo − I mo h  I mo − I mo +  4  +  
**  **

− ( I − I )   ( I − I ) +   +  − (T − T )   (T − T ) +   +  
** ** ** **
s  s h  s s  5  h 6

− ( E − E )   ( E − E ) +   +  − ( R − R )   ( R − R ) +   +  
** ** ** **
ar  ar h  ar 
ar 7  d d h d d 8

− ( Rw − Rw** )  h ( Rw − Rw** ) + 9  +  − P10 Vs − V  , ** 2

(38)
  s

where
1 = ( I − I ** ) + ( I ml − I ml** ) + ( I mo − I mo
**
) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** ) ,
 2 = ( S − S ** ) + ( I ml − I ml** ) + ( I mo − I mo
**
) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** ) ,
3 = ( S − S ** ) + ( I − I ** ) + ( I mo − I mo
**
) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** ) ,
 4 = ( S − S ** ) + ( I − I ** ) + ( I ml − I ml
**
) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** ) ,
5 = ( S − S ** ) + ( I − I ** ) + ( I ml − I ml** ) + ( I mo − I mo
**
) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** ) ,
 6 = ( S − S ** ) + ( I − I ** ) + ( I ml − I ml** ) + ( I mo − I mo
**
) + ( I s − I s** ) + ( Ear − Ear** ) + ( Rd − Rd** ) + ( Rw − Rw** ) ,
 7 = ( S − S ** ) + ( I − I ** ) + ( I ml − I ml
**
) + ( I mo − I mo** ) + ( I s − I s** ) + (T − T ** ) + ( Rd − Rd** ) + ( Rw − Rw** ) ,
8 = ( S − S ** ) + ( I − I ** ) + ( I ml − I ml** ) + ( I mo − I mo
**
) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rw − Rw** ) ,
9 = ( S − S ** ) + ( I − I ** ) + ( I ml − I ml** ) + ( I mo − I mo
**
) + ( I s − I s** ) + (T − T ** ) + ( Ear − Ear** ) + ( Rd − Rd** ) ,
 = 1 ( I − I ** ) +  2 ( I s − I s** ) .

© 2024 Department of Mathematics, Modibbo Adama University.

55 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

dG dG
 0 and = 0 if and only if S = S ** , I = I ** , I ml = I ml , I mo = I mo ,
** **
From equation (38) it is clear that
dt dt
I s = I s** , T = T ** , Ear = Ear** , Rd = Rd** , Rw = Rw** , Vs = Vs** . Furthermore, every solution of the model equation
given by system (12) approaches the snakebite-envenoming endemic equilibrium, E ( ) as t →  . Hence, the
**

largest positively invariant set in ( S , I , I , I , I , T , E = 0  is the singleton ( E ) . Hence,

dG  **
ml mo s , R , R ,V )   :
ar d w s
  dt
according to the Lasalle’s invariance principle [38], the snakebite-envenoming endemic equilibrium point is globally
asymptotically stable in the region Ω.

4. Numerical Simulations
The data in Tables 3 and 4 are used as baseline values for the state variable and the parameters of the model.
Numerical simulations to investigate the impacts of the diagnosis of SBE victims and the treatment with ASV on the
the burden of SBE were performed. The diagnosis- and treatment-related parameters were varied as low (25%),
moderate (50%), and high (75%). The results of the simulations of the model showing the impacts of the diagnosis of
SBE victims and the treatment with ASV on the burden of SBE in terms of SBE-induced deaths and disability are
presented in Figures. 11 and 12, respectively. The simulation result in Figure 11 clearly shows that an increase in the
rates of diagnosis of SBE victims combined with ASV treatment resulted in a reduction in the burden of SBE in terms
of SBE-induced deaths presented for one year. Similarly, the outcome of the simulation depicted in Figure 12 reveals
that when rates of diagnosis of SBE victims and treatment with ASV are increased, the number of SBE-induced
disabilities prevented increases. It is worth noting that when the rate of diagnosis of SBE victims and treatment with
ASV is set to zero, the number of SBE-induced deaths and disabilities increases. The results of these numerical
simulations indicate that the diagnosis of SBE victims for proper treatment with ASV will help in the struggle against
SBE in terms of reducing its burden. Thus, the affected community should be encouraged to always visit the hospital
for proper diagnosis of SBE victims and treatment with ASV.

Figure 11: Simulation showing the impact of diagnosis and treatment with ASV on SBE-induced deaths

© 2024 Department of Mathematics, Modibbo Adama University.

56 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

Figure 12: Simulation showing the impact of diagnosis and treatment with ASV on SBE-induced disability

5. Conclusion
In this work, a population-based mathematical model was developed to investigate the impact of the
diagnosis of SBE victims and treatment with ASV. The basic properties of the model in terms of positivity and
boundedness of the model solution were explored. The model is fitted using real data on SBE collected from
snakebite treatment and research hospitals in Kaltungo, Gombe State, in northeast Nigeria. The dynamical
behavior of the model in terms of equilibrium points and their associated stability analyses were carried out.
Furthermore, numerical simulations of the model were performed. Thus, the main qualitative and
epidemiological findings of this study are summarized as follows:

i. the outcome of the model fitting has shown that the model fitted well with the real data, suggesting that
the model can be used to study the dynamics of SBE in a given community;
ii. the qualitative study of the model has established that the model has two equilibrium points, namely,
snakebite-envenoming free and snakebite-envenoming endemic;
iii. the results of the stability analyses revealed that both the two equilibrium points, i.e., snakebite-
envenoming free and snakebite-envenoming endemic, are locally and globally asymptotically stable;
iv. numerical simulations of the model have shown that diagnosis of SBE victims and treatment with ASV
could reduce the burden of SBE in terms of preventing the number of SBE-induced deaths and
disabilities.

Therefore, based on the epidemiological findings in this work, we recommend that people living in a
high-risk community be encouraged to attend health facilities for appropriate diagnosis of snakebite victims
and treatment.

Acknowledgments
The authors appreciate the comments of reviewers and the editor, which were constructive to improve the
clarity of the manuscript. I sincerely express my gratitude to Prof. I. I. Adamu of the Department of Mathematics,
Modibbo Adama University Yola for thorough supervision of this work. A Special thanks to Abdullahi, S. A
(Ph.D) for assisting with the coding and general organization of this paper. We acknowledge Dr. Abubakar S.
Balla of Snakebite Treatment and Research Hospital Kaltungo, Gombe state for updating us (on consulting) with
relevant informations on management of snakebite envenoming patients. Also, we are thankful to the Snakebite
Control Network (S-CONET) for their professional inputs.

© 2024 Department of Mathematics, Modibbo Adama University.

57 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

References
Warrell, D. A. (2010). Snake bite. Lancet, 375(9708), 77–88. [Link] S0140-6736(09)61754-2

Williams D.J., Faiz M.A., Abela-Ridder B., Ainsworth S., Bulfone T.C., Nickerson A.D., Habib A.G., Junghanss T.,
Fan H.W., & Turner M. (2019). Strategy for a Globally Coordinated Response to a Priority Neglected Tropical
Disease: Snakebite Envenoming. PLoS Neglected Tropical Disease, 13:e0007059.

WHO (2019). Snakebite Envenoming: a Strategy for Prevention and Control. World Health Organization.
[Link]

Longbottom J., Shearer F. M., Devine M., Alcoba G., Chappuis F., & Weiss D. J (2018). Vulnerability to Snakebite
Envenoming: a Global Mapping of Hotspots. The Lancet, 392(10148), 673–684. [Link]
6736(18)31224-8 PMID: 30017551

Gampini S., Nassouri S., Chippaux J. P., & Semde R. (2016). Retrospective Study on the Incidence of Envenomation
and Accessibility to Antivenom in Burkina Faso. Journal of Venomous Animals and Toxins Including Tropical
Diseases, 22(1), 1–5. [Link]

Harrison, R. A., & Williams, D. J. (2019). Outlining progress since the first International Snakebite Awareness Day
and some key challenges for next year .Transactions of The Royal Society of Tropical Medicine and Hygiene,
113(10), 577–578. [Link] org/10.1093/trstmh/trz096

Gutierrez, J. M., Bruno, L., Guillermo, L., Alexandra, R., Fernando, C., & Yamileth, A. (2007). Trends in Snakebite
Envenomation Therapy: Scientific, Technological and Public Health Considerations. Current Pharmaceutical
Design, 13 (28), 2935−2950. doi: 10.2174/138161207782023784

Gutierrez, J. M., Calvete, J.J., Habib, A.G., Harrison, R.A., Williams, D.J., & Warrell, D.A. (2017). Snakebite
Envenoming. Nature Review Disease Primers, 3. e17063.

Habib, A. G & Brown, N. I. (2018). The Snakebite Problem and Antivenom Crisis from a Health-Economic
Perspective. Toxicon, 150, 115–123. [Link]/10.1016/[Link].2018.05.009

Gutiérrez, J. M. (2018). Global Availability of Antivenoms: The Relevance of Public Manufacturing

Laboratories. Toxins (Basel), 11(1) . December 24; pii: E5 10.3390/toxins11010005

WHO (2017). Prevalence of Snakebite Envenoming. World Health Organization.

[Link]

Gold, B. S., Dart, R. C., Barish, & R. A. (2002). Bites of Venomous Snakes. New England Journal of Medicine, 347,
347–356. DOI: 10.1056/NEJMra013477

Beri, D. Bhaumik, S. - Snakes, the Ecosystem, and Us: It’s Time we Change. The George Institute of Global Health.
July 2023. Available online at [Link]. The article is licensed under CC BY-NC 2.0.

Katie, V. and Kennymac (2013) Wildlife Epidemiology Laboratory. Why You Should Care about Snakes.
[Link] .Accessed September, 2023.

Abdullahi, Z. U., Musa, S. S., He, D. & Bello, U. M. (2021). Antiprotozoal Effect of Snake Venoms and their
Fractions. Pathogens, 10, 1632. [Link]/10.3390/pathogens10121632.

© 2024 Department of Mathematics, Modibbo Adama University.

58 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

Abubakar, S. B., Abubakar, I. S., Habib, A. G., Nasidi, A., Durfa, N., Yusuf, P. O., Larnyang, S., Garnvwa, J.,
Sokomba, E., Salako, L., Theakston, R. D., Juszczak, E., Alder, N., & Warrell, D. A. (2010). Randomized
Double Blind Comparative Trial of Two Antivenoms for the Treatment of Patients Envenomed by the Saw-Scaled
or Carpet Viper (Echisocellatus) in Northern Nigeria. PLoS Neglected Tropical Disease, 4, e767.
[Link]

Agom, D. I., Moses, B. A., & Priti, M. (2022, December, 3). Tackling Snakebite at the Community Level to
Mitigate Disability Worldwide. International Day of People with disability presented at the annual conference of
the Int’l Society for Neglected Tropical Disease. International Society for Neglected Tropical Diseases.
[Link]/isntd-connect. Availble on [Link]

Alams, M. I., & Gomes, A. (2003) Snake venom Neutralization by Indian Medicinal Plant (Vitexnegundo and
Emblicaofficinalis) Root Extracts. Journal of Ethnopharmcology, 86, 75-80.
[Link]

Zainal Abidin, S. A., Lee, Y. Q., Lekhsan, Y. O. & Rakesh, N. (2019). Malaysian Cobra Venom: a Potential
Source of Anti-Cancer. Toxins, 11.

Hussaini, N., Okuneye, K. & Gumel, A. B. (2017). Mathematical Analysis of a Model for Zoonotic Visceral
Leishmaniasis. Infectious Disease Modelling, 2, 455-474

Adamu, L. & Hussaini, N. (2019). An Epidemic Model of Zoonotic Visceral Leishmaniasis with Time Delay.
Journal for Nigeria Society of Physical Science, 1, 20−29.

Shigui, R. (2017). Modeling the Transmission Dynamics and Control of Rabies in China. Mathematical Biosciences,
286, 65 - 93.

Zhang, J., Jin, Z., Sun, G-Q., Zhou, T., Ruan, S. (2011). Analysis of Rabies in China: Transmission Dynamics and
Control. PLoS ONE, 6(7): e20891.

Yakob, L., Clements, A. C. A. (2013). A Mathematical Model of Chikungunya Dynamics and Control: The Major
Epidemic on Reunion Island. PLoS ONE, 8(3):e57448.

Sabiu, S.M. & Hussaini, N. (2017). Modeling the Dynamics of Dengue Virus in Human and Mosquito Populations.
Journal of Nigerian Association of Mathematical Physics, 43.

Bhunu, C. P & Mushayabasa, S. (2012). Transmission Dynamics of Lymphatic Filariasis:a Mathematical Approach.
ISRN Biomathematics, [Link]

Dumont, Y. & Chiroleu, F. (2010). Vector Control for the Chikungunya disease. Mathematical Biosciences and
Engineering, 7.

Omondi, E.O., Nyabadza, F., Bonyah, E. & Badu, K. (2017). Modeling the Infection Dynamics of Onchocerciasis
and its Treatment. Journal of Biological Systems, 25(2), 1−31. doi:10.1142/S0218339017500139.

Sabiu, S.M., Hussaini, N., Zhao, S. & Daihai H. (2020). Dynamical Analysis of Chikungunya and Dengue
Coinfection Model. Discrete and Continuous Dynamical System Series B, 25(5), 1907-1933.

Woolhouse, M.E.J. (1992). On the Application of Mathematical Models of Schistosome Transmission Dynamics. II
Control. Acta Trop. 50, 189-204.

© 2024 Department of Mathematics, Modibbo Adama University.

59 Joseph et al. (2024) International Journal of Development Mathematics Vol 1 Issue 3| 038 - 059

Bravo-Vega, C.A., Cordovez, J.M., Renjifo-Ibanez, C., Santos-Vega, M., & Sasa, M. (2019). Estimating Snakebite
Incidence from Mathematical Models: A Test in Costa Rica. PLoS Neglected Tropical Disease, 13(12):
e0007914.

Kim S. (2020). Introduction of a Mathematical Model to Characterize Relative Risk of Snakebite Envenoming, a
Neglected Tropical Disease. [Thesis]: the degree of Honors Baccalaureate of Science in Kinesiology: Oregon
State University, University Honors College.

Goldstein, E., Erinjery, J.J., Martin, G., Kasturiratne, A., Ediriweera, D.S. & de Silva, H.J. (2021).Integrating Human
Behavior and Snake Ecology with Agent-Based Models to Predict Snakebite in High Risk Landscapes. PLoS
Neglected Tropical Disease, 15(1): [Link]

Abdullahi, S.A., Habib, A.G. and Hussaini, N.: (2021). Control of Snakebite Envenoming: A Mathematical
Modeling Study. PLOS Neglected Tropical Diseases 15(8): e0009711.
[Link] (2021).

Bravo-Vega, C., Santos-Vega, M & Juan, Cordovez, J. M. (2021). Entangling Snakebite Dynamics: the
Spatiotemporal Role of Rainfall on Snake Envenoming in Colombia, 2021,Cold Spring Harbor Laboratory,
doi: 10.1101/2021.07.13.452215

Martin, G., Erinjery, J.J., Ediriweera, D., de Silva, H.J., Lalloo, D.G., Iwamura, T. & Murray, K.A. (2022). A
Mechanistic Model of Snakebite as a Zoonosis Envenoming Incidence is Driven by Snake Ecology,
Socioeconomics and its Impacts on Snakes. PLoS Neglected Tropical Disease 2022 May 12 (5) 16:e0009867.
doi: 10.1371/[Link].0009867.

Kumar,P. Vedat Suat Erturk, V. Govindaraj, & Baleanu, D. (2023). A Study on the Nonlinear Caputo-Type
Snakebite Envenoming Model with Memory. Computer Modeling in Engineering & Sciences, 136(3), 2487–2506.
[Link]

Abdullahi, S.A., Habib, A.G. & Hussaini, N. (2024). Mathematical Analysis for the Dynamics of Snakebite
Envenoming. Afrika Matematika, 35, 16 (2024). [Link]

Yañez-Arenas, C., Peterson, A. T., Mokondoko, P., Rojas-Soto, O. & Martínez-Meyer, E. (2014). The Use of
Ecological Niche Modeling to Infer Potential Risk Areas of Snakebite in the Mexican State of Veracruz. PLoS
ONE, 9(6). [Link]

Masoud, Y., Anooshe, K., and Ali, K., & Nima. N. (2020). Applying Species Distribution Models in Public Health
Research by Predicting Snakebite Risk Using Venomous Snakes’ Habitat Suitability as an Indicating Factor.
Nature Search, 10(1). [Link]

Casewell, N. R., Clevers, H. (2020). Snake Venom Gland Organoids. Cell, 80(2), 233-247.
doi: 10.1016/[Link].2019.11.038.

Menez, A (2003) The Subtle Beast. Snake, from Myth to Medicine, 1st ed, CRC Press.
[Link]

Warrel, D. A. (1996). Animal Toxins. In Gordon Cook. Manson's Tropical Diseases, 20th Ed. London: WB Saunders.

© 2024 Department of Mathematics, Modibbo Adama University.