Investigating the Role of Antimalarial Treatment and Mosquito Nets in Malaria Transmission and Control through Mathematical Modeling

1  Introduction

Malaria is a contagious disease caused by Plasmodium parasites, mainly transmitted to humans through the bite of a female mosquito [1]. Although there has been significant progress in decreasing the number of malaria cases worldwide, the disease continues to pose an important public health challenge, especially in tropical and subtropical regions. The groups at risk include children under the age of five, pregnant women, and travelers who lack immunity [1,2]. Travelers from countries without local malaria transmission are at significant risk of contracting malaria and its consequences when traveling abroad, as they lack immunity. Additionally, individuals who leave malaria-endemic countries and immigrate to malaria-free nations are also at risk when they return to their home countries to visit friends and family, as their immunity wanes or becomes nonexistent. However, Sub-Saharan Africa remains the most vulnerable region, representing a significantly high share of deaths attributed to malaria [3].

Female mosquitoes that transmit disease have salivary glands that contain sporozoites. In the human body, these sporozoites develop into two forms: erythrocytic and hepatocyte schizogony. During the hepatocyte phase, the sporozoite transforms into a hypnozoite, a dormant form, which later develops into a schizont [4]. The schizont undergoes rapid cell division, producing merozoites that infect red blood cells. Meanwhile, hypnozoites can remain dormant for an extended period. The merozoites produced in the liver invade red blood cells, thereby sustaining the cycle of infection [5]. Trophozoites form and associate with merozoites, developing into blood schizonts that produce merozoites to infect new red blood cells. In some plasmodium parasites, hypnozoites develop into liver schizonts that burst, releasing a large number of merozoites after 8 to 10 months. The number of sporozoites injected by mosquitoes into humans is a crucial factor in the recurrence of malaria [6].

To mitigate the impact of malaria on the global population, various scientific initiatives have been undertaken, including the development of mathematical models [7,8]. The SEIR model, introduced by Kermack and McKendrick, has been widely adopted for modeling malaria transmission [9–11]. Since then, numerous mathematical models have been developed to combat malaria and reduce its global mortality rate [12–14]. However, predicting the severity of malaria remains a challenging task despite these efforts. As climate change awareness grows, it’s essential to recognize the significant impact of climatic and environmental conditions on the spread and transmission of vector-borne diseases [15]. Environmental factors like temperature, rainfall, humidity, wind speed, and daylight duration play a crucial role in shaping the seasonal and daily rhythms of mosquito populations, influencing their ecological and behavioral aspects [16].

Mathematical modeling has become a powerful tool in understanding the transmission dynamics of malaria [17–21]. In [22], the primary focus was on the development of a mathematical model and the evaluation of two disease control strategies: human precautions and mosquito sprays. The study determined that the system achieves a disease-free state when at least 40% of susceptible humans undertake the necessary precautions. In another study [23], the authors proposed an effective disease control approach known as awareness-based intervention, which utilizes media for malaria management. The model incorporates time-dependent control measures for treatment, insecticides, and social media initiatives to minimize the costs associated with malaria control. A new mathematical model using partial differential equations was developed for malaria in [24]. The authors studied the effectiveness of spraying in malaria control and found that uniform spraying produced comparable results to a non-spatial model.

Most mathematical models overlook the mosquito life cycle, excluding eggs, larvae, and pupae from the transmission cycle. While this simplification is helpful, it leads to a significant limitation: these models often struggle to accurately predict the magnitude of malaria in regions where it is most prevalent. To comprehensively understand the complex processes driving malaria transmission and develop more accurate predictive models, it is essential to consider the mosquito life cycle and the impacts of seasonality [25]. Moulay et al. [26] recently developed a mathematical model that simulates the dynamics of mosquito populations, incorporating the self-regulatory mechanisms of the egg and larval stages. They identified a threshold that plays a crucial role in controlling mosquito population growth. However, these models often overlook important factors like the life cycle of mosquitoes and variations in the environment. This can lead to inaccuracies in predicting the severity of the disease, particularly in areas where malaria is widespread.

In comparison to earlier models [22–26] that either assume a constant mosquito population or neglect the impact of treatments, this study addresses these shortcomings by integrating the dynamics of mosquito populations and seasonal variations into a non-linear mathematical model. It also presents a more comprehensive model that incorporates both symptomatic and asymptomatic human infections, considers the different stages of the mosquito life cycle, and integrates treatment effects with vector control. This approach improves realism and facilitates the examination of integrated control strategies. Overall, the combination of stability analysis, detailed compartmental dynamics, and optimal control within realistic constraints represents a significant advancement that differentiates it from existing approaches in both scope and suitability. We will explore the most efficient approaches to managing the disease, emphasizing the utilization of anti-malaria medications, vector control methods to reduce mosquito populations, and the use of mosquito nets.

The sections of the article are arranged as follows. Section 2 outlines the formulation of the model, which divides human and mosquito populations into different compartments and uses a set of nonlinear ordinary differential equations (ODEs) to represent the spread of the disease. In Section 3, we present theoretical findings, including the existence, uniqueness, positivity, and boundedness of the model’s solutions. The stability analysis of disease-free equilibrium (DFE), and endemic equilibrium (EE) with respect to reproduction number R0 is also presented in this section. Section 4 focuses on sensitivity analysis, identifying the most crucial parameters influencing R0. This section also uses numerical simulations to demonstrate the system’s dynamics under various R0 scenarios, illustrating disease-free and endemic conditions. Section 5 expands the model to incorporate optimal control strategies, including the use of anti-malaria drugs, vector population control strategies, and mosquito nets. Section 6 formulates and solves an optimal control problem using Pontryagin’s Maximum Principle to minimize both malaria transmission and intervention costs. The paper concludes by identifying the most effective strategies for malaria control, supported by numerical results.

2  Mathematical Modeling of Malaria

We develop an extended SEIR-based mathematical model that incorporates asymptomatic human cases and vector populations, capturing the complex dynamics of malaria transmission more effectively than classical models. The SEIR type models are particularly effective for modeling malaria because they include a latency period in both human and mosquito populations. The model will analyze how human and mosquito populations change over time, with humans represented as Nh(t) and mosquitoes as Nv(t). We categorize mosquitoes into two populations: immature and mature. The mature mosquito population is divided into three compartments: susceptible (Sv), exposed (Ev), and infected (Iv). Similarly, the human population is divided into five compartments: susceptible (Sh), exposed (Eh), infected (Ih), asymptomatic (Ah), and recovered (Rh).

The infection period of mosquitoes is directly correlated with their relatively short life. Therefore, our model focuses exclusively on mature mosquitoes, as their maturity aligns with their infection period. At any given time t, the combined population consists of mature mosquitoes and humans, represented by:

Nv(t)=Sv(t)+Ev(t)+Iv(t),(1)

and

Nh(t)=Sh(t)+Eh(t)+Ih(t)+Ah(t)+Rh(t).(2)

The growth rate of mosquitoes is represented by Λv. When an uninfected female mosquito bites a human infected with a virus or parasite, the mosquito becomes a carrier of the infectious agent. The force of infection in this case is represented by:

σcSvIhNh,

and hence the susceptible mosquitoes move to the exposed class at this rate. The exposed mosquitoes become infectious at the rate νv. The mosquito population dies naturally at the rate μv.

The susceptible humans are recruited at a constant birth rate Λh. Once the infected mosquito bites a susceptible human, it can transmit the parasite to that human, potentially resulting in an infection. Thus, susceptible humans are exposed to the virus at the rate given as:

σbShIvNh.

The humans exposed to the virus become infectious, which can further be divided into two subgroups: symptomatic and asymptomatic. The rate at which symptomatic individuals develop symptoms is denoted by β1, and the rate at which asymptomatic individuals become infectious is denoted by β2. After individuals recover from infection, they gain immunity from the disease. Both symptomatic and asymptomatic individuals move to the recovered class, represented as Rh, after acquiring immunity at rates θ and γ, respectively. The rate of natural death is denoted by μh, while the rate of death due to the disease is denoted by δh. The above considerations are depicted by the flow diagram shown in Fig. 1.

Figure 1: Diagram showing the transmission of disease through the compartments

The malaria parasite transmission flow diagram is transformed into the following system of nonlinear ODEs:

dShdt=Λh−σbShIvNh−μhSh,(3a)

dEhdt=σbShIvNh−(β1+β2+μh)Eh,(3b)

dIhdt=β1Eh−(θ+μh+δh)Ih,(3c)

dAhdt=β2Eh−(μh+γ)Ah,(3d)

dRhdt=θIh+γAh−μhRh,(3e)

dSvdt=Λv−σcSvIhNh−μvSv,(3f)

dEvdt=σcSvIhNh−(νv+μv)Ev,(3g)

dIvdt=νvEv−μvIv,(3h)

with initial conditions that are non-negative, i.e.,

Sh(0)=Sh0>0,Eh(0)=Eh0≥0,Ih(0)=Ih0≥0,Ah(0)=Ah0≥0,Rh(0)=Rh0≥0,Sv(0)=Sv0>0,Ev(0)=Ev0≥0,Iv(0)=Iv0≥0.(3i)

Table 1 presents the model parameters’ descriptions and values. In compact form, we put the model (3) as follows:

dΨdt=H(Ψ(t)),Ψ(0)=Ψ0, 0<t<tf<∞,(4)

where Ψ(t):[0,tf]→R+8 and H:R+8→R+8 are function defined by:

Ψ(t)=(Sh(t),Eh(t),Ih(t),Ah(t),Rh(t),Sv(t),Ev(t),Iv(t))T,

with

Ψ0=(Sh(0),Eh(0),Ih(0),Ah(0),Rh(0),Sv(0),Ev(0),Iv(0))T,

and

H(Ψ(t))=(H1H2H3H4H5H6H7H8)=(Λh−σbShIvNh−μhShσbShIvNh−(β1+β2+μh)Ehβ1Eh−(θ+μh+δh)Ihβ2Eh−(μh+γ)AhθIh+γAh−μhRhΛv−σcSvIhNh−μvSvσcSvIhNh−(νv+μv)EvνvEv−μvIv).

3  Theoretical Results

In this section, we present the main theoretical results, demonstrating the physical viability of the proposed mathematical model and its potential for further exploration in numerical simulations. We establish the uniqueness of the solution of the model (3) and verify that the solution of each state equation remains positive and bounded for all t≥0. Furthermore, we determine the equilibrium points and the basic reproduction number (ℛ0) to assess the local and global stability of the proposed model at these points.

3.1 Positivity and Boundedness of the Solutions

We aim to establish that the variable Ψ=(Sh,Eh,Ih,Ah,Rh,Sv,Ev,Iv) remains positive and bounded for t≥0 in the analysis of the fundamental characteristics of the malaria transmission model.

Theorem 1: The solution Ψ=(Sv,Ev,Iv,Sh,Eh,Ih,Ah,Rh) to the model has been shown to maintain positivity under consideration of the initial condition provided in Eq. (3i).

Proof: The Eq. (3f) can be rewritten as:

dSvdt+(σcIhNh+μv)Sv=Λv,(5)

or

dSvdt+(Φ(t)+μv)Sv=Λv,(6)

where

Φ(t)=σcIhNh.

The differential Eq. (6) yields us an integrated factor exp(μvt+∫0tΦ(x)dx). As a result, the Eq. (6) can be written as:

ddt[exp(μvt+∫0tΦ(x)dx)Sv]=Λv exp(μvt+∫0tΦ(x)dx).

Integrating the above equation over the interval [0,τ], where τ=max{t>0, Ψ(t)>0}, we determine:

Sv(τ) exp(μvτ+∫0τΦ(x)dx)−Sv(0)=Λv∫0τexp(μvt+∫0tΦ(x)dx)dt.

After simplification, we obtain the following expression for Sv.

Sv(τ)=Sv(0) exp(−(μhτ+∫0τΦ(x)dx))+exp(−(μvτ+∫0τΦ(x)dx))(Λv∫0τexp(μvt+∫0tΦ(x)dx)dt).(7)

Given Sv(0)≥0, Eq. (7) implies that Sv(t)>0 for all t∈[0,τ]. Following a similar procedure, the positivity for all other variables can be demonstrated. □

Theorem 2: The solution Ψ(t) for the malaria disease model (3) is bounded for t≥0.

Proof: Differentiating Eqs. (1) and (2) with respect to time t and then incorporating the equations of disease model (3), we obtain the following ODEs for total vector and human populations.

dNvdt=Λv−μvNv,(8)

and

dNhdt=Λh−δhIh−μhNh≤Λh−μhNh,(9)

with

Nv(0)=Sv(0)+Ev(0)+Iv(0),

and

Nh(0)=Sh(0)+Eh(0)+Ih(0)+Ah(0)+Rh(0).

Suppose for any given initial condition, if

Nv(0)≤Λvμv,(10)

then

dNvdt≤Λv−μvNv.

Similarly if

Nh(0)≤Λhμh,(11)

then

dNhdt≤Λh−μhNh.

We employ the Grönwall’s inequality to determine the following solutions.

Nv(t)≤Λvμv+(Nv(0)−Λvμv)exp((−μv)t),

and

Nh(t)≤Λhμh+(Nh(0)−Λhμh)exp((−μh)t).

This implies that

Nv(t)≤Λvμv, for all t≥0,

and

Nh(t)≤Λhμh, for all t≥0.

Thus, we can write the following inequalities.

limt→∞Nv(t)≤Λvμv,limt→∞Nh(t)≤Λhμh,

which proves the boundedness of Ψ(t) for every t≥0. □

As a result of the above two theorems, the feasible region for the proposed model (3) is defined as follows.

Δ={Ψ∈R+80<Nv≤Λvμv, 0<Nh≤Λhμh, 0≤t≤τ<∞}.

3.2 Existence and Uniqueness

To demonstrate the existence and uniqueness of the solutions of the disease model (3), we present the following foundational results of [27,28].

Theorem 3: Consider the domain D={(t,ν1,…,νm)|a≤t≤b,−∞<νk<∞, for each k=1,2,…,m}. Suppose that each function 𝒢k(1,ν1,…,νm), for k=1,…,m, is continuous and satisfies the Lipschitz condition on D. Then, the system:

dνdt=𝒢(t,ν), ν(a)=ν0,

where ν=(ν1,…,νm)T, ν0=((ν0)1,…,(ν0)m)T, 𝒢=(𝒢1,…𝒢m)T has a unique solution ν for t∈[a,b].

Theorem 4: If 𝒢(t,ν) and ∂𝒢∂νk are in C(D) and if

|∂𝒢∂νk|<ℳ,

then for ℳ>0, 𝒢 satisfies the Lipschitz condition on D for each k = 1, …, m in D.

We now introduce the theorem that establishes the existence of a unique solution for IVP (4).

Theorem 5: Assume that the domain for the problem (4) is Δ⊂D. The system (4) has a unique solution bounded in Δ if H satisfies the Lipschitz condition on Δ.

Proof: Here, the function H(Ψ(t)) and its partial derivative are continuous since the state variable Ψ(t) is continuously differentiable for 0≤t≤τ<∞. Additionally, Theorems 1 and 2 state that H(Ψ(t)) is bounded. To demonstrate Lipschitz’s condition on Δ, all that is required is proving the boundedness of the derivatives ∂Hi∂Ψk, where i,k=1,…,8.

|∂H1∂Sh|=|−(σcIvNh+μh)|≤ℳ<∞,|∂H1∂Iv|=|−(σcShNh)|≤ℳ<∞,|∂H2∂Sh|=|(σcIvNh)|≤ℳ<∞,|∂H2∂Iv|=|(σcShNh)|≤ℳ<∞,|∂H2∂Eh|=|−(β1+β3+μh)|≤ℳ<∞,|∂H3∂Eh|=|(β1)|≤ℳ<∞,|∂H3∂Ih|=|−(θ+δh+μh)|≤ℳ<∞,|∂H4∂Eh|=|(β2)|≤ℳ<∞,|∂H4∂Ah|=|−(γ+μh)|≤ℳ<∞,|∂H5∂Ih|=|(θ)|≤ℳ<∞,|∂H5∂Ah|=|(γ)|≤ℳ<∞,|∂H5∂Rh|=|−(μh)|≤ℳ<∞.|∂H6∂Sv|=|−(σcIhNh+μv)|≤ℳ<∞,|∂H6∂Ih|=|(σcSvNh)|≤ℳ<∞,|∂H7∂Ih|=|(σcSvNh)|≤ℳ<∞,|∂H7∂Sv|=|(σcIhNh)|≤ℳ<∞,

|∂H7∂Ev|=|−(νv+μv)|≤ℳ<∞,|∂H8∂Iv|=|−(μv)|≤ℳ<∞,|∂H8∂Ev|=|(νv)|≤ℳ<∞,

All other first-order derivatives are equal to zero, i.e.,

|∂Hi∂Ψk|=0<ℳ<∞.

Here ℳ=k¯ℳ¯ where k¯>0 is the largest value of model parameters and ℳ¯=max{ℳ1,…,ℳ8}. Here, ℳ1,ℳ2,…,ℳ8 are respectively the upper bounds for each of the state variables Ψk, k=1,2,…,8.

Thus,

|∂Hi∂Ψk|≤ℳ, ∀i,k=1,…8.

This proves the existence of a unique solution of the disease model (4). □

3.3 Equilibrium Points

We determine the Disease-Free Equilibrium (DFE) point by setting Ih=0 and Iv=0 in the steady-state equation of the model (3). We solve the steady-state equations to find the following DFE point.

E0=(Sh0,Eh0,Ih0,Ah0,Rh0,Sv0,Ev0,Iv0)=(Λhμh,0,0,0,0,Λvμv,0,0).(12)

The endemic equilibrium point represents a stable state in which the disease continues to persist and circulate within the population with a constant number of infected individuals who coexist within the population. In this scenario, the disease persists in both the mosquito and human populations, which means Iv≠0 and Ih≠0. To determine the endemic equilibrium (EE) point, we solve the steady-state equations of the model (3) and determine the following EE point.

E1=(Sh1,Eh1,Ih1,Ah1,Rh1,Sv1,Ev1,Iv1),(13)

where

Sh1=[Λh(λh1+μh)],Eh1=[λh1Λhk1(λh1+μh)],Ih1=[β1λh1Λhk1k2(λh1+μh)],Ah1=[β1β2λh1Λhk1k2k3(λh1+μh)],Rh1=[β1λh1Λh(θk3+γβ2)k1k2k3μh(λh1+μh)],Sv1=[Λv(λv1+μv)],Ev1=[λv1Λvk4(λv1+μv)],Iv1=[νvλv1Λvk4μv(λv1+μv)],

and

k1=(β1+β2+μh), k2=(μh+δh+θ), k3=(μh+γ), k4=(νv+μv),λv1=σcIh1Nh, λh1=σbIv1Nh.

3.4 Reproduction Number

The reproduction number ℛ0 is a crucial mathematical statistic that is used to assess the spread of a disease. The disease will spread when ℛ0>1, while ℛ0<1 indicates that the disease will eventually die out. The next-generation matrix approach, introduced by Diekmann and Heesterbeek in 1970 [29,30], is widely used to calculate ℛ0. This approach defines ℛ0 as the spectral radius of the product of two matrices: F, the Jacobian of the new infection rate, and V−1, the inverse of the Jacobian of the remaining transition terms in the disease spread equation. In simpler terms, ℛ0 is the largest absolute eigenvalue of the next-generation matrix FV−1, which represents the potential for the spread of the disease.

To determine ℛ0, we analyze the equations for infectious states in both vector and host populations (Eh,Ih,Ah,Ev,Iv) and extract the matrices ℱ and 𝒱 to compute the next-generation matrix FV−1. The matrices ℱ and 𝒱 are given as follows.

ℱ=(σbShIvNh00σcSvIhNh0),𝒱=((β1+β2+μh)Eh−β1Eh+(θ+μh+δh)Ih−β2Eh+(γ+μh)Ah(νv+μv)Ev−νvEv+μvIv).

The jacobian matrices F and V are then determined as follows:

F=(∂ℱi∂yi)E0,i=1,2,3,4,5,V=(∂𝒱i∂yi)E0,i=1,2,3,4,5,

where (y1,y2,y3,y4,y5)=(Eh,Ih,Ah,Ev,Iv).

Thus, ℛ0, representing the spectral radius of FV−1, is given as:

ℛ0=σcΛvμhβ1Λhμv(θ+μh+δh)(β1+β2+μh).(14)

The reproduction number, ℛ0, determines the transmissibility of a disease. If ℛ0 is greater than one, an outbreak will occur, and the disease will spread. On the other hand, if ℛ0 is less than one, the disease will decline and eventually die out. The threshold value of ℛ0=1 marks the tipping point, where the disease is neither spreading nor declining.

3.5 Stability Characterization of Model

This section investigates the stability of the malaria mathematical model, analyzing both local and global stability at the equilibrium points of DFE and EE.

3.5.1 Evaluation of Local Stability

To evaluate the local stability of the model at the DFE point, the Jacobian matrix is calculated and evaluated at E0, i.e.,

JE0=(−μh000000−σb0−k100000σb0β1−k2000000β20−k3000000θγ−μh00000−σcΛvμhμvΛh00−μv0000σcΛvμhμvΛh000−k40000000νv−μv).(15)

We introduce the following theorem to analyze the local stability of model (3) at the DFE point E0.

Theorem 6: State model (3) is locally asymptotically stable at E0 when ℛ¯0<1, and is not stable for ℛ¯0>1.

Proof: The Jacobian matrix (15) has the following eigenvalues.

λ1=−μhλ2=−μhλ3=−k1λ4=−k2λ5=−k3λ6=−μvλ7=−k4λ8=−μv(1−ℛ¯0).

The eigenvalues (λ1 to λ8) are all negative when ℛ¯0<1, indicating local asymptotic stability at E0, where ℛ¯0=σνvbμv(μv+νv)ℛ0. □

3.5.2 Evaluation of Global Stability

This section presents a global stability analysis of the model, focusing on the disease-free equilibrium and the epidemic equilibrium points. We employed the Castillo-Chavez approach to investigate global stability at the DFE point (E0), utilizing Lyapunov theory and LaSalle’s invariance principle to examine global stability at the EE point (E1).

To establish global stability at the DFE point (E0), we employ the Castillo-Chavez approach, as described in [31]. This involves transforming the model into a suitable form, allowing us to apply this methodology and prove global stability. The model is reformulated in terms of two variables: 𝒳=(Sh,Sv), representing the susceptible individuals (mosquitoes and humans) who are not infected, 𝒴=(Eh,Ih,Ah,Ev,Iv), representing the individuals (mosquitoes and humans) in various infection states. Thus, we have the following reformulation of the model (3).

d𝒳dt=𝒢(𝒳,𝒴),d𝒴dt=ℒ(𝒳,𝒴),𝒴(𝒢,0)=0.(16)

Here, 𝒳 represents the susceptible individuals as a scalar, while 𝒴 is a 5-dimensional vector representing the infected individuals in various states. Note that the recovered class is excluded from the analysis, as it does not interact with the other compartments and therefore does not influence the disease’s dynamics. The DFE point is denoted as E0=(𝒳0,0)=(Λhμh,0,0,Λvμv,0,0,0,0). The following requirements are fulfilled to prove Global Asymptotic Stability (GAS) at the DFE point.

(G1)For d𝒳dt=𝒢(𝒳,0)=0, 𝒳0 is GAS.(17)

(G2)ℒ(𝒳,𝒴)=ℬ𝒴−ℒ¯(𝒳,𝒴) where ℒ¯(𝒳,𝒴)≥0 forall (𝒳,𝒴)∈Δ.(18)

Here, ℬ=D𝓏ℒ(𝒳0,0) represents an M-matrix and Δ denotes the feasible region of the model. The following lemma can be derived, based on the work of Castillo-Chavez et al. [31].

Lemma 1: If conditions G1 and G2 are met, the point E0=(𝒳0,0) is globally asymptotically stable (GAS) when ℛ0<1.

Now, we proceed to confirm the following theorem.

Theorem 7: The DFE point, E0, is GAS whenever the number ℛ0<1.

Proof: We define 𝒳=(Sh,Sv) to represent the number of susceptible individuals, comprising humans (Sh) and mosquitoes (Sv). 𝒴=(Eh,Ih,Ah,Ev,Iv) to represent the number of individuals in various infection states, including exposed humans (Eh), clinically infected humans (Ih), asymptomatic infected humans (Ah), exposed mosquitoes (Ev) and infected mosquitoes (Iv). The disease-free point, E0, is represented by E0=(𝒳0,0), where 𝒳0 denotes the number of susceptible individuals (both mosquitoes and humans) and zero represents the absence of infected individuals in all compartments. So

d𝒳dt=𝒢(𝒳,Y)=[Λh−σbShIvNh−μhShΛv−σcSvIhNh−μvSv].(19)

If Sh=Sh0 and Sv=Sv0, then 𝒢(𝒳,0)=0, i.e.,

d𝒳dt=[Λh−μhShΛv−μvSv].(20)

As t→∞, 𝒳→𝒳0. Therefore, 𝒳0=(Sh0,Sv0) is globally asymptotically stable.

Now,

ℬ𝒴−ℒ¯(𝒳,𝒴)=[−k1000σbSh0Nhβ1−k2000β20−k3000σcSv0Nh0−k40000νv−μv][EhIhAhEvIv]−[Υ0Π00],(21)

where, Υ=σbIvNh(Sh0−Sh) and Π=σcIhNh(Sv0−Sv), and

ℬ=[−k1000σbSh0Nhβ1−k2000β20−k3000σcSv0Nh0−k40000νv−μv],𝒴=[EhIhAhEvIv],ℒ¯(𝒳,𝒴)=[Υ0Π00].

The matrix ℬ is an M-matrix. As the equilibrium point is disease-free, we have Sh≤Sh0 and Sv≤Sv0, which implies that the matrix ℒ¯(𝒳,𝒴)≥0. Therefore, the disease-free equilibrium point, E0, is globally asymptotically stable. □

The following theorem illustrates the global stability of model (3) at the endemic equilibrium point E1.

Theorem 8: For disease model (3), if ℛ0>1, the epidemic equilibrium point E1 is globally asymptotically stable.

Proof: A Volterra-type Lyapunov function is a powerful tool for establishing global stability. The function is defined as:

L(Ψ)=[Sh−Sh1−Sh1log⁡ShSh1]+[Eh−Eh1−Eh1log⁡EhEh1]+[Ih−Ih1−Ih1log⁡IhIh1]+[Ah−Ah1−Ah1log⁡AhAh1]+[Rh−Rh1−Rh1log⁡RhRh1]+[Sv−Sv1−Sv1log⁡SvSv1]+[Ev−Ev1−Ev1log⁡EvEv1]+[Iv−Iv1−Iv1log⁡IvIv1].

The equilibrium point E1=(Sh1,Eh1,Ih1,Ah1,Rh1,Sv1,Ev1,Iv1) represents a state of equilibrium within the epidemic model. In terms of time t, we determine the following derivative of the Volterra-type Lyapunov function.

dLdt=[dShdt−Sh1ShdShdt]+[dEhdt−Eh1EhdEhdt]+[dIhdt−Ih1IhdIhdt]+[dAhdt−Ah1AhdAhdt]+[dRhdt−Rh1RhdRhdt]+[dSvdt−Sv1SvdSvdt]+[dEvdt−Ev1EvdEvdt]+[dIvdt−Iv1IvdIvdt],

or

dLdt=[Sh−Sh1Sh]dShdt+[Eh−Eh1Eh]dEhdt+[Ih−Ih1Ih]dIhdt+[Ah−Ah1Ah]dAhdt+[Rh−Rh1Rh]dRhdt+[Sv−Sv1Sv]dSvdt+[Ev−Ev1Ev]dEvdt+[Iv−Iv1Iv]dIvdt.

We substitute the derivatives in the equation with the right-hand side of the ODEs of the model (3), resulting in the following expression.

dLdt=[Sh−Sh1Sh][Λh−σbIvShNh−μhSh]+[Eh−Eh1Eh][σbIvShNh−(β1+β2+μh)Eh]+[Ih−Ih1Ih]×[β1Eh−(θ+μh+δh)Ih]+[Ah−Ah1Ah][β2Eh−(μh+γ)Ah]+[Rh−Rh1Rh][θIh+γAh−μhRh]+[Sv−Sv1Sv][Λv−σcIhSvNh−μvSv]+[Ev−Ev1Ev][σcIhSvNh−(νv+μv)Ev]+[Iv−Iv1Iv][νvEv−μvIv].