Electronic Journal of Differential Equations, Vol. 2025 (2025), No. 116, pp. 1–18.

ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, DOI: 10.58997/ejde.2025.116

PERIODIC SOLUTION AND STATIONARY DISTRIBUTION OF A

STOCHASTIC EPIDEMIC MODEL WITH TWO DIFFERENT EPIDEMICS

AND DIFFERENT EPIDEMIOLOGICAL FRAMEWORKS

SHIVAM KUMAR MISHRA, SYED ABBAS, JUAN JOSE NIETO

Abstract. In epidemiology, more than one infectious disease can pose a risk to the host pop-
ulation. This area of study has attracted researchers in recent times. In this article, we have

considered an epidemic model that incorporates two different transmission techniques, namely

SIR and SIRS. The considered deterministic model has been perturbed stochastically at trans-
mission rates. The analysis has been done for the resulting stochastic model. Firstly, we explore

the existence of the positive T -periodic solution for the stochastic system. We determine that

the non-autonomous periodic version of the system with white noise has a positive periodic solu-
tion using the Lyapunov function and Khasminskii theory. In addition, we analyze the positive

recurrence of the system. The results obtained in this article give the idea that reducing the
white noise in the stochastic model is critical for observing positive T -periodic solutions and

positive recurrence. Finally, we present some examples and perform numerical simulations to

validate the theoretical results established in this study.

1. Introduction

To understand how different diseases interact or spread in a community is tricky for the re-
searchers studying diseases. When multiple diseases spread at the same time, it gets even more
complicated. Therefore, it is evident to introduce some special types of modeling techniques to
acknowledge how two diseases spread together in a community.

Deterministic models, characterized by their use of fixed parameters and precise equations, have
been instrumental in laying the groundwork for our understanding of how diseases spread through
populations. These models, based on differential equations, have provided valuable insights into
disease dynamics, transmission rates, and the impact of interventions in controlling epidemics. For
instance, in [8], a discrete time model has been considered by Hamer et al. They investigated the
evolution of recurrent measles outbreaks. Kermack et al. [12] conducted a mathematical analysis
of the dynamics of the epidemics within a homogeneous population, establishing a threshold con-
dition for outbreaks and analyzing how an epidemic progresses over time. Their work established
the foundation for compartmental models such as the SIR model, which divided the population
into susceptible, infected and recovered categories. Further, they developed the threshold the-
ory and constructed the traditional SIS epidemic model [13]. Nieto investigated a deterministic
SIR framework enhanced with anti-infectious control strategies, such as vaccination and treat-
ment, and demonstrated how these interventions could significantly alter the epidemic trajectory
through numerical simulations and analytical insights [28]. For a long period, researchers have
relied on the deterministic model to investigate infectious diseases. Several studies have been done
to investigate the dynamical behavior; see [9, 17, 18, 19, 25].

While deterministic models offer a solid framework for studying disease spread, they often
simplify the inherent randomness and variability present in real-world epidemics, leading to an
increased interest in stochastic modeling approaches that better capture the complexities and un-
certainties of disease transmission dynamics. Stochastic mathematical models of epidemics include

2020 Mathematics Subject Classification. 45J05, 34K50, 92D25.

Key words and phrases. Double epidemic hypothesis; periodic solution; positive recurrence.
©2025. This work is licensed under a CC BY 4.0 license.
Submitted October 15, 2025. Published December 16, 2025.

1



2 S. K. MISHRA, S. ABBAS, J. J. NIETO EJDE-2025/116

uncertainty or unpredictability in the spread of infectious diseases. In the context of epidemiology,
the use of stochastic models developed as a response to the limitations of deterministic models in
capturing the inherent variability and random occurrences seen in actual disease outbreaks. The
inclusion of randomness into models, which yielded a more realistic portrayal of the dynamics of
disease spread, is what made the stochastic approach popular.

In the mid 19th century, researchers started finding challenges dealing with random events and
uncertainty in various domains, which is when stochastic modeling was first introduced. However,
in the early to mid 20th century, a significant increase in the widespread use and formalization
of stochastic modeling, especially in the field of sciences, can be seen. Several contributions have
been made by stochastic models. For instance, Liu et al. [20] studied an HBV infection model with
stochastic perturbation incorporating logistic growth. They established the necessary condition
for the disease’s extinction and the presence of the ergodic stationary distribution. Bao et al. [3]
conducted an analysis on stochastic SIRS model including interval parameters and established the
definition for the stochastic basic reproduction number. Kuang et al. [16] investigated a stochastic
SIRS epidemic model incorporating saturated incidence. They deduced that smaller white noise
is required for the persistence of the epidemic and large noise will prevent the epidemic from
spreading. Zinhi et al. [34] investigated a stochastic SIRS model and established conditions for the
existence of stationary distributions and disease extinction. Similarly, Hening et al. [10] analyzed
the long-term behavior of SIQRS models and provided rigorous stochastic stability results. For
vector-borne and zoonotic diseases, Dutta et al. [6] explored the dynamics of the Nipah virus,
focusing on equilibrium analysis, sensitivity, and uncertainty quantification. From a control-
theoretic perspective, Russell and Cunniffe [31] highlighted a counterintuitive outcome where
optimal control may hinder disease eradication in an environment with the stochastic domain.
Liu et al. [21] analyzed a stochastic prey-predator model and determined the sufficient conditions
for the extinction of the predator population. Therefore, in recent times, many studies have been
done on stochastic epidemic models and researchers have demonstrated how environmental noise
affects population model dynamics (see [11, 22, 7]).

Figure 1. Population flow for the disease without interaction effect

In epidemiology, a single disease outbreak is the main focus of majority of the epidemic models,
even though multiple diseases can generate an epidemic. There are instances when two major
health issues occur simultaneously. It is same as dealing with two difficult situations together,
which makes things even more difficult for everyone. Moreover, adding randomness to a two-
pathogen model helps to better capture real-world disease dynamics. Random effects, like en-
vironmental or demographic changes, can affect how diseases spread, persist, or die out. These
features cannot be fully explained by models with only one pathogen or by purely deterministic
models. Several analysis has been done previously on double epidemics. Chang et al. [5] proposed
a hypothesis for a twofold epidemic asymmetry and constructed an SIRS epidemic model with two



EJDE-2025/116 STOCHASTIC EPIDEMIC MODEL WITH TWO DIFFERENT EPIDEMICS 3

distinct non-linear saturation incidence rates. A double epidemic hypothesis-based delayed SIR
disease transmission model has been adopted by Meng [24] and Meng et al. [26]. In [33, 27, 30],
the authors have considered a stochastic epidemic model for SIS that includes double epidemic
hypothesis with a saturated incidence rate.

Motivated by the above-mentioned works, in our study, we aim to more accurately reflect
real population dynamics by integrating two different epidemiological frameworks: SIR and SIRS
models. To illustrate this, we consider two infectious diseases, for example, measles, which is well-
categorized by the SIR epidemic model due to the development of long-lasting immunity after
infection, and influenza, which is often modeled using the SIRS framework because individuals
can lose immunity over time and become susceptible again. These examples allow us to capture
different aspects of infection and immunity dynamics within a single modeling framework. In
Figure 1, we have used two distinct transmission techniques to illustrate the disease spread model.
The deterministic model resulting from the transformation of the model dynamics depicted in
Figure 1 is involving a system of first order differential equations that results from converting the
model dynamics shown in Figure 1 into a mathematical model

dS(t)
dt

= Λ− β1S(t)I1(t)− β2S(t)I2(t)− γS(t) + µR2(t),

dI1(t)
dt

= β1S(t)I1(t)− γ1I1(t)− δ1I1(t),

dI2(t)
dt

= β2S(t)I2(t)− γ2I2(t)− δ2I2(t),

dR1(t)

dt
= δ1I1(t)− γ3R1(t),

dR2(t)

dt
= δ2I2(t)− γ4R2(t)− µR2(t).

(1.1)

Here the number of people who are susceptible to the disease at time t is denoted by S(t), and
the number of people who are infected with the first and second diseases at time t are denoted by
I1(t) and I2(t), respectively. R1(t) andR2(t) denote the number of individuals who have recovered
from the infections with first and second disease at time t, respectively. It is important to observe
that all of the parameters of the system (1.1) are positive. The birth and immigration rate of the
population is denoted by Λ. The infection transmission coefficient rates by the first and second
diseases are denoted by β1 and β2, respectively, while γ, γ1, γ2, γ3, and γ4 represent the mortality
rates of susceptible individuals, individuals infected by the first disease, individuals infected by
the second disease, individuals recovered from the first disease, and individuals recovered from
the second disease, respectively. δ1 and δ2 are the recovery rates for first and second diseases
respectively. Lastly, µ stands for the immunity loss ratio for the second disease.

Note that all the parameters Λ, β1, β2, γ, γ1, γ2, γ3, γ4, δ1, δ2 and µ are constants.
Furthermore, in terms of biology, it makes sense to assume that

γ ≤ min{γ1, γ2, γ3, γ4}. (1.2)

Remark 1.1. In classical models, the competitive exclusion principle states that two infections
giving complete cross-immunity cannot co-exist. However, in our model, co-existence becomes
possible because we combine the SIR and SIRS structures. In the SIRS framework, recovered
individuals can lose immunity and become susceptible again, which helps both infections persist.
This is further supported by the stability analysis of the deterministic system given in Section 3.

It is assumed that environmental noise primarily perturbs infection transmission rates β1 and
β2 by β1 + σ1Ḃ1(t) and β2 + σ2Ḃ2(t), where Ḃi(t) for i = 1, 2 are standard Brownian motions

satisfying Bi(0) = 0 and are mutually independent. The term Ḃi(t) denotes the formal time
derivative of Bi(t), representing white noise, for i = 1, 2. Furthermore, the parameters σ2

ij > 0
represent the intensities of the respective white noise, i, j = 1, 2. Using these perturbations, we



4 S. K. MISHRA, S. ABBAS, J. J. NIETO EJDE-2025/116

derive the subsequent stochastic model that corresponds to model (1.1):

dS(t) =
[
Λ− β1S(t)I1(t)− β2S(t)I2(t)− γS(t) + µR2(t)

]
dt− σ1S(t)I1(t)dB1(t)

− σ2S(t)I2(t)dB2(t),

dI1(t) = [β1S(t)I1(t)− γ1I1(t)− δ1I1(t)]dt+ σ1S(t)I1(t)dB1(t),

dI2(t) = [β2S(t)I2(t)− γ2I2(t)− δ2I2(t)]dt+ σ2S(t)I2(t)dB2(t),

dR1(t) = [δ1I1(t)− γ3R1(t)]dt,

dR2(t) = [δ2I2(t)− γ4R2(t)− µR2(t)]dt.

(1.3)

Note that the stochastic perturbations are incorporated via multiplicative noise terms of the
form βj + σjḂj(t), where βj for j = 1, 2 remains a fixed, strictly positive parameter. This ensures
that the transmission rate remains biologically meaningful (non-negative) at all times.

To our knowledge, no study has been conducted regarding the existence of periodic solution
and stationary distribution of the stochastic system (1.3). Although, the existence and uniqueness
of solutions of the stochastic system (1.3) has been validated in [32]. Our goal in this article is to
study the existence of positive T -periodic solutions and the stationary distribution of the solution
to the system (1.3).

The remaining sections of the article are organized as follows. In Section 2, we present some
preliminaries required to establish our results. In Section 3, we perform a dynamical analysis of
the deterministic system, including the stability analysis of equilibria. Section 4 and Section 5 are
devoted to the main results of this article, where we establish sufficient conditions for the existence
of positive T -periodic solutions and positive recurrence for the stochastic system (1.3). Numerical
simulations are provided in Section 6 to validate the theoretical findings. Finally, conclusions are
drawn in Section 7.

2. Preliminaries

It is obvious to observe that the study on epidemic over any population starts when there is an
epidemic outbreak happens and hence

(S(0), I1(0), I2(0),R1(0),R2(0)) ∈ ∆, (2.1)

which guarantees that the population is positive. For the definition of ∆, we refer to the following
remark:

Remark 2.1 ([4]). We denote the total population as N(t). Thus N(t) is given as N(t) =
S(t) + I1(t) + I2(t) +R1(t) +R2(t).

We may infer from the system (1.3) that

dN(t)

dt
= [Λ− γS(t)− γ1I1(t)− γ2I2(t)− γ3R1(t)− γ4R2(t)].

Applying inequality (1.2), we obtain

dN(t)

dt
≤ [Λ− γN(t)].

Therefore,

N(t) ≤
(
N(0)− Λ

γ

)
e−γt +

Λ

γ
≤ max

(
N(0),

Λ

γ

)
:= M. (2.2)

where M is a positive constant.
Consider a set

∆ =
{(

(S(t), I1(t), I2(t),R1(t),R2(t)
)
∈ R+

5 : S(t) + I1(t) + I2(t) +R1(t) +R2(t) ≤
Λ

γ

}
. (2.3)

Thus, the bound N(t) ≤ Λ
γ is verified by the fact that (S(0), I1(0), I2(0),R1(0),R2(0)) ∈ ∆.

Hence, the total population N(t) is bounded.

Notation



EJDE-2025/116 STOCHASTIC EPIDEMIC MODEL WITH TWO DIFFERENT EPIDEMICS 5

(1) Throughout this article, ⟨X⟩T represents the time average of the process X(t) over the
interval [0, T ]. Formally

⟨X⟩T =
1

T

∫ T

0

X(t)dt.

(2) In this paper, “a.s.” or ”a.s.” refers to almost surely, indicating that the specified property
holds with probability one.

(3) For any two real numbers c and d, the symbol c∧ d represents their minimum, while c∨ d
represents their maximum.

Definition 2.2 ([14]). Let t1, t2, . . . , tn be any arbitrary finite sequence and X(t) be a stochastic
process. If X(t1 + h), X(t2 + h), . . . , X(tn + h), the joint distribution of random variables is
independent of h, where h = mT for (m = ±1,±2, . . . ), then X(t) is regarded periodic with
period T .

Lemma 2.3 ([14]). Consider the stochastic differential equation defined as follows:

dy(t) = f(t, y(t))dt+ σ(t, y(t))dB(t). (2.4)

The functions f(t, y(t)) and σ(t, y(t)) defined in (2.4) are periodic in t with periods T and y ∈ Rn.
Let us assume that (2.4) possesses a unique global solution. Suppose that U(t, y(t)) ∈ C2 is a
function that is periodic in t with period T as well as holds the properties mentioned below:

(i)

inf
|y|>R

U(t, y(t)) → ∞ as R → ∞, (2.5)

(ii) outside some compact set, we have

LU(t, y(t)) ≤ −1, (2.6)

where L is the operator

L =
∂

∂t
+

n∑
i=1

fi(t, y)
∂

∂yi
+

1

2

∑
i,j=1n

ci,j(t, y)
∂2

∂yiyj
.

Then, system (2.4) possesses a T -periodic solution.

We now present some useful information about Markov chains and stochastic differential equa-
tions. Let Q = {1, 2, . . . , N} represent the state space and ζ(t) be the continuous-time Markov
chain that switches between the N regimes. In [23], Mao and Yuan proposed the following de-
scription of the Markov process (y(t), ζ(t)) ∈ Rn

+ ×Q:

dy(t) = P (y(t), ζ(t))dt+ σ(y(t), ζ(t))dB(t),

y(0) = y0 ∈ Rn
+, ζ(0) = ū ∈ Q,

where P (y(t), ζ(t)) : Rn ×Q → Rn, σ(x(t), ζ(t)) : Rn ×Q → Rn×d.
Furthermore, we define a linear operator L for a function V (x, k) ∈ C2(Rn)×Q, by

LV (y, l) =

n∑
i=1

bi(y, l)
∂V (y, l)

∂yi
+

n∑
i,j=1

cij(y, l)

2

∂2V (y, l)

∂yiyj
+

n∑
i=1

qki(x)V (x, i), c(y, l)

= σ(y, l)σT (y, l).

Definition 2.4. Given a set D and τD is the first reach time of D for the process Xx
t . If

P (τD < ∞) = 1 for each x /∈ D, the Markov process Xx
t with X0 = x is recurrent on D.

Moreover, τD is defined as follows

τD = inf{t > 0, Xx
t ∈ D}.

Furthermore, if E(τD) < ∞ for every x /∈ D, the process Xx
t is called positive recurrent on D.



6 S. K. MISHRA, S. ABBAS, J. J. NIETO EJDE-2025/116

3. Dynamical analysis of the deterministic system

This section explores the existence and uniqueness of a global positive solution by employing the
Picard-Lindelöf theorem and related results of classical theory of ordinary differential equations
[29]. Additionally, we establish conditions that ensure both local and global asymptotic stability
of equilibrium points, using Lyapunov’s direct method.

3.1. Existence and uniqueness of the positive solution. This section focuses on analyzing
the positivity and boundedness of solutions for the deterministic system (1.1) under the given
initial condition

(S(0), I1(0), I2(0),R1(0),R2(0)) ∈ Ξ. (3.1)

Since system (1.1) has locally Lipschitz continuous coefficients, it guarantees the existence of
a unique maximal local solution (S(t), I1(t), I2(t),R1(t),R2(t)) for all t ∈ [0, τ), where τ denotes
the explosion time. Also, we ensure that the initial condition (3.1) holds.

To establish the global existence and positivity of the solution, we now present the following
theorem, which ensures that the solution (S(t), I1(t), I2(t),R1(t),R2(t)) remains well defined for
all t ≥ 0 without breaking down in a finite time.

Theorem 3.1. For the initial condition (3.1), the solution of the system (1.1) exists globally and
remains positive for all t ≥ 0.

Proof. Consider (S(t), I1(t), I2(t),R1(t),R2(t)) as a solution of the system (1.1) under the initial
condition (3.1). We first establish the positivity of S(t).

Suppose, for contradiction, that there exists a first time t1 such that S(t1) = 0. From the first
equation of system (1.1), we have

dS(t)
dt

∣∣
t=t1

= Λ > 0.

This indicates that S(t) becomes negative for some t ∈ (t1 − h, t1), where h > 0 is a sufficiently
small real number, which contradicts the assumption that S(t) > 0 for all t ∈ [0, t1). Hence, we
conclude that S(t) > 0 for all t ∈ [0, τ).

Next, we analyze the positivity of I1(t) and I2(t). From the second and third equations of the
system (1.1), it follows that

Ij(t) = Ij(0)exp
(∫ t

0

[
β1S(v)− (γ1 + δ1)

]
dv

)
, j = 1, 2.

Since the exponential function is always positive and Ij(0) > 0 from the initial condition (3.1),
it follows that Ij(t) > 0 for all t ∈ [0, τ) and j = 1, 2.

Lastly, we verify the positivity of R(t). From the fourth and fifth equations of system (1.1), we
obtain

dR1(t) ≥ −γ3R1(t) for all t ∈ [0, τ),

dR2(t) ≥ (γ4 + µ)R1(t) for all t ∈ [0, τ),

which can further be deduced as

R1(t) ≥ R1(0)e
−γ3t for all t ∈ [0, τ),

R2(t) ≥ R2(0)e
(−γ4+µ)t for all t ∈ [0, τ).

Since R1(0) > 0 and R2(0) > 0 from (3.1), we conclude that Rj(t) > 0 for all t ∈ [0, τ) and
j = 1, 2.

As the positivity of the solution (S(t), I1(t), I2(t),R1(t),R2(t)) has been established, it remains
to show that the solution is globally defined. The boundedness of the solution, as established
by inequality (2.2), implies that it cannot blow up in finite time, which confirms that τ = ∞.
Therefore, the solution (S(t), I1(t), I2(t),R1(t),R2(t)) remains positive and exists for all t ≥ 0.
This completes the proof. □



EJDE-2025/116 STOCHASTIC EPIDEMIC MODEL WITH TWO DIFFERENT EPIDEMICS 7

3.2. Equilibrium analysis and local stability of equilibrium. To begin our equilibrium
analysis, we first determine the equilibrium points of the system (1.1). By solving the system for
steady-state conditions, we obtain the following equilibrium points:

E∗
0 =

(Λ
γ
, 0, 0, 0, 0

)
,

E∗
1 =

(γ1 + δ1
β1

,
−γγ1 − γδ1 + β1Λ

β(γ1 + δ1)
, 0,

−δ1(γγ1 + γδ1 − β1Λ)

β1γ3(γ1 + δ1)
, 0
)

=
(γ1 + δ1

β1
,
γ(R1 − 1)

β
, 0,

δ1γ(R1 − 1)

β1γ3
, 0
)
,

E∗
2 =

(γ2 + δ2
β2

, 0,
−(γγ2 + γδ2 − β2Λ)(γ4 + µ)

β2(γ2γ4 + γ4δ2 + γ2µ)
, 0,

−δ2(γγ2 + γδ2 − β2Λ)

β2(γ2γ4 + γ4δ2 + γ2µ)

)
=

(γ2 + δ2
β2

, 0,
γ(γ2 + δ2)(γ4 + µ)(R2 − 1)

β2(γ2γ4 + γ4δ2 + γ2µ)
, 0,

γδ2(γ2 + δ2)(R2 − 1)

β2(γ2γ4 + γ4δ2 + γ2µ)

)
,

where R1 = β1Λ
γ(γ1+δ1)

and R2 = β2Λ
γ(γ2+δ2)

. Using the threshold parameters R1 and R2, we can

determine the local stability of the equilibria.

Remark 3.2. We now define the invasion reproduction numbers for the system (1.1) correspond-
ing to each disease as follows

Rinv
2 =

β2S
∗

γ2 + δ2
; S∗ =

γ1 + δ1
β1

,

Rinv
1 =

β1S
∗∗

γ1 + δ1
; S∗ =

γ2 + δ2
β2

.

Here, S∗ and S∗∗ represent the susceptible population at the equilibria E∗
1 and E∗

2 , respectively.
Biologically, Rinv

2 measures the average number of secondary infections caused by an individual
of second disease when introduced into a population where first disease is already endemic (and
similarly for Rinv

1 ).
If Rinv

2 < 1, second disease cannot invade, similarly Rinv
1 < 1 prevents invasion by first disease.

We summarize the stability results in the following theorem.

Theorem 3.3. Consider system (1.1), the local stability of its equilibrium points is characterized
as follows:

(1) The disease-free equilibrium E∗
0 =

(
Λ
γ , 0, 0, 0, 0

)
is locally asymptotically stable if R1 < 1

and R2 < 1.
(2) The equilibrium E∗

1 is locally asymptotically stable if R1 > 1 and Rinv
2 < 1.

(3) The equilibrium E∗
2 is locally asymptotically stable if R2 > 1 and Rinv

1 < 1.

Proof. We will present the proof of the theorem sequentially.

(1) To analyze the local stability of system (1.1) at the disease-free equilibrium E∗
0 =

(
Λ
γ , 0, 0, 0, 0

)
,

we compute the corresponding Jacobian matrix

J (E∗
0 ) =


−γ −β1

Λ
γ −β2

Λ
γ 0 µ

0 β1
Λ
γ − (γ1 + δ1) 0 0 0

0 0 β2
Λ
γ − (γ2 + δ2) 0 0

0 δ1 0 −γ3 0
0 0 δ2 0 −(γ4 + µ)

 .

The eigenvalues of J (E∗
0 ) are:

E1 = −γ, E2 = β1
Λ

γ
− (γ1 + δ1), E3 = −γ3, E4 = β2

Λ

γ
− (γ2 + δ2), E5 = −(γ4 + µ). (3.2)

.



8 S. K. MISHRA, S. ABBAS, J. J. NIETO EJDE-2025/116

From (3.2), we observe that E1 < 0, E3 < 0 and E5 < 0. Furthermore, since R1 < 1 and R2 < 1,
it follows that E2 < 0 and E4 < 0. Therefore, since all the eigenvalues possess negative real parts,

the equilibrium point E∗
0 =

(
Λ
γ , 0, 0, 0, 0

)
is asymptotically stable.

(2) The Jacobian matrix corresponding to the equilibrium point E∗
1 =

(
γ1+δ1
β1

, γ(R1−1)
β , 0, δ1γ(R1−1)

β1γ3
, 0
)

of the system (1.1) is expressed as

J (E∗
1 ) =


−γ(R1 − 1)− γ −(γ1 + δ1) 0 0 µ

γ(R1 − 1) 0 0 0 0

0 0 β2(γ1+δ1)
β1

− (γ2 + δ2) 0 0

0 δ1 0 −γ3 0
0 0 δ2 0 −(γ4 + µ)

 .

The eigenvalues of J(E∗
1 ) are:

E1 = −(γ4 + µ) < 0,

E2 =
β2(γ1 + δ1)

β1
− (γ2 + δ2) = (γ2 + δ2)(Rinv

2 − 1) < 0,

E3 = −γ3 < 0,

and the remaining two eigenvalues are the roots of the quadratic equation

E2 + (γ(R1 − 1) + γ) E + γ(γ1 + δ1)(R1 − 1) = 0.

By the Routh-Hurvitz criteria [15] for quadratic equations, both roots have negative real parts
if and only if the coefficients of the above quadratic equation satisfy γ(R1 − 1) > 0 and γ(γ1 +
δ1)(R1 − 1) > 0. Since γ, γ1, δ1 > 0, these conditions hold if and only if R1 > 1.

Thus, all eigenvalues have negative real parts if R1 > 1. Therefore, the equilibrium E∗
1 is

asymptotically stable.
(3) This case is similar to the previous one, so the analysis is omitted. This completes the

proof. □

3.3. Global stability of disease-free equilibrium E∗
0 . After establishing the local stability of

the disease-free equilibrium E∗
0 , we now aim to establish its global stability. We begin with the

following theorem.

Theorem 3.4. Under the conditions R1 < 1 and R2 < 1, the disease-free equilibrium E∗
0 of

system (1.1) is globally asymptotically stable.

Proof. We use Lyapunov’s direct method to prove the result. Consider the following positive
definite Lyapunov function,

V ∗ =
(
S − Λ

γ
+ I1 + I2

)2

+ a1I12 + a2I22 +R1
2 +R2

2,

where ai’s (i = 1, 2) are positive constants to be chosen later.
Differentiating V ∗ along the solutions of system (1.1), we obtain

dV ∗

dt
= 2

(
S − Λ

γ
+ I1 + I2

)(
Λ− γS + µR2 − (γ1 + δ1)I1 − (γ2 + δ2)I2

)
+ 2a1I1

(
β1SI1 − (γ1 + δ1)I1

)
+ 2a2I2

(
β1SI2 − (γ2 + δ2)I2

)
+ 2R1(δ1I1 − γ3R1) + 2R2(δ2I2 − γ4R2 − µR2)

= −2γ
(
S − Λ

γ

)2

− 2
(
(γ1 + δ1)(a1 + 1)

)
I12 − 2

(
(γ2 + δ2)(a2 + 1)

)
I22

− 2γ3R1
2 − 2(γ4 + µ)R2

2 + 2µ
(
S − Λ

γ

)
R2 − 2(γ1 + δ1)

(
S − Λ

γ

)
I1

− 2(γ2 + δ2)
(
S − Λ

γ

)
I2 − 2γ

(
S − Λ

γ

)
I1 + 2µI1R2 − 2(γ2 + δ2)I1I2



EJDE-2025/116 STOCHASTIC EPIDEMIC MODEL WITH TWO DIFFERENT EPIDEMICS 9

+ 2γ
(Λ
γ
− S

)
+ 2µI2R2 − 2(γ1 + δ1)I1I2 + 2a1β1SI12 + 2a2β2SI22

+ 2δ1I1R1 + 2δ2I2R2

≤ −2γ
(
S − Λ

γ

)2

− 2(γ1 + δ1)I12 − 2(γ2 + δ2)I22 − 2γ3R1
2 − 2(γ4 + µ)R2

2

+ 2µ
(
S − Λ

γ

)
R2 + 2(γ + γ1 + δ1)

(Λ
γ
− S

)
I1 + 2(γ + γ2 + δ2)

(Λ
γ
− S

)
I2

+ 2µI1R2 + 2a1

(β1Λ

γ
− (γ1 + δ1)

)
I12 + 2a2

(β2Λ

γ
− (γ2 + δ2)

)
I22

+ 2δ1I1R1 + 2(µ+ δ2)I2R2.

Using Young’s inequality for the following expressions

2
(Λ
γ
− S

)
I1 ≤ c1

(Λ
γ
− S

)2

+
1

c1
I2
1 ,

2
(Λ
γ
− S

)
I2 ≤ c1

(Λ
γ
− S

)2

+
1

c1
I2
2 ,

2
(
S − Λ

γ

)
R2 ≤ c1

(
S − Λ

γ

)2

+
1

c1
R2

2,

2I1R1 ≤ 1

c2
I12 + c2R1

2,

2I1R2 ≤ 1

c3
I12 + c3R2

2,

2I2R2 ≤ 1

c3
I22 + c3R2

2,

where

c1 =
γ

γ1 + γ2 + γ3 + γ4 + δ1 + δ2
, c2 =

γ1 + δ1
γ2 + γ3 + γ4

,

and c3 is a positive constant.
For R1 ≤ 1 and R2 ≤ 1, it is possible to choose constants a1 and a2 such that

a1

(β1Λ

γ
− (γ1 + δ1)

)
+

1

c1
(γ + γ1 + δ1) +

1

c2
δ1 +

1

c3
µ = 0,

a2

(β2Λ

γ
− (γ2 + δ2)

)
+

1

c1
(γ1 + γ2 + δ2) +

1

c3
(µ+ δ2) = 0.

Thus, the derivative of the Lyapunov function V ∗ satisfies the inequality

dV ∗

dt
≤ −γ

(
S − Λ

γ

)2

− (γ1 + δ1)I12 − (γ2 + δ2)I22 − γ3R1
2 − (γ4 + µ)R2

2.

Therefore, V ∗ is negative definite, implying that the disease-free equilibrium E∗
0 of system (1.1) is

globally asymptotically stable, which concludes the proof of the theorem. □

4. Existence of positive T -periodic solutions

Periodicity is a very important property for epidemic models to plan and respond to public
health outbreak. It may help us to predict the severity and timing of the outbreak. In this section,
we establish the existence of periodic solution for our model (1.3) under certain conditions.

Theorem 4.1. Assume that with respect to the strong kernel, the stochastic system (1.3) possesses
a solution (S(t), I1(t), I2(t),R1(t),R2(t)). Then, the system (1.3) possesses a positive T -periodic
solution if

G =
⟨2Λ(t)β1⟩1/2T

⟨γ1 + γ2 + δ1 + δ2 + γ + (σ2
1 + σ2

2)⟩T
> 1.



10 S. K. MISHRA, S. ABBAS, J. J. NIETO EJDE-2025/116

Proof. Let us define a C2-function

V1 = −lnS − lnI1 − lnI2 + κ(t).

Using Ito’s formula, we have

LV1 =
−1

S
(
Λ− β1S(t)I1(t)− β2S(t)I2(t)− γS(t) + µR2(t)

)
+

1

2
σ2
1I2

1 (t) +
1

2
σ2
2I2

2 (t)

− 1

I1
(
β1S(t)I1(t)− γ1I1(t)− δ1I1(t)

)
+

1

2
σ2
1I2

1 (t)−
1

I2
(
β2S(t)I2(t)− γ2I2(t)− δ2I2(t)

)
+

1

2
σ2
2I2

2 (t) + κ̇(t)

= σ2
1I2

1 (t) + σ2
2I2

2 (t) + γ1 + γ2 + δ1 + δ2 + γ − Λ

S
+ β1I1 + β2I2 − µ

R2

S
− β1S − β2S + κ̇(t)

≤ (σ2
1 ∨ σ2

2) + γ1 + γ2 + δ1 + δ2 + γ − (2Λβ1)
1/2 − (2µβ2R2)

1/2 + β1I1 + β2I2 + κ̇(t).

Let us assume that

κ̇(t) = ⟨(σ2
1 ∨ σ2

2) + γ1 + γ2 + δ1 + δ2 + γ⟩T − 2β
1/2
1 ⟨Λ1/2⟩T

− (σ2
1 ∨ σ2

2) + γ1 + γ2 + δ1 + δ2 + γ + (2Λβ1)
1/2.

Using the above assumption, we obtain

LV1 ≤ ⟨(σ2
1 ∨ σ2

2) + γ1 + γ2 + δ1 + δ2 + γ⟩T − 2β
1/2
1 ⟨Λ1/2⟩T − (2µβ2R2)

1/2 + β1I1 + β2I2
≤ −⟨(σ2

1 ∨ σ2
2) + γ1 + γ2 + δ1 + δ2 + γ⟩T (H − 1) + β1I1 + β2I2 − (2µβ2R2)

1/2

≤ −A(H − 1) + β1I1 + β2I2 − (2µβ2R2)
1/2,

(4.1)

where A = inf{(σ2
1 ∨ σ2

2) + γ1 + γ2 + δ1 + δ2 + γ}. Let us now define

V2 = − lnS, V3 = − lnR1 − lnR2, V4 = − ln I1 − ln I2,

V5 =
(S + I1 + I2 +R1 +R2)

θ+1

θ + 1
, 0 < θ < 1.

Supposing N = S + I1 + I2 +R1 +R2, thus

V5 =
N θ+1

θ + 1
, 0 < θ < 1.

Using the similar computation as above, we obtain

LV2 ≤ γ +
1

2
(σ2

1 ∨ σ2
2)−

Λ

S
− µ

R2

S
.

By denoting B = sup{γ + 1
2 (σ

2
1 ∨ σ2

2)}, we obtain

LV2 ≤ B − Λ

S
− µ

R2

S
. (4.2)

Similarly, we have

LV3 ≤ C − δ1
I1
R1

− δ2
I2
R2

, (4.3)

where C = sup{γ3 + γ4 + µ}. Moreover, we have

LV4 ≤ D − (η1 − η2)S, (4.4)

where D = sup{γ1 + γ2 + δ1 + δ2 +
1
2 (σ

2
1 ∨ σ2

2)}. Applying Ito’s formula on V5, we obtain

LV5 = N θ(Λ− γS − γ1I1 − γ2I2 − γ3R1 − γ4R2) +
1

2
θN θ−1(2σ2

1S2I2
1 + 2σ2

2S2I2
2 ).

Suppose α = min{γ, γ1, γ2, γ3, γ4}. Then
LV5 ≤ N θ

(
Λ− αN

)
+ θN θ−1(σ2

1 ∨ σ2
2)

≤ −
(
α− θ

2
(σ2

1 ∨ σ2
2)
)(

Sθ+1 + Iθ+1
1 + Iθ+1

2 +Rθ+1
1 +Rθ+1

2

)
+ ΛN θ,



EJDE-2025/116 STOCHASTIC EPIDEMIC MODEL WITH TWO DIFFERENT EPIDEMICS 11

where α− θ
2 (σ

2
1 ∨ σ2

2) > 0 i.e. 0 < θ < 2α
σ2
1∨σ2

2
and

E = sup
{
− 1

2

(
α− θ

2
(σ2

1∨σ2
2)
)(

Sθ+1+Iθ+1
1 +Iθ+1

2 +Rθ+1
1 +Rθ+1

2

)}
+β(S+I1+I2+R1+R2)

θ.

After some computations, we obtain

LV5 ≤ −1

2

(
α− θ

2
(σ2

1 ∨ σ2
2)
)(

Sθ+1 + Iθ+1
1 + Iθ+1

2 +Rθ+1
1 +Rθ+1

2

)
+ E. (4.5)

We now consider a C2-function, which is defined as V : [0,∞)× R5
+ → R by

V (t,S, I1, I1, ,R1,R2) = MV1 + V2 + V3 + V4 + V5,

where the positive constant M satisfies −MA(H − 1) +B + C +D + E ≤ −2.

Note that V (t,S, I1, I1, ,R1,R2) → ∞ as |(S, I1, I1, ,R1,R2)| → ∞, ensuring that condition
(i) of Lemma 2.3 is satisfied. Now, we are required to prove condition (i) of Lemma 2.3.

Equations (4.1), (4.2), (4.3), (4.4) and (4.5) when combined yield

LV ≤ −MA(H − 1) +Mβ1I1 +Mβ2I2 −M(2µβ2R2)
1/2 +B − Λ

S
− µ

R2

S
+ C − δ1

I1
R1

− δ2
I2
R2

+D − (η1 − η2)S − 1

2

(
α− θ

2
(σ2

1 ∨ σ2
2)
)(

Sθ+1 + Iθ+1
1 + Iθ+1

2 +Rθ+1
1 +Rθ+1

2

)
+ E

≤ −2 +Mβ1I1 +Mβ2I2 −M(2µβ2R2)
1/2 − Λ

S
− µ

R2

S
− δ1

I1
R1

− δ2
I2
R2

− (η1 − η2)S

− 1

2

(
α− θ

2
(σ2

1 ∨ σ2
2)
)(

Sθ+1 + Iθ+1
1 + Iθ+1

2 +Rθ+1
1 +Rθ+1

2

)
. (4.6)

WE define a closed and bounded set

Θ =
{
(S, I1, I2,R1,R2) ∈ R5

+ : ϵ1 < S <
1

ϵ1
, ϵ2 < I1 <

1

ϵ2
, ϵ3 < I2 <

1

ϵ3
,

ϵ4 < R1 <
1

ϵ4
, ϵ5 < R2 <

1

ϵ5

}
,

where 0 < ϵ1, ϵ2, ϵ3, ϵ4, ϵ5 < 1 are sufficiently small real numbers. Clearly, Θ is a closed and
bounded subset of R5

+, hence compact. Furthermore, ϵ1, ϵ2, ϵ3, ϵ4 and ϵ5 satisfy the following
condition on the set R5

+\Θ:

0 < ϵ1 ≤ min
{[ 1

2

(
α− θ

2 (σ
2
1 ∨ σ2

2)
)
− (β1 + β2)

O + 1

] 1
θ+1

,
Λ

O + 1

}
,

0 < ϵ2 ≤ min
{[α− θ

2 (σ
2
1 ∨ σ2

2

4(O + 1)

]1/θ+1

,
1

Mβ1

}
,

0 < ϵ3 ≤ min
{[α− θ

2 (σ
2
1 ∨ σ2

2)

4(O + 1)

]1/θ+1

,
1

Mβ2

}
,

0 < ϵ4 ≤ min
{[α− θ

2 (σ
2
1 ∨ σ2

2)

2(O + 1)

]1/θ+1

,
δ1ϵ2
O + 1

}
,

0 < ϵ5 ≤ min
{[1

2

(
α− θ

2 (σ
2
1 ∨ σ2

2)
)
− µ

ϵ1

O + 1

]1/θ+1

,
δ2ϵ3
O + 1

}
,

where

O = −2 +Mβ1I1 +Mβ2I2 −M(2µβ2R2)
1/2 − 1

4

(
α− θ

2

(
σ2
1 ∨ σ2

2

))(
Iθ+1
1 + Iθ+1

2

)
.

To facilitate the further discussion, we divide Θc = R5
+\Θ into ten domains Θi; i = 1, 2, 3 . . . 10

such that (S, I1, I2,R1,R2) ∈ Θi for every i = 1, 2, 3 . . . 10 and

Θ1 =
{
S >

1

ϵ1

}
, Θ2 =

{
I1 >

1

ϵ2

}
,

Θ3 =
{
I2 >

1

ϵ3

}
, Θ4 =

{
R1 >

1

ϵ4

}
,



12 S. K. MISHRA, S. ABBAS, J. J. NIETO EJDE-2025/116

Θ5 =
{
0 < S < ϵ1

}
, Θ6 =

{
0 < S < ϵ1,R2 >

1

ϵ5

}
,

Θ7 =
{
0 < I1 < ϵ2

}
, Θ8 =

{
0 < I2 < ϵ3

}
,

Θ9 =
{
ϵ1 ≤ S ≤ 1

ϵ1
, ϵ2 ≤ I1 ≤ 1

ϵ2
, ϵ3 ≤ I2 ≤ 1

ϵ3
, 0 < R1 < ϵ4

}
,

Θ10 =
{
ϵ1 ≤ S ≤ 1

ϵ1
, ϵ2 ≤ I1 ≤ 1

ϵ2
, ϵ3 ≤ I2 ≤ 1

ϵ3
, ϵ4 ≤ R1 ≤ 1

ϵ4
, 0 < R2 < ϵ5

}
.

It can be observed here that Θc =
10⋃
i=1

Θi. Further, we prove that LV ≤ −1 for any (t,S, I1, I2,R1,R2) ∈

R+×Θc, i.e. we have to prove that LV ≤ −1 on all the ten domains mentioned above. We discuss
it in the following cases:

Case 1. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ1,

LV ≤ O − 1

2

(
α− θ

2
(σ2

1 ∨ σ2
2)
)
Sθ+1 − (β1 + β2)S

≤ O − 1

2

(
α− θ

2 (σ
2
1 ∨ σ2

2)
)

ϵθ+1
1

− (β1 + β2)

ϵ1
.

Since 0 < ϵ1 < 1, we have − 1
ϵ1

≤ − 1
ϵθ1
; θ > 0. Thus, we obtain

LV ≤ O − 1

2

(
α− θ

2 (σ
2
1 ∨ σ2

2)
)

ϵθ+1
1

− (β1 + β2)

ϵθ+1
1

≤ −1.

Case 2. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ2, we have

LV ≤ O − 1

4

(
α− θ

2
(σ2

1 ∨ σ2
2)
)
Iθ+1
1 ≤ O − 1

4

(
α− θ

2 (σ
2
1 ∨ σ2

2)
)

ϵθ+1
2

≤ −1.

Case 3. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ3, we obtain

LV ≤ O − 1

4

(
α− θ

2
(σ2

1 ∨ σ2
2)
)
Iθ+1
2 ≤ O − 1

4

(
α− θ

2 (σ
2
1 ∨ σ2

2)
)

ϵθ+1
3

≤ −1.

Case 4. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ4, we have

LV ≤ O − 1

4

(
α− θ

2
(σ2

1 ∨ σ2
2)
)
Rθ+1

1 ≤ O − 1

4

(
α− θ

2 (σ
2
1 ∨ σ2

2)
)

ϵθ+1
4

≤ −1.

Case 5. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ6, we obtain

LV ≤ O − Λ

S
< O − Λ

ϵ1
< −1.

Case 6. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ5, we obtain

LV ≤ O − 1

4

(
α− θ

2
(σ2

1 ∨ σ2
2)
)
Rθ+1

2 − µ
R2

S

≤ O − 1

4

(
α− θ

2 (σ
2
1 ∨ σ2

2)
)

ϵθ+1
5

− µ

ϵ1ϵ5
.

Since 0 < ϵ5 < 1, we have − 1
ϵ5

≤ − 1
ϵθ5
; θ > 0. Thus, we obtain

LV ≤ O − 1

4

(
α− θ

2 (σ
2
1 ∨ σ2

2)
)

ϵθ+1
5

− µ

ϵ1ϵ
θ+1
5

≤ −1.



EJDE-2025/116 STOCHASTIC EPIDEMIC MODEL WITH TWO DIFFERENT EPIDEMICS 13

Case 7. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ7, we have

LV ≤ −2 +Mβ1I1 ≤ −2 +Mβ1ϵ2 < −1.

Case 8. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ8, we obtain

LV ≤ −2 +Mβ2I2 ≤ −2 +Mβ2ϵ3 < −1.

Case 9. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ9, we obtain

LV ≤ O − δ1
I1
R1

≤ O − δ1
ϵ2
ϵ4

< −1.

Case 10. For any (t,S, I1, I2,R1,R2) ∈ R+ ×Θ10, we have

LV ≤ O − δ2
I2
R2

≤ O − δ2
ϵ3
ϵ5

< −1.

Thus, from the cases mentioned above, we can conclude that

LV ≤ −1 ∀ (t,S, I1, I2,R1,R2) ∈ R+ ×Θc,

where Θ ∈ R5
+ is a closed set, i.e., LV ≤ −1 outside a compact set R+ × Θ. Therefore, the

condition (2.6) of Lemma 2.3 holds. Hence, system (1.3) possesses a positive T -periodic solution
according to Lemma 2.3. □

5. Positive recurrence of the system

Positive recurrence is an essential property to study the long term behavior and persistence of
the epidemics. A theorem which is required to show that the system (1.3) is positive recurrent is
stated below:

Theorem 5.1. Consider (S, I1, I2,R1,R2) as a solution of system (1.3). Suppose that

H =
2
√
β1Λ

(σ2
1 ∨ σ2

2) + γ1 + γ2 + δ1 + δ2 + γ
> 1.

Then, the solution (S, I1, I2,R1,R2) is positive recurrent with respect to the domain

Dν =
{
(S, I1, I2,R1,R2) ∈ R5

+ : ν1 ≤ S ≤ 1

ν1
, ν2 ≤ I1 ≤ 1

ν2
, ν3 ≤ I2 ≤ 1

ν3
,

ν4 ≤ R1 ≤ 1

ν4
, ν5 ≤ R2 ≤ 1

ν5

}
,

where νi, {i = 1, 2, 3, 4, 5} are sufficiently small real numbers.

Proof. Let (X1, X2, X3, X4, X5) = (S, I1, I2,R1,R2). Define

V1 = − lnS − ln I1 − ln I2.
By applying Itô’s formula to V1, we obtain

LV1 = − 1

S
(
Λ− β1SI1 − β2SI2 − γS + µR2

)
− 1

I1
(
β1SI1 − γ1I1 − δ1I1

)
− 1

I2
(
β2SI2 − γ2I2 − δ2I2

)
+

1

2

(
σ2
1 + σ2

2

)
.

Simplifying gives

LV1 ≤ (σ2
1 ∨ σ2

2) + γ1 + γ2 + δ1 + δ2 + γ − (2Λβ1)
1/2 − (2µβ2R2)

1/2 + β1I1 + β2I2.
Hence,

LV1 ≤ −A′(H − 1) + β1I1 + β2I2 − (2µβ2R2)
1/2, (5.1)

where
A′ = (σ2

1 ∨ σ2
2) + γ1 + γ2 + δ1 + δ2 + γ.



14 S. K. MISHRA, S. ABBAS, J. J. NIETO EJDE-2025/116

Next, we define

V2 = − lnS, V3 = − lnR1 − lnR2, V4 = − ln I1 − ln I2,

V5 =
1

θ + 1
(S + I1 + I2 +R1 +R2)

θ+1,

where 0 < θ < 1. For V2,

LV2 ≤ B′ − Λ

S
− µ

R2

S
, (5.2)

where B′ = γ + 1
2 (σ

2
1 ∨ σ2

2). Similarly,

LV3 ≤ C ′ − δ1
I1
R1

− δ2
I2
R2

, (5.3)

where C ′ = γ3 + γ4 + µ, and
LV4 ≤ D′ − (η1 − η2)S, (5.4)

where D′ = γ1 + γ2 + δ1 + δ2 +
1
2 (σ

2
1 ∨ σ2

2).

Applying Itô’s formula to V5, we have

LV5 = N θ
(
Λ− αN

)
+

θ

2
N θ−1(σ2

1 ∨ σ2
2),

where N = S + I1 + I2 +R1 +R2 and α = min{γ, γ1, γ2, γ3, γ4}. Thus,

LV5 ≤ −1

2

(
α− θ

2
(σ2

1 ∨ σ2
2)
)∑

i

Xθ+1
i + E′, (5.5)

for some constant E′.
Now, we consider the C2 function

V = MV1 + V2 + V3 + V4 + V5,

where M > 0 is chosen such that

−MA′(H − 1) +B′ + C ′ +D′ + E′ ≤ −2.

Combining (5.1)-(5.5), we obtain

LV ≤ −2 +Mβ1I1 +Mβ2I2 −M(2µβ2R2)
1/2 − Λ

S
− µ

R2

S
− δ1

I1
R1

− δ2
I2
R2

− (η1 − η2)S

− 1

2

(
α− θ

2
(σ2

1 ∨ σ2
2)
)∑

i

Xθ+1
i .

Following the same argument as in Theorem 4.1, we conclude that

LV ≤ −1, for (S, I1, I2,R1,R2) ∈ R5
+ \Dν .

Hence, the expected return time satisfies

E(τDν
) ≤ V (S(0), I1(0), I2(0),R1(0),R2(0)) < ∞.

Therefore, the solution (S, I1, I2,R1,R2) of (1.3) is positive recurrent in the domain R5
+ \ Dν .

Hence, by Definition 2.4, the assertion is proved. □

Remark 5.2. We discuss a few special cases by restricting σ1 and σ2. The following are the
situations:

Case 1: When σ1 ̸= 0 and σ2 = 0. This represents that only the transmission rate for the first
disease has been perturbed stochastically. In this case, although the transmission rate for the
second disease has not been perturbed, changes in the dynamics of the first disease will indirectly
influence the second disease. For instance, the behavior or interactions of individuals afflicted
with the first disease (e.g., isolation, enhanced precautions) may change and impact the dynamics
of the second disease’s transmission. Over time, the stochastic perturbation in the first disease
transmission can accumulate and lead to population level changes. These changes may have
an indirect effect on susceptible, infected or recovered population from both the diseases. This
impacts the overall epidemic dynamics and changes the course of the second disease.



EJDE-2025/116 STOCHASTIC EPIDEMIC MODEL WITH TWO DIFFERENT EPIDEMICS 15

Case 2: When σ1 = 0 and σ2 ̸= 0. This represents that only the transmission rate for the second
disease has been perturbed stochastically. This case is similar to case 1. Thus, perturbing only the
second disease transmission can still have an indirect and complex effects on the overall epidemic
dynamics.

Case 3: When σ1 = 0 and σ2 = 0. This represents the deterministic case, which is system (1.1).

Remark 5.3. We have constructed similar Lyapunov functional for proving both the existence of
T -periodic solution and stationary distribution. The function has been chosen to satisfy the specific
requirements of both the results. Although, different Lyapunov functional can be constructed for
each problem. This choice was made in order to maintain uniformity and clarity, which provides
the reader a unified framework for study.

6. Numerical simulations

This section includes discussion for numerical simulation results of the proposed stochastic
models corresponding to the theoretical results proved in the previous sections.

Example 6.1. We have chosen the following parameter values in the system (1.3): Λ(t) = 0.1,
β1 = 0.1, β2 = 0.1, γ = 0.02, γ1 = 0.02, γ2 = 0.02, γ3 = 0.1, γ4 = 0.1, δ1 = 0.02, δ2 = 0.02,
σ1 = 0.1 and σ2 = 0.1. Further, we compute the condition

G =
(2Λβ1)

1
2(

γ1 + γ2 + δ1 + δ2 + γ + (σ2
1 + σ2

2)
) = 1.1785 > 1.

0 1000 2000 3000 4000 5000 6000 7000 8000

Time

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

P
o
p
u
la

ti
o
n

S

I1

I2

R1

R2

Figure 2. Trajectories of the stochastic system (1.3) obtained by the parameter
values given in Example 6.1 and showing the existence of a T -periodic solution.
The colored curves represent stochastic sample paths. Their oscillatory behavior
around a stable mean over long time intervals suggests the presence of a nontrivial
T -periodic solution in the presence of randomness.

To illustrate the existence of a positive T -periodic solution, we simulate the stochastic system
over a long time horizon. As shown in Figure 2, all compartments exhibit bounded fluctuations
and indicate recurrent oscillatory behavior, which suggests a form of periodicity in the presence of
stochastic perturbations. In particular, the infected compartment I1 exhibits significant periodic
oscillations, while the susceptible and recovered classes demonstrate bounded variations that align



16 S. K. MISHRA, S. ABBAS, J. J. NIETO EJDE-2025/116

with this periodicity. These results visually support our theoretical result based on Lyapunov
function techniques and Khasminskii’s theory.

The purpose of this example is to confirm the existence of a positive T -periodic solution in
accordance with Theorem 4.1, rather than to assess long-term coexistence or extinction of either
infection.

Example 6.2. In this example, we choose the following parameter values for system (1.3): Λ =
0.07, β1 = 0.1, σ2

1 = 0.02, σ2
2 = 0.04, so σ2

1 ∨ σ2
2 = 0.04, γ1 = 0.02, γ2 = 0.02, δ1 = 0.02, δ2 = 0.02,

γ = 0.02. We compute

H =
2
√
β1Λ

(σ2
1 ∨ σ2

2) + γ1 + γ2 + δ1 + δ2 + γ
= 1.19523 > 1.

0 20 40 60 80 100 120 140 160 180 200

Time

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

P
o
p
u
la

ti
o
n

S

I1

I2

R1

R2

Figure 3. Sample trajectories of the stochastic system (1.3) under parameter
values given in Example 6.2 and satisfying the positive recurrence condition given
in Theorem 5.1. The compartments S, I1, I2, R1, and R2 are observed to
fluctuate within bounded regions over a long period which indicates that the
stochastic system admits a stationary distribution.

Figure 3 illustrates the behavior of the stochastic system when the parameter values (given
in Example 6.2) fulfill the threshold condition for positive recurrence, as discussed in 5.1. The
trajectories of all compartments remain within a stable range over time, which indicates the
existence of a unique stationary distribution. It is important to note that, the infected class I2
stabilizes around a relatively higher mean compared to other classes, while I1 and the recovered
classes fluctuate mildly around their mean values.

7. Discussions and concluding remarks

This article introduces a stochastic epidemic model that incorporates two different epidemi-
ological frameworks, SIR and SIRS each governed by different nonlinear incidence rates. The
randomness in disease transmission, an important aspect in real-world epidemics, is modeled by
introducing stochastic perturbations in the infection rates. Specifically, in the system (1.1), the
deterministic transmission rates βjdt are replaced by βjdt + σjdBj(t) for j = 1, 2 to account for
environmental fluctuations.

Our primary objective was to explore the long-term behavior of this model (1.3) under stochastic
effects. We established sufficient conditions for the existence of a nontrivial positive T -periodic



EJDE-2025/116 STOCHASTIC EPIDEMIC MODEL WITH TWO DIFFERENT EPIDEMICS 17

solution, using Lyapunov functions and Khasminskii’s method for stochastic periodic systems. We
further showed the existence of a unique stationary distribution under an appropriate threshold
condition, signifying long-term disease persistence under randomness.

Compared to previous studies, this work deals with the additional complexity arising from hav-
ing two pathogens with different reinfection dynamics (SIRS and SIR), and multiplicative noise
terms affecting transmission. Proving the existence of stationary distribution and periodic solu-
tions under such a setting is significantly more involved due to the nonlinearity and the stochastic
interactions between the compartments.

This model offers valuable insights into the dynamics of two interacting infectious diseases;
however, it does have its limitations. For instance, we have assumed that stochasticity only enters
through the transmission rates, while recovery and mortality processes are kept deterministic.
The analysis is also restricted to specific noise structures (multiplicative Brownian motions), and
future work could explore more general Lévy noise or correlated perturbations.

A potential extension of this work can involve the development of digital twins, which serve as
real-time, data-integrated virtual representations of epidemic processes. Given that the present
model accounts for two distinct epidemic dynamics, a digital twin approach could enhance its
applicability by allowing real-time updates, scenario testing, and adaptive intervention strate-
gies. The concept and mathematical foundation for such digital twins using Stieltjes differential
equations, as proposed by Area et al. [1, 2], offer a novel perspective for extending the current
theoretical basis toward real-time and adaptive disease modeling.

From an epidemiological perspective, our model captures the interaction between co-circulating
pathogens in a fluctuating environment. However, for direct applicability to specific diseases, fur-
ther refinement is required, including parameter estimation from data and possibly incorporating
age structure, spatial effects, or vaccination strategies.

In summary, this work contributes to the mathematical understanding of stochastic epidemic
models involving multiple transmission routes. It lays the foundation for more realistic multi-
pathogen models and highlights the analytical challenges and richness of such systems.

Acknowledgments. Shivam Kumar Mishra was supported by the University Grants Commission
(UGC) with NTA ref. no. 201610181884. The authors want to thank the anonymous reviewers
for their valuable comments that improved the quality of the manuscript.

Author contributions. All authors contributed equally. The first draft was written by Shivam
Kumar Mishra. All authors read and approved the final version.

References

[1] I. Area, F. J. Fernández, J. J. Nieto, F. A. F. Tojo; Application of Stieltjes parabolic partial differential
equations to the population dynamics of Vespa Velutina, Nonlinear Analysis Real World Applications 86(2025),
104388.

[2] I. Area, F. J. Fernández, J. J. Nieto, F. A. F. Tojo, Concept and solution of digital twin based on a Stieltjes

differential equation, Mathematical Methods in the Applied Sciences, 45(12) (2022), 7451-7465.
[3] K. Bao, L. Rong, Q. Zhang; Analysis of a stochastic SIRS model with interval parameters, Discrete and

Continuous Dynamical Systems-B, 24(9) (2019), 4827-4849.
[4] B. Boukanjime, M. El Fatini, A. Laaribi, R. Taki; Analysis of a deterministic and a stochastic epidemic model

with two distinct epidemics hypothesis, Physica A: Statistical Mechanics and its Applications, 534(122321)
(2019), 0378-4371.

[5] Z. Chang, X. Meng, X. Lu; Analysis of a novel stochastic SIRS epidemic model with two different saturated
incidence rates, Physica A: Statistical Mechanics and its Applications, 472 (2017), 103-116.

[6] P. Dutta, G. Samanta, J. J. Nieto; Nipah virus transmission dynamics: equilibrium states, sensitivity and
uncertainty analysis, Nonlinear Dynamics, 113(9) (2025), 10617-10657.

[7] A. Gray, D. Greenhalgh, L. Hu, X. Mao, J. Pan, A Stochastic Differential Equation SIS Epidemic Model,
SIAM Journal on Applied Mathematics 71(3)(2011), 876-902.

[8] W. H. Hamer, M. A., M. D. Cantab; The Milroy lectures on epidemic disease in England—The evidence of

variability and persistence of type, The Lancet, 167(4305) (2006), 569-574.
[9] K. Hattaf, N. Yousfi, A. Tridane, Mathematical virus dynamics model with general incidence rate and cure

rate, Nonlinear Analysis: Real World Applications, 13(4) (2012), 1866-1872.
[10] A. Hening, D. H. Nguyen, T. Ta, S. C. Ungureanu; Long-term behavior of stochastic SIQRS epidemic models,

Journal of Mathematical Biology, 90(4) (2025), 41.



18 S. K. MISHRA, S. ABBAS, J. J. NIETO EJDE-2025/116

[11] N. T. Hieu, N. H. Du, P. Auger, N. H. Dang; Dynamical behavior of a stochastic SIRS epidemic model,

Mathematical Modelling of Natural Phenomena, 10(2) (2015), 56-73.

[12] W. O. Kermack, A. G. McKendrick; Contributions to the mathematical theory of epidemics, Bulletin of
Mathematical Biology 53(1-2) (1991), 33-55.

[13] W. O. Kermack, A. G. Mckendrick; Contributions to the mathematical theory of epidemics. II. The problem

of endemicity, Proceedings of the Royal Society A, 138(834) (1932), 55-83.
[14] R. Khasminskii; Stochastic Stability of Differential Equations, 2nd ed., Springer-Verlag, Heidelberg, 2012.

[15] G. A. Korn, T. M. Korn; Mathematical handbook for scientists and engineers: definitions, theorems, and

formulas for reference and review, Courier Corporation, 2000.
[16] D. Kuang, Q. Yin, J. Li; The threshold of a stochastic SIRS epidemic model with general incidence rate under

regime-switching, Journal of the Franklin Institute, 360(17) (2023), 13624-13647.

[17] T. Li, D. Acosta-Soba, A. Columbu, G. Viglialoro; Dissipative gradient nonlinearities prevent δ-formations in
local and nonlocal attraction-repulsion chemotaxis models, Studies in Applied Mathematics, 154(2) (2025),

e70018.
[18] T. Li, S. Frassu, G. Viglialoro; Combining effects ensuring boundedness in an attraction-repulsion chemotaxis

model with production and consumption, Zeitschrift für angewandte Mathematik und Physik, 74(3) (2023),

109.
[19] G. Li, Z. Jin; Global stability of an SEI epidemic model with general contact rate, Chaos, Solitons & Fractals,

23(3) (2005), 997-1004.

[20] Q. Liu, D. Jiang, N. Shi, T. Hayat, A. Alsaedi; Dynamical behavior of a stochastic HBV infection model with
logistic hepatocyte growth, Acta Mathematica Scientia, 37(4) (2017), 927-940.

[21] Q. Liu, D. Jiang, T. Hayat, A. Alsaedi; Dynamical behavior of stochastic predator-prey models with distributed

delay and general functional response, Stochastic Analysis and Applications, 38(3) (2019), 403-426.
[22] Q. Lu; Stability of SIRS system with random perturbations, Physica A: Statistical Mechanics and its Appli-

cations, 388(18) (2009), 3677-3686.

[23] X. Mao, C. Yuan; Stochastic Differential Equations with Markovian Switching, Imperial College Press, London,
2006.

[24] X. Meng; Stability of a novel stochastic epidemic model with double epidemic hypothesis, Applied Mathematics
and Computation, 217(2) (2010), 506-515.

[25] X. Meng, Z. Wu, T. Zhang; The dynamics and therapeutic strategies of a SEIS epidemic model, International

Journal of Biomathematics, 6(5) (2013), 1793-5245.
[26] X. Meng, S. Zhao, T. Feng, T. Zhang; Dynamics of a novel nonlinear stochastic SIS epidemic model with

double epidemic hypothesis, Journal of Mathematical Analysis and Applications, 433(1) (2016), 227-242.

[27] A. Miao, X. Wang, T. Zhang, W. Wang, B. G. Sampath Aruna Pradeep; Dynamical analysis of a stochastic
SIS epidemic model with nonlinear incidence rate and double epidemic hypothesis, Advances in Difference

Equations, 226(1) (2017).

[28] J. J. Nieto, SIR Epidemic Model: Anti-Infectious, J. Innovation Sciences and Sustainable Technologies,
4(2025), 283-295.

[29] L. Perko, Differential equations and dynamical systems, Springer Science & Business Media, 7 , 2013.

[30] H. Qi, L. Liu, X. Meng; Dynamics of a non-autonomous stochastic SIS epidemic model with double epidemic
hypothesis, Complexity 14(1)(2017), 1076-2787.

[31] R. Russell, N. J. Cunniffe, Optimal control prevents itself from eradicating stochastic disease epidemics, PLOS
Computational Biology, 21(2) (2025), e1012781.

[32] T. Tamil Selvan, M. Kumar; Analysis of a stochastic epidemic model driven by bilinear incidence rate with

two different transmission mechanisms, The Journal of Analysis, 32(1) (2024), 509-527.
[33] X. Zhang, D. Jiang, B. Ahmad, T. Hayat; Dynamics of a stochastic SIS model with double epidemic diseases

driven by Lévy jumps, Physica A: Statistical Mechanics and its Applications, 471(2017), 767-777.
[34] A. Zinihi, M. R. Sidi Ammi, M. Ehrhardt; Dynamical behavior of a stochastic epidemiological model: stationary

distribution and extinction of a SIRS model with stochastic perturbations, SeMA Journal (,2025), 1-22.

Shivam Kumar Mishra

School of Mathematical and Statistical Sciences, Indian Institute of Technology Mandi, Mandi, H.P.,

175005, India
Email address: shivammishra1807@gmail.com

Syed Abbas
School of Mathematical and Statistical Sciences, Indian Institute of Technology Mandi, Mandi, H.P.,

175005, India

Email address: abbas@iitmandi.ac.in, sabbas.iitk@gmail.com

Juan Jose Nieto

CITMAga, Departamento de Estat́ıstica, Análise Matemática e Optimización, Universidade de Santiago
de Compostela, 15782, Santiago de Compostela, Spain

Email address: juanjose.nieto.roig@usc.es