EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 3, Article Number 6522
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Stability Analysis and Numerical Simulation of Fractional-Order
Models for Wolbachia Transmission in Aedes aegypti Mosquitoes

with Seasonal Effects

Mohammad Yar1,2,∗, Shumaila Javeed1, Tanveer Abbas Khan1

1 Department of Mathematics, COMSATS University Islamabad, Park Road, Chak Shahzad
Islamabad, 45550, Pakistan
2 Department of Mathematics, Kabul Polytechnic University, Kabul 1001, Afghanistan

Abstract. Mosquito-borne diseases have historically impacted people and remain a global problem.
Wolbachia represents an innovative vector management strategy capable of diminishing mosquito pop-
ulations and reducing the threat of mosquito-borne diseases. In order to provide optimal and efficient
control, Wolbachia should be released at each step of the mosquito life cycle. Employing fractional
calculus enables the capture of memory effects and hereditary characteristics of this process. In this
study, we develop and analyze a fractional-order mathematical model to investigate Wolbachia trans-
mission dynamics, which accounts for imperfect maternal transmission and infection loss. The model
considers Wolbachia-infected and uninfected subpopulations of Aedes aegypti mosquitoes, assuming
equal numbers of adult males and females. We establish the positivity and boundedness of solutions
for non-negative initial conditions. The invasive reproduction number R0w|w̄ has been found to de-
termine whether the Wolbachia infection spreads. We consider two fractional-order models: one that
ignores the impacts of seasons on the mosquito populations and another that incorporates these ef-
fects. Furthermore, the stability of the model is analyzed using Lyapunov functions and Ulam-Hyers
stability theories. The models are numerically solved using the Adams-Bashforth-Moulton technique,
demonstrating the effects of model parameters and fractional-order values. Based on the results, we
find that the fractional order α = 0.5 is optimal, and the corresponding conditions could be applied
in real-world experiments to increase the population of Wolbachia-infected mosquitoes. These results
highlight the importance of fractional-order modeling.

2020 Mathematics Subject Classifications: 34A08, 92D30, 34D20

Key Words and Phrases: Wolbachia, Aedes aegypti, Fractional calculus, Lyapunov function, Ulam-
Hyers stability, Seasonal effects, Adams-Bashforth-Moulton method

1. Introduction

Mosquito-borne diseases are predominantly transmitted by female mosquitoes while feed-
ing on the blood of living organisms, including humans, animals, and birds. A female mosquito
infected with a parasite, virus, or bacterium can transmit these pathogens to humans [1].
Dengue fever is widespread globally, with approximately 3.9 billion individuals at risk and
390 million new cases annually [2]. Various strategies exist for managing Aedes aegypti

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6522

Email addresses: myar@kpu.edu.af (M. Yar),
shumailajaveed@comsats.edu.pk (S.Javeed), tanveerabbas947@gmail.com (T.A.Khan)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 2 of 23

mosquitoes, including insecticide spraying, sterile insect techniques, incompatible insect tech-
niques, combined approaches, and genetic engineering [3, 4].

Integrated vector management is a promising control strategy currently under investiga-
tion [5]. Substituting the Aedes aegypti mosquito population with alternative vector agents
that cannot transmit viruses, rather than only preventing human-vector interactions, has
shown a decrease in dengue distribution, and this method seems effective in managing other
mosquito-borne diseases like Zika, malaria, West Nile virus, and Chikungunya [6–8]. Wol-
bachia is a genus of bacteria that reduces mosquitoes’ capacity for transmitting viruses
[9]. Wolbachia-infected (WI) mosquitoes have less ability to transmit diseases compared
to Wolbachia-uninfected (WU) mosquitoes. Moreover, Wolbachia infection frequently causes
cytoplasmic incompatibility, resulting in premature embryonic mortality when infected males
mate with uninfected females. When a Wolbachia-infected mosquito bites a virus-infected
human, it may become infectious, but it is unable to transmit the virus to an uninfected
human. The Wolbachia method would significantly diminish the mosquito population density
and facilitate the elimination of mosquito-borne diseases [10, 11].

On the other hand, imperfect maternal transmission can also prevent WI mosquitoes from
effectively controlling disease transmission [12–14]. To implement the method extensively,
comprehensive details about introducing Wolbachia into wild mosquitoes are required. The
application of mathematical modeling in the decision-making process has facilitated the use
of numerous widely recognized control policies. It plays an important role in assessing the
impacts of various factors on the dynamics of infectious diseases. Several mathematical models
have been established on the introduction of Wolbachia in mosquitoes, each illustrating the
conditions that facilitate the rapid expansion of WI mosquitoes [15–21].

Additionally, fractional calculus (FC) is an emerging field of mathematical analysis that
is developing rapidly. Numerous scientists have established several models of real processes
using FC. A realistic representation of a physical phenomenon requires not just the current
time but also the history of the preceding time, which FC can provide [22, 23]. In other words,
fractional derivative definitions naturally allow for the incorporation of memory and heredi-
tary characteristics in the modeling of diverse materials and processes [24]. Several models of
infectious diseases based on FC have already been discussed in the literature [25–28].

Recently, researchers have developed mathematical models to investigate the key fac-
tors that influence Wolbachia’s effectiveness in managing viral infections. Wan and Xu [29]
created a mathematical model using optimal control theory to understand how Wolbachia
infections and dengue spread work together, aiming to find the best ways to release young
mosquitoes that die quickly. They found that increasing larval mortality prevented Wol-
bachia establishment and dengue treatments efficiency. Dianavinnarasi et al. [30] formu-
lated a fractional-order Wolbachia invasion model, utilising impulsive control approaches to
maintain infection in Aedes aegypti mosquitoes. The results showed that quick releases of
Wolbachia can lower the number of local mosquitoes while keeping the infected ones alive. A
fractional-order dengue transmission model by Vijayalakshmi and Ariyanatchi [31] looked at
how both Wolbachia-infected and uninfected mosquitoes behave. They found that Wolbachia
reduces dengue transmission and improves mosquito control. Dianavinnarasi et al. [32] used
a fractional-order mathematical model to study how different Wolbachia strains spread in
Aedes aegypti mosquitoes, with the goal of finding the best strain for long-term use. They
demonstrated that Wolbachia is an efficacious approach for managing mosquito-borne ill-



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 3 of 23

nesses. Ufuktepe [33] created a system that looks at how Wolbachia spreads in mosquito
groups over time, considering how the presence of fewer wild insects can affect this spread.
They examined the rivalry between released mosquitoes and indigenous mosquitoes. The in-
vestigation revealed that Wolbachia-infected mosquitoes can surpass native populations under
certain conditions, resulting in possible mosquito control.

Additionally, it is expected that climate change will alter the distribution and seasonal
behaviors of mosquitoes, significantly affecting the persistence and seasonality of vector-borne
diseases [34]. This is expected to become increasingly important in the coming years as envi-
ronmental conditions are anticipated to become more variable.

The novelty of this study is in the use of a fractional-order mathematical model to inves-
tigate the dynamics of Wolbachia invasion and the Effects of changes in the seasons on the
development and spread of Wolbachia within the Aedes aegypti mosquito population. The
model is constructed using the Caputo fractional derivative, which provides benefits including
the ability to capture memory-dependent effects and enhance model accuracy. We incorporate
the consequences of incomplete maternal transmission and the absence of Wolbachia infection.
A mosquito population comprising two types, WU and WI, is analyzed. The model is evalu-
ated for positivity and boundedness to ensure biological feasibility. The invasive reproduction
number is determined to evaluate the efficacy of Wolbachia-infected mosquitoes in dengue
control. The stability of the equilibrium points is examined using Lyapunov’s direct method.
Numerical simulations, performed using the fractional Adams–Bashforth–Moulton technique,
illustrate the influence of model parameters and fractional-order values on the spread and
control of Wolbachia infection. Furthermore, we investigate the impact of seasonal fluctua-
tions on the introduction and transmission of Wolbachia within the Aedes aegypti mosquito
population, emphasizing their effects across various life stages of the mosquito. The paper
is structured as follows: An introductory section sets the context. Section 2 introduces the
model formulation, fundamental definitions, and essential properties of the fractional oper-
ator in the Caputo sense. Section 3 discusses the positivity and boundedness of solutions.
In Section 4, we examine the model by identifying the disease-free and endemic equilibrium
points and derive the expression for the control reproduction number. Section 5 focuses on
the stability analysis of equilibrium points and provides some graphical representations of the
model’s dynamics. Section 6 presents the Ulam-Hyers stability of the model, and Section 7
examines the seasonal effect of the proposed model. Finally, Section 8 concludes the study
with a brief discussion.

2. Description of model framework

We develop an innovative fractional-order mathematical model that represents the
life stages of Aedes aegypti mosquitoes. The Aedes aegypti mosquito population is categorized
into two subpopulations: WI (w) and WU (w̄). Assuming an equal distribution of adult
male and female mosquitoes, we designate adult female mosquitoes free of Wolbachia as Fw̄
and adult WI female mosquitoes as Fw [18, 35, 36]. The fractional mathematical model
for Wolbachia invasion, considering imperfect maternal transmission and infection loss, is as



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 4 of 23

follows: 

CDαQw̄ =

[
ϕw̄F

2
w̄ + ρ1ϕwF

2
w + ρ2ϕwFwFw̄

Fw̄ + Fw

](
1− Q

K

)
− (µa + ψ)Qw̄,

CDαQw =

[
(1− ρ1)ϕwF

2
w + (1− ρ2)ϕwFwFw̄
Fw̄ + Fw

](
1− Q

K

)
− (µa + ψ)Qw,

CDαFw̄ =
ψ

2
Qw̄ + σFw − µw̄Fw̄,

CDαFw =
ψ

2
Qw − σFw − µwFw.

(1)

Here, Qw̄(0) = 200000, Qw(0) = 500000, Fw̄(0) = 900000, and Fw(0) = 600000 are the initial
conditions. In the above equation, CDα denotes the Caputo fractional derivative where α is
the order of derivative, 0 < α ≤ 1. K stands for the aquatic phase’s carrying capacity, ϕw̄
denotes the egg-laying rate for Wolbachia-free mosquitoes and, ϕw is the egg-laying rate for
mosquitoes with Wolbachia infection. The ρ1 is the fraction of eggs that are Wolbachia free
when adult female and male mosquitoes with Wolbachia infection mate, ρ2 is the fraction of
eggs that are Wolbachia free when adult Wolbachia free males and females with Wolbachia
infection mate, σ denotes loss of Wolbachia infection, ψ is the maturation rate, µa represents
the aquatic death rate, µw̄ is the death rate of Wolbachia-free mosquitoes, and µw denotes the
death rate of mosquitoes with Wolbachia infection. The following definitions and properties
of fractional operators will be useful throughout our work.

Definition 1. The fractional integral for a function f : (t0,∞) → R is defined by

Iα
t f(t) :=

1

Γ(q)

∫ t

t0

(t− x)α−1f(x)dx,

where α is the order of integral with 0 < α ⩽ 1, n ∈ N and Γ(.) denotes the gamma function.

Definition 2. The Caputo fractional derivative of order α with n− 1 < α ⩽ n is defined as:

CDα
t f(t) :=

1

Γ(n− α)

∫ t

0

fn(x)

(t− x)α−n+1
dx.

Theorem 1. Let Re(α) > 0, n = [Re(v)] + 1. Then

Iα
t D

α
t f(t) = f(x)−

n∑
k=1

(Dk
t f)

k!
xk. (2)

Lemma 1. [37] Let x(t) be a positive real, continuous and differentiable function. Then for
any t ≥ t0 and 0 < α ⩽ 1, and x∗ > 0, we have

CDα
t

[
x(t)− x∗ − x∗ ln

(
x(t)

x∗

)]
≤
(
1− x∗

x(t)

)
Dα

t x(t). (3)

We select the Caputo derivative; however, various fractional derivatives, like the
Riemann-Liouville derivative and the Caputo-Fabrizio derivative, are also relevant. A key
advantage of adopting the Caputo fractional derivative is its ability to facilitate the solution



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 5 of 23

of problems with standard initial and boundary conditions. Real-world problems seem to be
simple to physically understand using the Caputo derivative [19].

Description Parameter Estimate[Range] Unit References
Aquatic stage’s carrying capacity K 106[104, 108] Aquatic mosquito [38]
Egg laying rate for WU mosquitoes ϕw̄ 13[12-18] Eggs per day [39, 40]
Egg laying rate for WI mosquitoes ϕw 11[8-12] Eggs per day [12, 39]
Fraction of eggs which is WU when
adult female and male mosquitoes with
WI mate

ρ1 0.05 [0-0.11] Dimensionless [12]

Fraction of eggs which is WU when
adult WU male and WI female
mosquitoes mate

ρ2 0.05[0-0.1] Dimensionless [12]

Per capita reduction of Wolbachia in-
fection

σ 0.04[0-0.1] Per day [38]

Fraction of eggs that develop as fe-
males

a 0.5[0.34-0.6] Dimensionless [36, 41]

Per capita maturation rate ψ 0.11 [0.1-0.12] Per day [12, 39]
Per capita aquatic death rate µa 0.02 Per day [18]
Per-capita mortality rate of WU
mosquitoes

µw̄ 0.061[0.02-0.09] Per day [12, 39]

Per-capita mortality rate of WI
mosquitoes

µw 0.068[0.03-0.14] Per day [12, 42]

Table 1: Overview of model parameters and their corresponding descriptions

3. Positivity and boundedness of solutions

Establishing the positivity and boundedness of all solutions of the model will be
significant because the problem is a biological model, which means that all of the solutions
of the system (1) should remain non-negative and bounded ∀t ≥ 0. First, from a biological
perspective, the initial data Qw̄(0), Qw(0), Fw̄(0), and Fw(0) must be larger than or equal to
zero. For the positivity, consider system (1) in the form:

CDαQw̄|Qw̄=0 =

[
ϕw̄F

2
w̄ + ρ2ϕwFwFw̄
Fw̄ + Fw

](
1− Q

K

)
⩾ 0,

CDαQw|Qw=0 =

[
ϕwF

2
w + (1− ρ2)ϕwFwFw̄

Fw̄ + Fw

](
1− Q

K

)
⩾ 0,

CDαFw̄|Fw̄=0 =
ψ

2
Qw̄ ⩾ 0,

and

CDαFw|Fw=0 =
ψ

2
Qw ⩾ 0.

According to Lemmas 5 and 6 in [43], all the solutions of the model remain non-negative. Next,
we show the boundedness of the solutions. Following [44], we analyze the entire population
and define the total population as:

M(t) = Qw̄(t) +Qw(t) + Fw̄(t) + Fw(t).

By adding all equations of the system (1), we will have

CDαM(t) =CDαQw̄(t) +
CDαQw(t) +

CDαFw̄(t) +
CDαFw(t)



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 6 of 23

=

[
ϕw̄F

2
w̄ + ϕwFwFw̄ + ϕwF

2
w

Fw̄ + Fw

](
1− Q

K

)
− µa(Qw̄ +Qw)

− ψ

2
(Qw̄ +Qw)− µw̄Fw̄ − µwFw. (4)

Since Qw̄ < K, Qw < K, it follows from model(1), Fw̄ ⩽ ψK
2µ1

and Fw ⩽ ψK
2µ1

, where µ1 =
min(µw̄, µw, µa). Thus, Eq. (4) becomes

CDαM(t) ⩽ Λ3 − µ1M(t),

where, Λ3 =
ψK(ϕw̄+2ϕw)

4µ1
. So we deduce

M(t) ⩽M(0)Eα(−µ1 · tα) +
Λ3

µ1
(1− Eα(−µ1 · tα)).

Where Eα is the Mittag-Leffler function of parameter α. Since 0 < Eα(−µ1 · tα) ⩽ 1 and
1 − Eα(−µ1 · tα) ⩽ 1. We obtain M(t) ⩽ M(0) + Λ3

µ1
. Therefore, the solutions of model (1)

are uniformly bounded.

4. Equilibrium points

The equilibrium points of the non-integer order system (1) are obtained by solving
the corresponding nonlinear algebraic system under the conditions CDαQw̄ = 0, CDαQw = 0,
CDαFw̄ = 0 and CDαFw = 0. The model (1), with ρ1 = 0 and σ = 0, yields four equilibrium
points: E1 = (0, 0, 0, 0) - demonstrating the absence of mosquitoes; E2 = (Q∗

w̄, 0, F
∗
w̄, 0) - it il-

lustrates the dominance of WU mosquitoes; E3 = (0, Q∗
w, 0, F

∗
w) - demonstrating the existence

of mosquitoes with Wolbachia infection and E4 = (Q∗
w̄, Q

∗
w, F

∗
w̄, F

∗
w) - showing coexistence of

both WU and WI mosquitoes. Determining the nature of stability points is important for
controlling arboviral infections transmitted by Aedes aegypti mosquitoes.

4.1. No mosquitoes

The equilibrium point E1 is trivial but uninteresting since it is biologically un-
realistic. However, by looking at a particular example in which there is no contact between
WI and WU mosquitoes can offer additional insight into the dynamics of this steady-state
solution. We derived

R0w̄ =
ϕw̄ψ

2µw̄(µα + ψ)
, (5)

and

R0w =
ϕwψ

2µw(µa + ψ)
, (6)

which represent threshold values determining whether each group will survive or go extinct in
the absence of contact. In the absence of contact between infected and uninfected mosquitoes,
the thresholds in Eqs. (5) and (6) are obtained from the stability criteria of the appropriate
Jacobian matrix under conditions of no contact among the two groups of mosquitoes. In other



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 7 of 23

words, each group is autonomous from the others. Equivalent methods for the dynamics
explicitly including the male mosquito states were presented in [18]. Thus, for the model
(1), both populations are extinguished when R0w < 1 and R0w̄ < 1, as the reproductive
parameters are insufficient to maintain the populations. Furthermore, given that the solutions
remain non-negative for non-negative initial conditions, the solutions tend to the no-mosquito
equilibrium point. However, and except for the biological implications of using insecticides,
applying insecticides and destroying breeding sites have been an effective method in reducing
mosquito populations.

4.2. Wolbachia uninfected mosquitoes only

The equilibrium point of WU mosquitoes is specified as E2 = (Q∗
w̄, 0, F

∗
w̄, 0) where

Q∗
w̄ = K

(
1− 1

R0w̄

)
, (7)

and

F ∗
w̄ =

ψK

2µw̄

(
1− 1

R0w̄

)
. (8)

Hence, for this equilibrium point, it is necessary that R0w̄ > 1; otherwise, WU mosquitoes
will not persist. We now establish the invasive reproduction number R0w|w̄ using the next-
generation matrix method. WI populations can be divided into two categories: one char-
acterized by the emergence rate of new WI mosquitoes (F) and the other transition rates,
including mortality and development into adult mosquitoes infected with Wolbachia (V) [12].

F =

(
(ϕwF

2
w+(1−ρ2)ϕwFwFw̄

Fw̄+Fw
)(1− Q

K )

0

)
,

V =

(
(µa + ψ)Qw

−ψQw

2 + µwFw

)
.

Furthermore, we introduce the matrices

Fij =
∂Fi
∂xj

|E2 , (9)

and

Vij =
∂Vi
∂xj

|E2 , (10)

where xj show the infected population Qw and Fw. Thus,

F =

(
0

ϕw(K−Q∗
w̄)(1−ρ2)
K

0 0

)
,

V =

(
µa + ψ 0

−ψ
2 µw

)
,



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 8 of 23

the next generation matrix is

FV −1 =

(
ψϕw(K−Q∗

w̄)(1−ρ2)
2(µa+ψ)µwK

ϕw(K−Q∗
w̄)(1−ρ2)

µwK

0 0

)
.

Hence, the invasive reproduction number is

R0w|w̄ = λ(FV −1) =
ϕwµw̄(1− ρ2)

µwK
=
R0w(1− ρ2)

R0w̄
. (11)

Here, λ(FV −1) denotes the spectral radius of FV −1. The number R0w|w̄ represents the
expected number of infected offspring produced by a single WI mosquito introduced into a
population of WU mosquitoes at equilibrium. Biologically, this metric quantifies the potential
for Wolbachia to invade and establish in a WU dominated population, under the assumption
that infected mosquitoes do not lose their infection. The factor (1−ρ2) in this equation reflects
the effect of imperfect maternal transmission: some offspring from WI females mating with
WU males are uninfected, which diminishes the chances of successful invasion. This factor
thus demonstrates how the ratio of aquatic-stage mosquitoes that are WI affects the likelihood
of WI mosquitoes displacing WU mosquitoes when WU males mate with WI females. If
R0w|w̄ < 1, the infection is unable to establish itself, even with cytoplasmic incompatibility
present. Therefore, cytoplasmic incompatibility by itself is not enough for successful invasion;
the reliability of maternal transmission is also critical.

4.3. Wolbachia infected mosquitoes-only

In this instance, the equilibrium point is E3 = (0, Q∗
w, 0, F

∗
w) where,

Q∗
w = K

(
1− 1

R0w

)
, (12)

and

F ∗
w =

ψK

2µw

(
1− 1

R0w

)
. (13)

4.4. Coexistence of Wolbachia infected and Wolbachia uninfected mosquitoes

The presence of both WI andWUmosquitoes within the Aedes aegypti population
is a significant situation. In this case, it is desirable for the majority of the mosquito population
to be Wolbachia-infected, as this can reduce the transmission of arboviral diseases. The
coexistence equilibrium point of system (1) is defined as

E4 = (d1F
∗
w̄, d2F

∗
w, d3F

∗
w, F

∗
w), (14)

where F ∗
w = Kψ

2(µw̄d3+µw)

[
R0w(1+(1−ρ2)d3)−(1+d3)

R0w(1+(1−ρ2)d3)

]
, d1 =

2µw̄
ψ , d2 =

2µw
ψ and, d3 =

[
R0w|w̄(µw̄−ρ2µw)

µw̄(1−ρ2)(1−R0w|w̄)

]
.



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 9 of 23

5. Global stability of equilibrium points

Lyapunov’s direct method, also known as the second Lyapunov method, gives
us an effective way to analyze the global behavior of a system without explicitly solving
it. Classical Lyapunov functions in mathematical biology are linear combinations of linear,
common quadratic, and Volterra-type functions [27, 37, 45]. We now examine the global
asymptotic stability of the WU equilibrium point E2, the WI mosquito equilibrium point E3,
and the combined WU and WI equilibrium point E4 of system (1). The equilibrium point for
adult female WI mosquitoes derived from Eq. (11) can be written as

F ∗
w =

ψK

2µw

(
1− 1− ρ2

R0w|w̄R0w̄

)
. (15)

The above expression clearly shows that when R0w|w̄ < 1, the WI mosquitoes-only equilibrium
may exist. The conventional infectious diseases modeling work indicates a backward bifur-
cation in the presence of endemic equilibria for R0w|w̄ < 1 [46]. This equilibrium is unstable

when R0w|w̄ < 1−ρ2
R0w̄

(=⇒R0w < 1) and locally asymptotically stable even when R0w|w̄ < 1,
since R0w̄ > 1 and µw̄ > ρ2µw. Whenever R0w|w̄ < 1, both equilibrium points E2 and E3 are
locally asymptotically stable since R0w̄ > 1 for E2, R0w > 1 and µw̄ > ρ2µw for E3.

Theorem 2. The equilibrium point E2 = (Q∗
w̄, 0, F

∗
w̄, 0) is globally asymptotically stable when-

ever R0w|w̄ > 1 and R0w̄ > 1.

Proof. Let us consider the following Lyapunov function:

V1(Qw̄(t), Fw̄(t)) =
ψϕw̄

2µw̄(µa + ψ)

∫ Qw̄

Q∗
w̄

(
1− Q∗

w̄

x

)
dx+

1

µw̄

∫ Fw̄

F ∗
w̄

(
1− F ∗

w̄

x

)
dx

=
ψϕw̄

2µw̄(µa + ψ)

[
Qw̄ −Q∗

w̄ −Q∗
w̄ ln

(
Q∗
w̄

Qw̄

)]
+

1

µw̄

[
Fw̄ − F ∗

w̄ − F ∗
w̄ ln

(
Fw̄
Fw̄

)]
.

(16)
The V1 is a well-defined, continuous, and positive definite function. Utilizing the property of
fractional derivatives as defined in the lemma 1, we compute the derivative of Eq. (16) with
respect to time along the solution of system (1) and we demonstrate that the corresponding
fractional derivative is negative definite, using Lemma 3.1 in [45], we have

CDαV1(Qw̄(t), Fw̄(t)) ≤
ψϕw̄

2µw̄(µa + ψ)

(
1− Q∗

w

Qw̄

)
CDαQw̄(t) +

1

µw̄

(
1− F ∗

w̄

Fw̄

)
CDαFw̄(t).

(17)
Substituting the expression for the model (1) we have,

ψϕw̄
2µw̄(µa + ψ)

(
1− Q∗

w̄

Qw̄

)
CDαQw̄(t) =

ψϕw̄
2µw̄(µa + ψ)

(
1− Q∗

w̄

Qw̄

)
[
Λ1

(
1− Q

K

)
− (µa + ψ)Qw̄

]
=

ψϕw̄Λ1

2µw̄(µa + ψ)

(
1− Q∗

w̄

Qw̄

)(
1− Q

K

)
− ψQw̄

2µw̄
+
ψQ∗

w̄

2µw̄
, (18)



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 10 of 23

where Λ1 =
(
ϕw̄F 2

w̄+ρ2ϕwFwFw̄

Fw̄+Fw

)
and,

1

µw̄

(
1− F ∗

w̄

Fw̄

)
CDαFw̄(t) =

1

µw̄

(
1− F ∗

w̄

Fw̄

)(
ψ

2
Qw̄ − µw̄Fw̄

)
=
ψQw̄
2µw̄

− ψQw̄F
∗
w̄

2µw̄Fw̄
− Fw̄ + F ∗

w̄. (19)

Adding Eqs.(18) and (19) yields

CDαV1(Qw̄(t), Fw̄(t)) ≤
ψϕw̄Λ1

2µw̄(µa + ψ)

(
1− Q∗

w̄

Qw̄

)(
1− Q

K

)
+
ψQw̄
2µw̄

− ψQw̄F
∗
w̄

2µw̄Fw̄
− Fw̄ + F ∗

w̄.

(20)

Using the relations at the steady state and some algebraic manipulations, we get

CDαV1(Qw̄(t), Fw̄(t)) ≤
ψϕw̄Λ1

2µw̄(µa + ψ)

(
1− Q∗

w̄

Qw̄

)(
1− Q

K

)
+ F ∗

w̄

(
2− Qw̄F

∗
w̄

Q∗
w̄Fw̄

− Q∗
w̄Fw̄

Qw̄F ∗
w̄

)
− Fw̄

(
1− Q∗

w̄

Qw̄

)
.

(21)

Thus,

CDαV1(Qw̄(t), Fw̄(t)) ≤Fw̄
(
1− Q∗

w̄

Qw̄

)[
Row̄(Fw̄ϕw̄ + ρ2ϕwFw)

Fw̄ + Fw

(
1− Q

K

)
− 1

]
+ F ∗

w̄

(
2− Qw̄F

∗
w̄

Q∗
w̄Fw̄

− Q∗
w̄Fw̄

Qw̄F ∗
w̄

)
.

(22)

Eq. (22) can be reformulated as

CDαV1(Qw̄(t), Fw̄(t)) ≤F ∗
w̄

(
F ∗
w̄

Fw̄
+
Qw̄F

∗
w̄

Q∗
w̄Fw̄

− 2

)[
Row̄(Fw̄ϕw̄ + ρ2ϕwFw)

Fw̄ + Fw

(
1− Q

K

)
− 1

]
+ F ∗

w̄

(
2− Qw̄F

∗
w̄

Q∗
w̄Fw̄

− Q∗
w̄Fw̄

Qw̄F ∗
w̄

)[
Row̄(Fw̄ϕw̄ + ρ2ϕwFw)

Fw̄ + Fw

(
1− Q

K

)]
.

(23)
The second component on the right side of Eq. (23) is less than or equal to zero due to(

2− Qw̄F
∗
w̄

Q∗
w̄Fw̄

− Q∗
w̄Fw̄

Qw̄F ∗
w̄

)
≤ 0.

Additionally, Q∗
w̄ ≤ Qw̄ ≤ K results in the first term either less than or equal to zero. When

0 ≤ Qw̄ ≤ Q∗
w̄, Eq. (23) is less than zero because

(
1− Q∗

w̄
Qw̄

)
< 0. Hence, CDαV1(Qw̄(t), Fw̄(t)) ≤

0. So, the function V1(Qw̄(t), Fw̄(t)) is negative definite for 0 < α < 1. Then the equilibrium
state E2 is globally asymptotically stable.

Theorem 3. The equilibrium point E3 = (0, Q∗
w, 0, F

∗
w) is globally asymptotically stable when-

ever R0w|w̄ > 1, R0w > 1 and µw̄ > ρ2µw.



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 11 of 23

Proof. When R0w|w̄ > 1 then R0w > R0w̄
1−ρ2 and this implies R0w > R0w̄. We

consider the following Volterra-type Lyapunov function as

V2(Qw(t), Fw(t)) =
ψ

2µw(µa + ψ)

∫ Qw

Q∗
w

(
1− Q∗

w

x

)
dx+

1

µw

∫ Fw

F ∗
w

(
1− F ∗

w

x

)
dx. (24)

By calculating the α-order derivative of V2(Qw(t), Fw(t)), one has,

CDαV2(Qw(t), Fw(t)) ≤
ψ

2µw(µa + ψ)

(
1− Q∗

w

Qw

)
CDαQw(t) +

1

µw

(
1− F ∗

w

Fw

)
CDαFw(t).

(25)
Substituting the expression for the model (1) we have,

ψ

2µw(µa + ψ)

(
1− Q∗

w

Qw

)
CDαQw(t) =

ψ

2µw(µa + ψ)

(
1− Q∗

w

Qw

)
[
Λ2

(
1− Q

K

)
− (µa + ψ)Qw

]
,

(26)

where Λ2 =
[
(1−ρ1)ϕwF 2

w+(1−ρ2)ϕwFwFw̄

Fw̄+Fw

]
and

1

µw

(
1− F ∗

w

Fw

)
CDαFw(t) =

1

µw

(
1− F ∗

w

Fw

)(
ψ

2
Qw − µwFw

)
. (27)

Form Eq. (26),

ψ

2µw(µa + ψ)

(
1− Q∗

w

Qw

)
CDαQw(t) =

ψΛ2

2µw(µa + ψ)

(
1− Q∗

w

Qw

)(
1− Q

K

)
− ψQw

2µw
+
ψQ∗

w

2µw
,

(28)

and from Eq. (27),

1

µw

(
1− F ∗

w

Fw

)
CDαFw(t) =

ψQw
2µw

− ψQwF
∗
w

2µwFw
− Fw + F ∗

w. (29)

By adding Eqs. (28) and (29), and after rearrangements, we obtain

CDαV2(Qw(t), Fw(t)) ≤ Fw

(
1− Q∗

w

Qw

)[
R0w(Fw + (1− ρ2)Fw̄)

Fw + Fw̄

(
1− Q

K

)
− 1

]
+ F ∗

w

(
2− QwF

∗
w

Q∗
wFw

− Q∗
wFw

QwF ∗
w

)
.

(30)

Eq. (30), which can then be simplified to

CDαV2(Qw(t), Fw(t)) ≤ F ∗
w

(
F ∗
w

Fw
+
QwF

∗
w

Q∗
wFw

− 2

)[
R0w(Fw + (1− ρ2)Fw̄)

Fw + Fw̄

(
1− Q

K

)
− 1

]
+ F ∗

w

(
2− QwF

∗
w

Q∗
wFw

− Q∗
wFw

QwF ∗
w

)[
R0w(Fw + (1− ρ2)Fw̄)

Fw + Fw̄

(
1− Q

K

)]
.

(31)
Thus, CDαV2(Qw(t), Fw(t)) ≤ 0. So, the function V2(Qw(t), Fw(t)) is negative definite for
0 < α < 1. Then the equilibrium state E3 is globally asymptotically stable.



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 12 of 23

Theorem 4. The equilibrium point E4 = (d1F
∗
w̄, d2F

∗
w, d3F

∗
w, F

∗
w) of the non-integer order

system (1) is globally asymptotically stable if R0w|w̄ > 1, R0w > 1 and R0w̄ > 1.

Proof. The following positive definite Lyapunov function is used to analyze the
global asymptotic stability of the equilibrium point E4:

V3(t) = T1Q(Qw̄(t)) + T2Q(Qw(t)) + T3F (Fw̄(t)) + T4F (Fw(t)), (32)

where,

Q(Qw(t)) = Qw −Q∗
w −Q∗

w ln

(
Q∗
w

Qw

)
, Q(Qw̄(t)) = Qw̄ −Q∗

w̄ −Q∗
w̄ ln

(
Q∗
w̄

Qw̄

)
,

F (Fw̄(t)) = Fw̄ − F ∗
w̄ − F ∗

w̄ ln

(
F ∗
w̄

Fw̄

)
, F (Fw(t)) = Fw − F ∗

w − F ∗
w ln

(
F ∗
w

Fw

)
,

T1 =
ψϕw̄

2µw̄(µa + ψ)
, T2 =

ψ

2µw(µa + ψ)
, T3 =

1

µw̄
and T4 =

1

µw
.

Differentiating V3(t) with respect to time and by applying Lemma 1, one has

CDαV3(t) ≤ ψϕw̄
2µw̄(µa + ψ)

(
1− Q∗

w

Qw̄

)
CDαQw̄(t) +

ψ

2µw(µa + ψ)

(
1− Q∗

w

Qw

)
CDαQw(t)

+
1

µw̄

(
1− F ∗

w̄

Fw̄

)
CDαFw̄(t) +

1

µw

(
1− F ∗

w

Fw

)
CDαFw(t). (33)

Substituting the expressions for the derivatives into Eq. (33), and subsequently rearranging
and manipulating the terms, gives

CDαV3(t) ≤ F ∗
w̄

(
F ∗
w̄

Fw̄
+
Qw̄F

∗
w̄

Q∗
w̄Fw̄

− 2

)[
Row̄(Fw̄ϕw̄ + ρ2ϕwFw)

Fw̄ + Fw

(
1− Q

K

)
− 1

]
+F ∗

w̄

(
2− Qw̄F

∗
w̄

Q∗
w̄Fw̄

− Q∗
w̄Fw̄

Qw̄F ∗
w̄

)[
Row̄(Fw̄ϕw̄ + ρ2ϕwFw)

Fw̄ + Fw

(
1− Q

K

)]
+F ∗

w

(
F ∗
w

Fw
+
QwF

∗
w

Q∗
wFw

− 2

)[
R0w(Fw + (1− ρ2)Fw̄)

Fw + Fw̄

(
1− Q

K

)
− 1

]
+F ∗

w

(
2− QwF

∗
w

Q∗
wFw

− Q∗
wFw

QwF ∗
w

)[
R0w(Fw + (1− ρ2)Fw̄)

Fw + Fw̄

(
1− Q

K

)]
. (34)

Hence, we conclude by saying that since V3(t) is negative definite for 0 < α < 1, then the
equilibrium state E4 is globally asymptotically stable.



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 13 of 23

0 500 1000 1500 2000 2500 3000 3500 4000

Days

0

1

2

3

4

5

6

7

8

9

N
um

be
r o

f A
e.

ae
gy

pt
i m

os
qu

ito
es

×105 α  = 0.7

Aquatic WU

Aquatic WI

Adult WU

Adult WI

(a)

0 500 1000 1500 2000 2500 3000 3500 4000

Days

0

1

2

3

4

5

6

7

8

9

N
um

be
r o

f A
e.

ae
gy

pt
i m

os
qu

ito
es

×105 α  = 0.8

Aquatic WU

Aquatic WI

Adult WU

Adult WI

(b)

0 50 100 150 200 250 300 350 400

Days

0

1

2

3

4

5

6

7

8

9

N
um

be
r o

f A
e.

ae
gy

pt
i m

os
qu

ito
es

×105 α  = 0.9

Aquatic WU

Aquatic WI

Adult WU

Adult WI

(c)

0 50 100 150 200 250 300 350 400

Days

0

1

2

3

4

5

6

7

8

9

N
um

be
r o

f A
e.

ae
gy

pt
i m

os
qu

ito
es

×105 α  = 1

Aquatic WU

Aquatic WI

Adult WU

Adult WI

(d)

Figure 1: This figure illustrates the numerical simulation of aquatic-stage WI and WU
mosquitoes, as well as adult WI and WU mosquitoes, across various fractional order values of
α. The plots illustrate the impact of the memory effect, governed by the fractional order, on
the temporal dynamics of these populations as α approaches 1. We set ϕw = 0.01, ϕw̄ = 0.01,
K = 2×106, Qw̄(0) = 2×105, Qw(0) = 5×105, Fw̄(0) = 9×105, and Fw(0) = 6×105. Other
parameters used for these graphs simulations are provided in Table 1.

Hence, for R0w < 1 and R0w̄ < 1 both populations die out (see Figure 1) because
the populations can not be sustained in reproductive terms. Furthermore, the solutions con-
verge to the non-mosquito equilibrium point, as they remain non-negative with non-negative
initial conditions. All populations take much time to die out as the order of α decreases (see
Figure 1).



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 14 of 23

0 50 100 150 200 250 300 350 400

Days

2

3

4

5

6

7

8

9

10

11

12

N
um

be
r 

of
 A

e.
ae

gy
pt

i m
os

qu
ito

es

×105
α  = 0.5

Aquatic WU

Aquatic WI

Adult WU

Adult WI

(a)

0 50 100 150 200 250 300 350 400

Days

2

3

4

5

6

7

8

9

10

11

12

N
um

be
r 

of
 A

e.
ae

gy
pt

i m
os

qu
ito

es

×105
α  = 0.7

Aquatic WU

Aquatic WI

Adult WU

Adult WI

(b)

0 50 100 150 200 250 300 350 400

Days

2

3

4

5

6

7

8

9

10

11

12

N
um

be
r 

of
 A

e.
ae

gy
pt

i m
os

qu
ito

es

×105
α  = 0.8

Aquatic WU

Aquatic WI

Adult WU

Adult WI

(c)

0 50 100 150 200 250 300 350 400

Days

2

3

4

5

6

7

8

9

10

11

12

N
um

be
r 

of
 A

e.
ae

gy
pt

i m
os

qu
ito

es

×105
α  = 0.9

Aquatic WU

Aquatic WI

Adult WU

Adult WI

(d)

0 50 100 150 200 250 300 350 400

Days

2

3

4

5

6

7

8

9

10

11

12

N
um

be
r 

of
 A

e.
ae

gy
pt

i m
os

qu
ito

es

×105
α  = 1

Aquatic WU

Aquatic WI

Adult WU

Adult WI

(e)

Figure 2: This figure presents the numerical simulation of the model under different values of
the fractional order α. In these graphs we set ϕw = 2, ϕw̄ = 1, K = 2× 106, Qw̄(0) = 2× 105,
Qw(0) = 5× 105, Fw̄(0) = 9× 105, and Fw(0) = 600000. Table 1 contains a list of additional
parameters that are used for the simulations of these graphs.

The α = 0.5 is the most suitable value for this model. The number of WI
mosquitoes rises as the value of α falls from 1, while the number of WU mosquitoes de-
creases slowly. The WI mosquito population continues to increase, but at a very slow rate



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 15 of 23

when the value of α falls below 0.5. At this point, however, the WU mosquito population
continues to decrease, as shown in Figure 2. WU mosquito populations started increasing
when α < 0.5.

6. Ulam-Hyers stability analysis

The stability theory is a fundamental aspect of the qualitative analysis of dif-
ferential equations. The concept of Ulam stability was introduced by Ulam [47]. in 1940.
He proposed to investigate the degree of approximation between the approximate solution
and the exact solution of the equations. In 1941, Hyers [48] addressed Ulam’s notion and
established the concept of Ulam-Hyers stability for equations. It addresses short-term per-
turbations and analyzes how minor alterations in the equation and initial conditions affect
the solution inside a local neighborhood. Furthermore, given that stability is necessary for
an approximate solution, we focus on the Ulam-Hyers stability for the model (1) by using the
technique of nonlinear functional analysis. We will first present some lemmas essential for
stability. We now reformulate the model (1) in the following manner.

CDαQw̄ = Θ1(t, Qw̄, Qw, Fw̄, Fw),
CDαQw = Θ2(t, Qw̄, Qw, Fw̄, Fw),
CDαFw̄ = Θ3(t, Qw̄, Qw, Fw̄, Fw),
CDαFw = Θ4(t, Qw̄, Qw, Fw̄, Fw).

(35)

Where

Θ1(t, Qw̄, Qw, Fw̄, Fw) =

[
ϕw̄F

2
w̄ + ρ1ϕwF

2
w + ρ2ϕwFwFw̄

Fw̄ + Fw

](
1− Q

K

)
− (µa + ψ)Qw̄,

Θ2(t, Qw̄, Qw, Fw̄, Fw) =

[
(1− ρ1)ϕwF

2
w + (1− ρ2)ϕwFwFw̄
Fw̄ + Fw

](
1− Q

K

)
− (µa + ψ)Qw,

Θ3(t, Qw̄, Qw, Fw̄, Fw) =
ψ

2
Qw̄ + σFw − µw̄Fw̄,

Θ4(t, Qw̄, Qw, Fw̄, Fw) =
ψ

2
Qw − σFw − µwFw.

(36)
Thus, the proposed model (1) takes the form{

Dα
t ν(t) = R(t,ν(t)); t ∈ [0, b], 0 < α ⩽ 1,

ν(0) = ν0.
(37)

In view of this the problem (37) is given by,

ν(t) = ν0 + Iα
0R(t,ν(t))

= ν0 +
1

Γ(α)

∫ t

0
(t− k)α−1R(k,ν(k))dk. (38)

Let ε > 0 and consider the inequality given below

|Dα
t ν(t)−R(t,ν(t))| ≤ ε, t ∈ [0, b], (39)

where ε = max(εj)
T , j = 1, 2, 3, 4.



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 16 of 23

Definition 3. The proposed model (37) which is equivalent to system (1) is Ulam-Hyers stable
if there exists CR > 0, such that for every ε > 0 and a solution ν(t) ∈ B satisfies Eq. (39),
there exists a unique solution ν(t) ∈ B of system (1), with

|ν(t)− ν(t)| ≤ XRε, t ∈ J, (40)

where XR = max(CRj)
T .

Remark 1. A function ν(t) ∈ B satisfies inequality (39) if and only if there exists a function
h ∈ B with properties below:

(i) |h(t)| ≤ ε, t ∈ J .

(ii) Dα
t ν(t) = R(t,ν(t)) + h(t), t ∈ J .

Lemma 2. Assume that ν(t) ∈ B satisfies inequality (39), then ν(t) satisfies the integral
inequality described by

|ν(t)− ν0 −
1

Γ(α)

∫ t

0
(t− k)α−1R(k,ν(k))dk| ≤ Ωε. (41)

Proof. Thanks to (ii) of Remark 1,

Dα
t ν(t) = R(t,ν(t)) + h(t),

and Eq. (38), gives

ν(t) =ν0 +
1

Γ(α)

∫ t

0
(t− k)α−1R(k,ν(k))dk

+
1

Γ(α)

∫ t

0
(t− k)α−1h(k)dk.

(42)

Using (i) of Remark 1, we get

|ν(t)− ν0 −
1

Γ(α)

∫ t

0
(t− k)α−1R(k,ν(k))dk| ⩽ 1

Γ(α)

∫ t

0
(t− k)α−1|h(k)|dk ⩽ Ωε, (43)

where, Ω = bα

Γ(α+1) , t ∈ [a, b].

Theorem 5. Suppose that R : J ×R4 → R is continuous for every ν(t) ∈ B then the system
(37) which is the equivalence of system (1) is Ulam-Hyers stable.

Proof. Suppose that ν(t) ∈ B holds (39) and ν(t) is a unique solution of system
(37). Therefore, ∀ε > 0, t ∈ J and Lemma 2, it is obtained.

|ν(t)− ν(t)| ≤ max
t∈J

|ν(t)− ν0 −
1

Γ(α)

∫ t

0
(t− k)α−1R(k,ν(k))dk|

≤ max
t∈J

|ν(t)− ν0 −
1

Γ(α)

∫ t

0
(t− k)α−1R(k,ν(k))dk|



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 17 of 23

+max
t∈J

1

Γ(α)

∫ t

0
(t− k)α−1|R(k,ν(k)−R(k,ν(k))|dk

≤ |ν(t)− ν0 −
1

Γ(α)

∫ t

0
(t− k)α−1R(k,ν(k)dk|

+

[
LR
Γ(α)

] ∫ t

0
(t− k)α−1|ν(k)− ν(k)|dk

≤ Ωε+ΩLR|ν(k)− ν(k)|.

So,
||ν(t)− ν(t)|| ≤ XRε.

From Definition 3, the problem (1) has the Ulam-Hyers stability on t ∈ J .

7. Seasonal Effect

The environmental conditions, including temperature, humidity, and rainfall, sig-
nificantly affect the mortality of adult WI mosquitoes [34], leading to a sinusoidally forced
death rate.

CDαQw̄ =

[
ϕw̄F

2
w̄ + ρ1ϕwF

2
w + ρ2ϕwFwFw̄

Fw̄ + Fw

](
1− Q

K

)
− (µa + ψ)Qw̄, (44)

CDαQw =

[
(1− ρ1)ϕwF

2
w + (1− ρ2)ϕwFwFw̄
Fw̄ + Fw

](
1− Q

K

)
− (µa + ψ)Qw, (45)

CDαFw̄ =
ψ

2
Qw̄ + σFw −

[
µw̄

(
1− η cos

(
2π(t+ ω)

365

))]
Fw̄, (46)

CDαFw =
ψ

2
Qw − σFw −

[
µw

(
1− η cos

(
2π(t+ ω)

365

))]
Fw, (47)

Here, η represents the level of seasonal influence on the adult mortality rate, µw and µw̄ denote
the average mortality rates of adult WI and WU mosquitoes, respectively; t denotes time, and
ω indicates phase shift, which adjusts the seasonal phase in the cosine function. Because the
mosquito population is explicitly modeled, it is not necessary to introduce external seasonal
forcing terms. Since the mosquito population size is highly sensitive to the death rate, this
parameter was selected for seasonal forcing. As a result, the adult mosquito population
fluctuates seasonally as needed. This leads to seasonal changes in the aquatic population
which are appropriate since the mating function depends on population size. Seasonal effects
are shown in Figures 3 and 4.



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 18 of 23

0 100 200 300 400 500 600 700 800 900 1000

Days

2

3

4

5

6

7

8
105 Aquatic WU mosquitoes for   =0

=0.1

=0.3

=0.5

=0.7

=1

(a)

0 100 200 300 400 500 600 700 800 900 1000

Days

2

4

6

8

10

12

14

16
105 Aquatic WU mosquitoes for   =0.6228

=0.1

=0.3

=0.5

=0.7

=1

(b)

0 100 200 300 400 500 600 700 800 900 1000

Days

5

6

7

8

9

10

11

12
105 Aquatic WI mosquitoes for   =0

=0.1

=0.3

=0.5

=0.7

=1

(c)

0 100 200 300 400 500 600 700 800 900 1000

Days

2

4

6

8

10

12

14
105 Aquatic WI mosquitoes for   =0.6228

=0.1

=0.3

=0.5

=0.7

=1

(d)

Figure 3: Comparison between WI and WU mosquitoes with and without seasonal effects at
different values of α.



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 19 of 23

0 100 200 300 400 500 600 700 800 900 1000

Days

8.8

9

9.2

9.4

9.6

9.8

10
105 Adult WU mosquitoes for   =0

=0.1

=0.3

=0.5

=0.7

=1

(a)

0 100 200 300 400 500 600 700 800 900 1000

Days

0

0.5

1

1.5

2

2.5

3

3.5
106 Adult WU mosquitoes for   =0.6228

=0.1

=0.3

=0.5

=0.7

=1

(b)

0 100 200 300 400 500 600 700 800 900 1000

Days

5

5.1

5.2

5.3

5.4

5.5

5.6

5.7

5.8

5.9

6
105 Adult WI mosquitoes for   =0

=0.1 =0.3 =0.5 =0.7 =1

(c)

0 100 200 300 400 500 600 700 800 900 1000

Days

1

2

3

4

5

6

7

8

9
105 Adult WI mosquitoes for   =0.6228

=0.1

=0.3

=0.5

=0.7

=1

(d)

Figure 4: Comparison between adult WI and WU mosquitoes with and without seasonal
effect at different values of α.

Figures 3 and 4 show the time based dynamics of aquatic-stage and adult-stage
mosquitoes, both WI and WU, influenced by seasonally varying parameters. In this context,
η denotes the impact of seasonal forcing applied on the adult mortality rate µw, whereas α
indicates the fractional-order derivative integrating memory effects inside the system. As η
varies, seasonal variations in mortality cause periodic changes in mosquito population levels,
simulating authentic environmental dynamics. Furthermore, when α decreases, the oscillation
amplitude is reduced, signifying a slower system reaction and enhanced memory retention.
This combined influence of seasonal forcing and fractional-order dynamics suggests that the
model (44)–(47) effectively captures the recurrent and long-term behaviour of mosquito pop-
ulations under realistic ecological conditions.



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 20 of 23

8. Conclusions

This study developed and analyzed fractional-order mathematical models to in-
vestigate the introduction of Wolbachia into the Aedes aegypti mosquito population, including
imperfect maternal transmission and infection loss. The positivity and boundedness of the
solutions were demonstrated, validating the biological feasibility of the system. The invasive
reproduction number R0w|w̄ was derived to determine the conditions for Wolbachia persis-
tence. If R0w|w̄ < 1, the Wolbachia infection is not likely to spread and the WI population
will eventually die out after introduction into the wild mosquito population. If R0w|w̄ > 1,
then the WI mosquito population persists after being introduced into the wild population.
The equilibrium points were derived, and global stability was analyzed using Volterra-type
Lyapunov functions. In addition, the proposed model was shown to be Ulam-Hyers stable,
indicating that the approximate solution remains valid under small perturbations. The nu-
merical simulations, conducted using the Adams-Bashforth-Moulton approach, further aid in
understanding the dynamics across different fractional orders. The results indicate that as the
fractional order α decreases, the WI mosquito population increases, while the WU mosquito
population decreases. Additionally, the seasonal term was included to study the seasonal
effects on mosquito populations and observed fluctuations in the curves for different orders of
α. The study suggests that α = 0.5 yields optimal outcomes for increasing the WI mosquito
population, and the conditions used for α = 0.5 can be applied in real-life experiments to
increase the population of WI mosquitoes. This study offers important insights that may
assist in the examination of other mosquito-borne diseases, including Yellow fever, malaria,
Zika virus, and West Nile virus.

References

[1] R. V. Gibbons. Dengue: an escalating problem. BMJ, 324(7353):1563–1566, June 2002.
[2] M. B. Khan et al. Dengue overview: An updated systematic review. Journal of Infection

and Public Health, 16(10):1625–1642, August 2023.
[3] L. Alphey et al. Sterile-insect methods for control of mosquito-borne diseases: An anal-

ysis. Vector-Borne and Zoonotic Diseases, 10(3):295–311, September 2009.
[4] J. Bouyer and T. Lefrançois. Boosting the sterile insect technique to control mosquitoes.

Trends in Parasitology, 30(6):271–273, April 2014.
[5] World Health Organization (WHO). Global strategy for dengue prevention and control

2012-2020. Geneva, 2012.
[6] E.-E. Ooi, K.-T. Goh, and D. J. Gubler. Dengue prevention and 35 years of vector control

in singapore. Emerging Infectious Diseases, 12(6):887–893, June 2006.
[7] A. A. Hoffmann et al. Successful establishment of wolbachia in aedes populations to

suppress dengue transmission. Nature, 476(7361):454–457, August 2011.
[8] H. L. C. Dutra et al. Wolbachia blocks currently circulating zika virus isolates in brazilian

aedes aegypti mosquitoes. Cell Host and Microbe, 19(6):771–774, May 2016.
[9] L. A. Moreira et al. A wolbachia symbiont in aedes aegypti limits infection with dengue,

chikungunya, and plasmodium. Cell, 139(7):1268–1278, December 2009.
[10] A. P. Turley, L. A. Moreira, S. L. O’Neill, and E. A. McGraw. Wolbachia infection re-

duces blood-feeding success in the dengue fever mosquito, aedes aegypti. PLoS Neglected



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 21 of 23

Tropical Diseases, 3(9):e516, September 2009.
[11] T. H. Ant, C. S. Herd, V. Geoghegan, A. A. Hoffmann, and S. P. Sinkins. The wolbachia

strain wau provides highly efficient virus transmission blocking in aedes aegypti. PLoS
Pathogens, 14(1):e1006815, January 2018.

[12] T. Walker et al. The wmel wolbachia strain blocks dengue and invades caged aedes
aegypti populations. Nature, 476(7361):450–453, August 2011.

[13] P. A. Ross, S. A. Ritchie, J. K. Axford, and A. A. Hoffmann. Loss of cytoplasmic in-
compatibility in wolbachia-infected aedes aegypti under field conditions. PLoS Neglected
Tropical Diseases, 13(4):e0007357, April 2019.

[14] M. Turelli. Cytoplasmic incompatibility in populations with overlapping generations.
Evolution, 64(1):232–241, August 2009.

[15] D. E. Campo-Duarte, O. Vasilieva, D. Cardona-Salgado, and M. Svinin. Optimal con-
trol approach for establishing wmelpop wolbachia infection among wild aedes aegypti
populations. Journal of Mathematical Biology, 76(7):1907–1950, February 2018.

[16] P. Saha, G. C. Sikdar, J. K. Ghosh, and U. Ghosh. Disease dynamics and optimal
control strategies of a two-serotypes dengue model with co-infection. Mathematics and
Computers in Simulation, 209:16–43, February 2023.

[17] M. Rafikov, M. E. M. Meza, D. P. F. Correa, and A. P. Wyse. Controlling aedes aegypti
populations by limited wolbachia-based strategies in a seasonal environment. Mathemat-
ical Methods in the Applied Sciences, 42(17):5736–5745, March 2019.

[18] L. Xue, C. A. Manore, P. Thongsripong, and J. M. Hyman. Two-sex mosquito model
for the persistence of wolbachia. Journal of Biological Dynamics, 11(sup1):216–237,
September 2016.

[19] S. Javeed, D. Baleanu, A. Waheed, M. S. Khan, and H. Affan. Analysis of homotopy
perturbation method for solving fractional order differential equations. Mathematics,
7(1):40, January 2019.

[20] J. Dianavinnarasi, R. Raja, J. Alzabut, J. Cao, M. Niezabitowski, and O. Bagdasar.
Application of caputo–fabrizio operator to suppress the aedes aegypti mosquitoes via
wolbachia: An lmi approach. Mathematics and Computers in Simulation, 201:462–485,
February 2021.

[21] A. I. Adekunle, M. T. Meehan, and E. S. McBryde. Mathematical analysis of a wolbachia
invasive model with imperfect maternal transmission and loss of wolbachia infection.
Infectious Disease Modelling, 4:265–285, January 2019.

[22] K. S. Miller and B. Ross. An introduction to the fractional calculus and fractional dif-
ferential equations. Wiley, 1993.

[23] I. Podlubny. Fractional differential equations. Academic Press, San Diego, CA, USA,
1999.

[24] A. A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo. Theory and applications of fractional
differential equations. Elsevier Science Limited, 2006.

[25] M. A. Taneco-Hernández and C. Vargas-De-León. Stability and lyapunov functions for
systems with atangana–baleanu caputo derivative: An hiv/aids epidemic model. Chaos
Solitons and Fractals, 132:109586, January 2020.

[26] H. Khan, N. Kamran, S. Maqsood, D. K. Almutairi, and J. Alzabut. A logistic growth
epidemiological seir model with computational and qualitative results. European Journal
of Pure and Applied Mathematics, 18(2):5944, May 2025.

[27] A. Boukhouima, K. Hattaf, E. M. Lotfi, M. Mahrouf, D. F. M. Torres, and N. Yousfi.
Lyapunov functions for fractional-order systems in biology: Methods and applications.



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 22 of 23

Chaos Solitons and Fractals, 140:110224, August 2020.
[28] A. S. Rashed, M. M. Mahdy, S. M. Mabrouk, and R. Saleh. Fractional order mathematical

model for predicting and controlling dengue fever spread based on awareness dynamics.
Computation, 13(5):122, May 2025.

[29] H. Wan and J. Xu. How does the enhanced mortality of wolbachia-infected immature
mosquitoes affect dengue transmission? International Journal of Biomathematics, April
2024.

[30] J. Dianavinnarasi, R. Raja, J. Alzabut, M. Niezabitowski, and O. Bagdasar. Control-
ling wolbachia transmission and invasion dynamics among aedes aegypti population via
impulsive control strategy. Symmetry, 13(3):434, March 2021.

[31] G. M. Vijayalakshmi and M. Ariyanatchi. Fractional order modelling of wolbachia-
carrying mosquito population dynamics for dengue control. Modeling Earth Systems and
Environment, 11(4), May 2025.

[32] D. Joseph, R. Ramachandran, J. Alzabut, S. A. Jose, and H. Khan. A fractional-order
density-dependent mathematical model to find the better strain of wolbachia. Symmetry,
15(4):845, April 2023.

[33] U. Ufuktepe. Discrete wolbachia diffusion in mosquito populations with allee effects.
European Journal of Pure and Applied Mathematics, 15(4):1613–1622, October 2022.

[34] H. M. Yang, M. L. G. Macoris, K. C. Galvani, M. T. M. Andrighetti, and D. M. V.
Wanderley. Assessing the effects of temperature on the population of aedes aegypti, the
vector of dengue. Epidemiology and Infection, 137(8):1188–1202, February 2009.

[35] B. Zheng, M. Tang, J. Yu, and J. Qiu. Wolbachia spreading dynamics in mosquitoes
with imperfect maternal transmission. Journal of Mathematical Biology, 76(1–2):235–
263, June 2017.

[36] J. Arrivillaga and R. Barrera. Food as a limiting factor for aedes aegypti in water-storage
containers. PubMed, 29(1):11–20, June 2004.

[37] A. Boukhouima, H. Zine, E. M. Lotfi, M. Mahrouf, D. F. M. Torres, and N. Yousfi.
Lyapunov functions and stability analysis of fractional-order systems. In Elsevier eBooks,
pages 125–136. 2022.

[38] S. T. Ogunlade, A. I. Adekunle, M. T. Meehan, D. P. Rojas, and E. S. McBryde. Model-
ing the potential of wau-wolbachia strain invasion in mosquitoes to control aedes-borne
arboviral infections. Scientific Reports, 10(1), October 2020.

[39] A. A. Hoffmann et al. Stability of the wmel wolbachia infection following invasion into
aedes aegypti populations. PLoS Neglected Tropical Diseases, 8(9), September 2014.

[40] C. J. McMeniman and S. L. O’Neill. A virulent wolbachia infection decreases the via-
bility of the dengue vector aedes aegypti during periods of embryonic quiescence. PLoS
Neglected Tropical Diseases, 4(7):e748, July 2010.

[41] L. P. Lounibos and R. L. Escher. Sex ratios of mosquitoes from long-term censuses of
florida tree holes. Journal of the American Mosquito Control Association, 24(1):11–15,
March 2008.

[42] C. J. McMeniman et al. Stable introduction of a life-shortening wolbachia infection into
the mosquito aedes aegypti. Science, 323(5910):141–144, January 2009.

[43] A. Boukhouima, K. Hattaf, and N. Yousfi. Dynamics of a fractional order hiv infection
model with specific functional response and cure rate. International Journal of Differ-
ential Equations, 2017:1–8, January 2017.

[44] H.-L. Li, L. Zhang, C. Hu, Y.-L. Jiang, and Z. Teng. Dynamical analysis of a fractional-
order predator-prey model incorporating a prey refuge. Journal of Applied Mathematics



M. Yar, S. Javeed, T. Abbas Khan / Eur. J. Pure Appl. Math, 18 (3) (2025), 6522 23 of 23

and Computing, 54(1–2):435–449, May 2016.
[45] C. Vargas-De-León. Volterra-type lyapunov functions for fractional-order epidemic sys-

tems. Communications in Nonlinear Science and Numerical Simulation, 24(1–3):75–85,
January 2015.

[46] J. Dushoff, W. Huang, and C. Castillo-Chavez. Backwards bifurcations and catastrophe
in simple models of fatal diseases. Journal of Mathematical Biology, 36(3):227–248,
February 1998.

[47] S. M. Ulam. A Collection of Mathematical Problems. Interscience Publishers, New York,
1960.

[48] D. H. Hyers. On the stability of the linear functional equation. Proceedings of the National
Academy of Sciences, 27(4):222–224, April 1941.