Next Article in Journal
Solvability of Mixed Problems for a Fourth-Order Equation with Involution and Fractional Derivative
Next Article in Special Issue
First-Passage Times and Optimal Control of Integrated Jump-Diffusion Processes
Previous Article in Journal
Existence Results for a Differential Equation Involving the Right Caputo Fractional Derivative and Mixed Nonlinearities with Nonlocal Closed Boundary Conditions
Previous Article in Special Issue
On the Estimation of the Persistence Exponent for a Fractionally Integrated Brownian Motion by Numerical Simulations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Simulation of Nonlinear Stochastic Analysis for Measles Transmission: A Case Study of a Measles Epidemic in Pakistan

1
School of Mathematic and Statistics, Jishou University, Jishou 416000, China
2
School of Computer Science and Cyber Engineering, Guangzhou University, Guangzhou 510006, China
3
Department of Mathematics, Sun Yat-sen University, Guangzhou 510275, China
*
Authors to whom correspondence should be addressed.
Fractal Fract. 2023, 7(2), 130; https://doi.org/10.3390/fractalfract7020130
Submission received: 14 December 2022 / Revised: 11 January 2023 / Accepted: 25 January 2023 / Published: 30 January 2023
(This article belongs to the Special Issue Stochastic Modeling in Biological System)

Abstract

:
This paper presents a detailed investigation of a stochastic model that rules the spreading behavior of the measles virus while accounting for the white noises and the influence of immunizations. It is hypothesized that the perturbations of the model are nonlinear, and that a person may lose the resistance after vaccination, implying that vaccination might create temporary protection against the disease. Initially, the deterministic model is formulated, and then it has been expanded to a stochastic system, and it is well-founded that the stochastic model is both theoretically and practically viable by demonstrating that the model has a global solution, which is positive and stochastically confined. Next, we infer adequate criteria for the disease’s elimination and permanence. Furthermore, the presence of a stationary distribution is examined by developing an appropriate Lyapunov function, wherein we noticed that the disease will persist for R 0 s > 1 and that the illness will vanish from the community when R 0 s < 1 . We tested the model against the accessible data of measles in Pakistan during the first ten months of 2019, using the conventional curve fitting methods and the values of the parameters were calculated accordingly. The values obtained were employed in running the model, and the conceptual findings of the research were evaluated by simulations and conclusions were made. Simulations imply that, in order to fully understand the dynamic behavior of measles epidemic, time-delay must be included in such analyses, and that advancements in every vaccine campaign are inevitable for the control of the disease.

1. Introduction

Measles is still a major worldwide health issue, particularly in the less developed countries. Measles (also known as Rubella or morbilli) is an extremely contagious illness caused mostly by the Morbillivirus genus in the Paramyxovirus [1,2]. Although efficient vaccines against this severe illness are commonly accessible, still measles is a leading cause of death among children below five of years of age [3]. The disease infecting hundreds of millions of children each year and resulting in a high mortality and morbidity in the child population, owing primarily to complicated conditions that exist side by side with the disease like malnutrition, diarrhea, and pneumonia [4]. Sneezing and coughing, touching the nasal or aerosol fluids, or close physical touch with an infected person are all ways to spread measles. It can stay extremely contagious for a maximum of two hours in the atmosphere and on the surfaces. Clinical manifestations include soar throat, cough, nasal congestion, blurred vision, and small white patches in the mouth; these signs and symptoms often develop within 10–12 days post-infection. Subsequently, a rash appears, extending downward out of the nose. The period of peak infectivity (virus shed) starts four days before that and four days just after commencement of the rashes. The usual incubation time is 14 days; however, it can range from 7–18 days [5]. In reality, even vaccinated people may still be susceptible if the immunization fails or existing immunity from the vaccine wears off. Despite the fact that vaccination has cut worldwide measles fatalities by 73% between 2000–2018, measles continues to remain a widespread in many underdeveloped nations, particularly in Asia and Africa. Around 140,000 individuals died from the measles in 2018. During 2000–2018, worldwide measles immunization results in an 85 percent drop in measles transmission [6,7]. Despite the abundance of a safe and effective vaccination in 2017, around 110,000 deaths occurred from measles, primarily children under the age of six, as per the report of the World Health Organization (WHO) [8]. Vaccination is amongst the most successful health promotion strategy for reducing death rates and the spreading of epidemics; it has been demonstrated that vaccination saves over 3 million people only in Nigeria every year. The vaccination will work with the immune system of the body to establish the body’s natural defenses, reducing the likelihood of relapse [9]. The MMR vaccination can protect against measles, which is a vaccine-preventable illness. The MMR vaccines are highly effective at protecting both adults and children from the epidemic measles. Only one dose of the MMR vaccines is roughly 92% successful in suppressing the measles, whereas two doses are around 95% effective. The MMR vaccine also protects against rubella and mumps [10]. This disease is an endemic one in Nigeria, with epidemics occurring at regular intervals. Measles is present across Nigeria at all seasons; however, it is more widespread in the summer months.
Pakistan is one of the most measles-affected nations in Asia [11]. Every 8–10 years, the nation has a recurring measles epidemic. In fact, 2845 identified measles infections were reported in Pakistan during the year 2016. This figure increased to 6791 in 2017 and, in the year 2018, 33,007 cases were reported. These results represent about 44, 20, and 51 percent of the total number of cases recorded in the East Mediterranean, which includes 22 nations. In 2017, over 130 children lost their lives due to this infection, with the figure rising to nearly 300 in 2018 [12].
It is strongly advised to employ the methods of mathematical modeling to examine the transmission process and prevention of an epidemic disease [7,13,14,15], modeling with fractional differential equations also have several applications in all fields of science [16,17,18]. While depicting the natural history of an infectious disease, the tools of modeling can create a balance among the data and its real biology. Thus far, models responsible for describing the dynamics of measles both from population to outbreak levels have demonstrated a wide variety of disease patterns. External noise is usually the main significant feature of physical processes and bio-systems. It has been discovered that environmental variables have a significant impact on the dynamics of measles disease spread [19]. Because of the uncertainty of person-to-person interactions or inherent characteristics of the population, outbreak onset and propagation are fundamentally unpredictable. As a result, the condition of the disease is influenced by the environment’s heterogeneity and uncertainty.
Changes in the environment likewise have a significant impact on the parasites’ persistence and distribution. Because the stochasticity of parameters and states depicts the exact dynamic behavior of an infectious disease, it is regarded as an essential part in epidemiological studies. Even though the perturbations are varying, these should be autocorrelated in a positive way. Furthermore, the perturbations may be estimated theoretically using the linked problem’s probability density function [20,21,22,23]. There are two main techniques to epidemic modeling: the deterministic modeling and stochastic modeling. Stochastic differential equation (SDE) models are recommended over deterministic models for mathematical modeling of biological functions because they may provide a higher level of reality than their deterministic equivalents [23,24,25,26]. We may choose to use SDEs to generate a distribution of the predicted output(s); for example, the number of infected individuals over time t. Moreover, when tested numerous times, a stochastic model produces more useful outputs than a deterministic one. A deterministic model, on the other hand, will produce only one outcome irrespective of the number of experiments. To explain the viral evolution of COVID-19, numerous deterministic infectious disease models have been suggested; for example, see [27,28].
The rest of the manuscript is organized in the following manner: Section 2 deals with the formulation of the proposed stochastic model for the spreading of the epidemic measles. The uniqueness and existence problem for obtaining a positive global solution is presented in Section 3. In Section 4 and Section 5, we characterized sufficient criteria for the stationary distribution and extinction of the disease. We optimize the proposed model using data from Pakistan compiled in the first ten months of 2019 in Section 6. We quantitatively compared the obtained analytical results, and graphical illustrations were presented in Section 7. We concluded the work in Section 8 by presenting the future directions and a comprehensive summary.

2. Model Formulation

Olumuyiwa et al. [29] have recently developed a model of the transmission of rubella disease by using the differential equations. Keeping in view the work of Olumuyiwa et al., here we intend to extend the model to a stochastic model. Furthermore, by considering different stages during the measles epidemic, we stratified the total population into six disjoint classes, namely: susceptible, vaccinated, exposed, infectious, hospitalized and recovered individuals whose sizes in mathematical terms are, respectively, S ( t ) , V ( t ) , E ( t ) , I ( t ) , H ( t ) , and R ( t ) .
The entrance of new persons through this population is captured by the rate ϕ and will be kept in the susceptible compartment. People in the vulnerable group start receiving a vaccination at a rate τ and setback immune function at a vaccine wane rate ω . The contact rate of susceptible persons is α , and thus the term force of infection becomes α S I , with the transition again from exposures to infection stages occurring at a rate of β . People who are infected with the measles seek medical attention at a rate of ρ and recover from the infection after the successful treatment supplemented at a rate of γ . We consider two types of death rates: the natural μ (that is constant for all classes) and the disease-induced mortality rate δ . This study assumes the recovery from measles that is possible due to the treatment only, that is, the study considering no natural recovery. The movements between the compartments are depicted via a flowchart in Figure 1. The above assumptions will lead to the following deterministic system:
d S ( t ) t = ϕ α S ( t ) I ( t ) + ω V ( t ) ( τ + μ ) S ( t ) , d V ( t ) t = τ S ( t ) ( μ + ω ) V ( t ) , d E ( t ) t = α S ( t ) I ( t ) ( μ + β ) E ( t ) , d I ( t ) t = β E ( t ) ( μ + δ + ρ ) I ( t ) , d H ( t ) t = ρ I ( t ) ( δ + γ + μ ) H ( t ) , d R ( t ) t = γ H ( t ) μ R ( t ) .
The threshold parameter has the following expression for model (1) as
R 0 = ( μ + ω ) ϕ β α ( μ + β ) ( μ + δ + ρ ) ( μ + ω + τ ) μ .
In order to consider the noises due the environment in the study (i.e., model (1)), we shall take into account the standard Brownian motions W i ( t ) for i = 1 , , 6 with W i ( 0 ) = 0 . Furthermore, system (1) is modified by considering the incidence rate of bi-linear form α S ( t ) I ( t ) N ( t ) . For the noise intensities, we have taken ξ 1 , ξ 2 , ξ 3 , ξ 4 , ξ 5 , ξ 6 as the relative weights. By considering these stochastic terms, the deterministic model (1) takes the form
d S = ϕ α S ( t ) I ( t ) N ( t ) + ω V ( t ) ( τ + μ ) S ( t ) d t + ξ 1 S ( t ) d W 1 ( t ) , d V = τ S ( t ) ( μ + ω ) V ( t ) d t + ξ 2 V d W 2 ( t ) , d E = α S ( t ) I ( t ) N ( t ) ( μ + β ) E ( t ) d t + ξ 3 E ( t ) d W 3 ( t ) , d I = β E ( t ) ( μ + δ + ρ ) I ( t ) d t + ξ 4 I ( t ) d W 4 ( t ) , d H = ρ I ( t ) ( μ + δ + γ ) H ( t ) d t + ξ 5 I ( t ) d W 5 ( t ) , d R = γ H ( t ) μ R ( t ) d t + ξ 6 R ( t ) d W 6 ( t ) .
Keeping in view system (3), the authors have a keen interest to find the possible answers to the following questions:
Q1:
What role do the random noises play in the transmission measles?
Q2:
What role contaminated vaccination in the spreading of measles disease?
Q3:
Under what condition(s) will the disease tend to go extinct?
Q4:
Under what condition(s) will the epidemic measles persist in the population?

3. Stochastic Model Analysis

This section investigates the existence/uniqueness of solutions, global asymptotic behavior, derivation of conditions under which the disease tends to go extinct, and the presence of an ergodic stationary distribution for the proposed stochastic model.

Positive Global Solution of the Model

The very first crucial question in studying the dynamic behavior is whether there is a possibility of the existence of a global solution to the model. Furthermore, for a system modeling the population dynamics, the nature of the solution’s value is of major relevance. In addition, we demonstrate in this section that the solution of randomized system (3) is global in nature and positive. It is well established that, for every given initial amount, the coefficients of a stochastic equation must fulfill the normal growth constraint and the local Lipschitz criterion in order to have a unique global solution (i.e., without any explosion in a limited period).
Theorem 1.
There exists a unique solution ( S ( t ) , V ( t ) , E ( t ) , I ( t ) , H ( t ) , R ( t ) ) of system (3) for t 0 under the initial conditions from the state R + 6 . Moreover, the solution remains in the same space (i.e., R + 6 ) surely for t 0 .
Proof. 
Keeping in view the Lipschitzness of the coefficients used in the model and from the fact ( ( S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) ) R + 6 , we can say that the proposed system has a unique local solution in [ 0 , τ e ) and t 0 . The term τ e stands for the explosion time, and readers are referred to [30,31] for a detailed analysis. By showing that, actually τ e = , we reach the conclusion that such a solution is global in nature. To do so, it is necessary that we assume a large k 0 > 0 in such a way that [ 1 k 0 , k 0 ] contains all parts of the solution. Assume k 0 k and let us define
τ k = i n f { t [ 0 , τ e ) : 1 k min { S ( t ) , V ( t ) , E ( t ) , I ( t ) , H ( t ) , R ( t ) } or k max { S ( t ) , V ( t ) , E ( t ) , I ( t ) , H ( t ) , R ( t ) } .
Whenever ϕ represents the empty set, then we shall write inf ϕ = . By increasing the value of k, one can notice that it increases τ k . We apply the limit k and assume that the τ k τ and a.s. τ e τ . Thus, if we show τ = a.s., it ensures τ e = . Proving all these guarantees that ( S ( t ) , V ( t ) , E ( t ) , I ( t ) , H ( t ) , R ( t ) ) R + 6 for any time t 0 . Let us assume the contrary case that τ e ; then, there must exist positive real numbers T and ϵ ( 0 , 1 ) in such a way
ϵ < P { τ T } .
Thus, for an integer k 0 k 1 , we have
P { T τ k } ϵ , k 1 k .
To begin, first let us establish a Lyapunov function of the type
G = ( S 1 log S ) + ( E log E 1 ) + ( I log I 1 ) + ( H log H 1 ) + ( R log R 1 ) ,
By utilizing the formula due to I t o ^ , letting k 0 k and assuming a very large non-negative real T, Equation (6) can be written in the form of
d G ( S , V , E , I , H , R ) = L G ( S , V , E , I , H , R ) d t + ξ 1 ( S 1 ) d W 1 ( t ) + ξ 2 ( V 1 ) d W 2 ( t ) + ξ 3 ( E 1 ) d W 3 ( t ) + ξ 4 ( I 1 ) d W 4 ( t ) + ξ 5 ( H 1 ) d W 5 ( t ) + ξ 6 ( R 1 ) d W 6 ( t ) .
In Equation (23), the L V operator is from the space R + 6 to R + . □
The remaining parts of the proof are merely similar to Theorem 2.1 in [26,30]. Thus, it is very simple for the reader to follow the result and and hence, is not completely proved here.

4. Extinction

While modeling the dynamical aspects of any epidemic diseases, it is important to investigate the situations under which the disease will become exterminated or tend to become extinct in the community. Within this section, we demonstrate that, when the white noises are large enough, the solution of the associated stochastic model (3) surely vanishes. Let us define
B ( t ) = 1 t 0 t B ( s ) d s .
Lemma 1.
(Strong Law)[32,33]Let Z = { Z } 0 t be continuous and real valued along with a local martingale, which vanishes as t 0 , then
lim t Z , Z t = , a . s . , lim t Z t Z , Z t = 0 , a . s . lim t sup Z , Z t t < 0 , a . s . , lim t Z t t = 0 , a . s .
Lemma 2.
Assume that ( S , V , E , I , H , R ) corresponds to initial data S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) in the space R + 6 and is a solution of model (3). Then,
lim sup t ln S ( t ) t = 0 , lim sup t V ( t ) t = 0 , lim sup t ln E ( t ) t = 0 , lim sup t ln I ( t ) t = 0 , lim sup t ln H ( t ) t = 0 , lim sup t ln R ( t ) t = 0 , a . s .
Furthermore, if μ > ξ 1 2 ξ 2 2 ξ 3 2 ξ 4 2 2 , and d > ξ 5 2 2 , then
lim t 0 t S ( s ) d W 1 ( s ) t = 0 , lim t 0 t V ( u ) d W 2 ( u ) t = 0 , lim t 0 t E ( u ) d W 3 ( u ) t = 0 , lim t 0 t I ( s ) d W 4 ( s ) t = 0 , lim t 0 t H ( s ) d W 5 ( s ) t = 0 , lim t 0 t R ( s ) d W 6 ( s ) t = 0 , a . s .
Then, the solution of (3)
lim sup t S ( t ) = ( μ + ω ) ϕ ( μ + ω + τ ) μ , lim sup t V ( t ) = ϕ τ ( μ + ω + τ ) μ , lim sup t E ( t ) = 0 , lim sup t I ( t ) = 0 , lim sup t H ( t ) = 0 , lim sup t R ( t ) = 0 , a . s .
To prove Lemma 2, we follow almost the same techniques as performed in proving Lemmas 2.1 and 2.2 carried out in the work of Zhao [32], and thus the readers can simply prove the results.
Next, to develop the extinction theory related to model (3), we first define the threshold quantity R s for the proposed stochastic model which can be written in the form of
R s = α ( μ + β ) ( μ + δ + ρ ) + ξ 3 2 2 + ξ 4 2 2 .
Theorem 2.
The exposed population E ( t ) and infected population I ( t ) of system ((3) tends to zero exponentially almost surely if R s < 1 , where R s is given by Equation (12).
Proof. 
Let ( S ( t ) , V ( t ) , E ( t ) , I ( t ) , H ( t ) , R ( t ) ) correspond to initial data S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) in the space R + 6 , being a solution of model (3).
Define
G 1 ( t ) = β E ( t ) + ( μ + β ) I ( t ) .
Differentiating Equation (13) following Ito’s formula, one can obtain
d ( l n G 1 ( t ) ) = 1 G 1 α β S I N ( μ + β ) ( μ + δ + ρ ) I β 2 E 2 ξ 3 2 + ( μ + β ) 2 ξ 4 2 I 2 2 ( G 1 ) 2 + β ξ 3 [ β E ( t ) + ( μ + β ) I ] E d W 3 ( t ) + ( μ + β ) ξ 4 [ E + ( μ + β ) I ] I d W 4 ( t ) , 1 G 1 α β I ( μ + β ) ( μ + δ + ρ ) I β 2 E 2 ξ 3 2 2 ( G 1 ) 2 ( μ + β ) 2 ξ 4 2 I 2 2 ( G 1 ) 2 + β ξ 3 [ β E ( t ) + ( μ + β ) I ] E d W 3 ( t ) + ( μ + β ) ξ 4 [ E + ( μ + β ) I ] I d W 4 ( t ) , [ S N ] 1 ( μ + β ) α ( μ + β ) ( μ + δ + ρ ) ξ 3 2 2 ξ 4 2 2 + β ξ 3 [ β E ( t ) + ( μ + β ) I ] E d W 3 ( t ) + ( μ + β ) ξ 4 [ E + ( μ + β ) I ] I d W 4 ( t ) . [ I I + β E ( μ + β ) ]
By taking the integral of either sides of the above inequality over the interval [ 0 , t ] , we have
l n G 1 ( t ) 1 ( μ + β ) α ( μ + β ) ( μ + δ + ρ ) + ξ 3 2 2 + ξ 4 2 2 + 0 t β ξ 3 E d W 3 ( t ) [ β E ( t ) + ( μ + β ) I ] + 0 t ( μ + β ) ξ 4 I d W 4 ( t ) [ E + ( μ + β ) I ] , 1 ( μ + β ) α ( μ + β ) ( μ + δ + ρ ) + ξ 3 2 2 + ξ 4 2 2 + 0 t β ξ 3 E d W 3 ( t ) [ β E ( t ) + ( μ + β ) I ] + 0 t ( μ + β ) ξ 4 I d W 4 ( t ) [ E + ( μ + β ) I ] , ( μ + β ) ( μ + δ + ρ ) + ξ 3 2 2 + ξ 4 2 2 ( μ + β ) R s 1 + 0 t β ξ 3 E d W 3 ( t ) [ β E ( t ) + ( μ + β ) I ] + 0 t ( μ + β ) ξ 4 I d W 4 ( t ) [ E + ( μ + β ) I ] .
By assuming the lim sup as t and multiplying both sides of relation (15) by 1 t while considering Lemma 2, we obtain
lim sup t ( l n G 1 ( t ) ) ( μ + β ) ( μ + δ + ρ ) + ξ 3 2 2 + ξ 4 2 2 ( μ + β ) R s 1 .
If R s < 1 , then lim t G 1 = lim t [ β E ( t ) + ( μ + β ) I ( t ) ] = 0 , a . s if R s < 1 . Again, β > 0 , ( μ + β ) > 0 , and we assert that lim t [ β E ( t ) + ( μ + β ) I ( t ) ] = 0 lim t E = lim t I = 0 —hence the result. □

5. The Stationary Distribution of the Disease

We understand that there are no endemic states in stochastic models. As a result, the stability of the system cannot be used to predict how long an illness would endure. As a result, one must concentrate on the existence/uniqueness assumption of the stationary distribution. In certain aspects, this assists the illness with survival. For this purpose, we employ a well-known method due to Hasminskii [34].
Assume a homogeneously time-dependent Markov process X ( t ) in the space R + n that satisfies the below stochastic model
d X ( t ) = b ( X ) d t + r k σ r d B r ( t ) .
The diffusion matrix can be demonstrated as:
A ( X ) = [ a i j ( ϰ ) ] , a i j ( ϰ ) = r = 1 k σ r i ( ϰ ) σ j r ( ϰ ) .
Lemma 3.
[34]. The process X ( t ) has a one and only one stationary distribution m ( . ) whenever there is a bounded domain having a regular boundary in such a way that U ¯ R d U ¯ closure U ¯ R d , and having the below characteristics
  • The smallest eigenvalue of the matrix A ( t ) is bounded below from the origin in the open domain U and in its neighborhood.
  • For ϰ R d U , the average time period τ (at which a path starts from ϰ reaching the set U) is bounded, and for every compact subset K R n , we have S u p ϰ k E ϰ τ < . When f ( . ) is an integrable function with the measure p i , then we have
P lim T 1 T 0 T f ( X ϰ ( t ) ) d t = R d f ( ϰ ) π ( d x ) = 1 ,
for all ϰ R d .
Let us define the following parameter for our future use:
R 0 s = μ β α ξ + μ + σ 2 2 2 α + μ + σ 3 2 2 δ + μ + σ 4 2 2 .
Theorem 3.
For R 0 s > 1 , then a solution ( S ( t ) , V ( t ) , E ( t ) , I ( t ) , H ( t ) , R ( t ) ) of system (3) is ergodic and has one and only one stationary distribution π ( . ) .
Proof. 
For verifying condition (2) of Lemma (3), we must develop a positive C 2 function G 2 : R + 6 R + . To do so, we first formulate
G 2 = S + V + E + I + H + R a 1 ln S a 2 ln E a 3 ln I ,
where a 1 , a 2 and a 3 are all real and positive constants, and must be calculated later on. By assuming the formula due to Itô’s and keeping in view system (3), we obtain
L ( V + S + I + E + H + R ) = ϕ μ ( V + S + I + E + H + R ) δ ( I + H ) , L ( ln S ) = Π S + α I N ω V S + ( τ + μ ) + ξ 1 2 2 , L ( ln V ) = τ S V + ( μ + ω ) + μ + ξ 2 2 2 , L ( ln E ) = α S I E N + ( μ + β ) + ξ 3 2 2 , L ( ln I ) = β E I + ( μ + δ + ρ ) + ξ 4 2 2 , L ( ln H ) = ρ I H + ( γ + δ + μ ) + ξ 5 2 2 , L ( ln R ) = γ H R + μ + ξ 6 2 2 .
Hence,
L G 2 = ϕ μ ( V + S + I + E + H + R ) δ ( I + H ) a 1 Π S + a 1 α I N a 1 ω V S + a 1 ( τ + μ ) + a 1 ξ 1 2 2 a 2 α S I E N + a 2 ( μ + β ) + a 2 ξ 3 2 2 a 3 β E I + a 3 ( μ + δ + ρ ) + a 3 ξ 4 2 2 .
The preceding calculation indicates that
L G 2 4 μ N × a 1 ϕ S × a 2 α S I E N × a 3 β E I 1 4 + a 1 ( τ + μ + ξ 1 2 2 ) + a 2 ( μ + β + ξ 3 2 2 ) + a 3 ( μ + δ + ρ + ξ 4 2 2 ) + ϕ + a 1 α I N a 1 ω V S δ ( I + H ) .
Suppose
a 1 ( τ + μ + ξ 1 2 2 ) = a 2 ( μ + β + ξ 3 2 2 ) = a 3 ( μ + δ + ρ + ξ 4 2 2 ) = ϕ .
Namely,
a 1 = ϕ ( τ + μ + ξ 1 2 2 ) , a 2 = ϕ ( μ + β + ξ 3 2 2 ) , a 3 = ϕ ( μ + δ + ρ + ξ 4 2 2 ) .
As a result,
L G 2 4 ϕ 4 μ α β ( τ + μ + ξ 1 2 2 ) ( μ + β + ξ 3 2 2 ) ( μ + δ + ρ + ξ 4 2 2 ) 1 4 ϕ + a 1 α I N a 1 ω V S δ ( I + H ) ,
L G 2 4 Π ( R 0 s ) 1 / 4 1 + a 1 α I N a 1 ω V S δ ( I + H ) .
Furthermore, we obtain
G 2 = a 4 ( E + S + V + H + I + R a 1 ln S a 2 ln E a 3 ln I ) ln S ln V ln H ln R E + S + V + H + I + R , = ( a 4 + 1 ) ( E + S + V + H + I + R ) ( a 1 a 4 + 1 ) ln S a 2 a 4 ln E a 3 a 4 ln I ln V ln H ln R ,
where a 4 > 0 is an unknown real number and must be determined later. It is very helpful to state
lim inf ( S , V , E , I , H , R ) R + 6 U k G 3 ( S , V , E , I , H , R ) = + , as k ,
where U k = ( 1 k , k ) × ( 1 k , k ) × ( 1 k , k ) .
In the below steps, we show that actually the function G 3 ( S , V , E , I , H , R ) has the minimum value, G 3 ( S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) . The partial derivatives of the function G 3 with respect to the state variables are given by:
G 3 ( S , V , E , I , H , R ) S = 1 + a 4 1 + a 1 a 4 S , G 3 ( S , V , E , I , H , R ) V = 1 + a 4 1 V , G 3 ( S , V , E , I , H , R ) E = 1 + a 4 a 2 a 4 E , G 3 ( S , V , E , I , H , R ) I = 1 + a 4 a 3 a 4 I , G 3 ( S , V , E , I , H , R ) H = 1 + a 4 c 3 c 4 H , G 3 ( S , V , E , I , H , R ) R = 1 + a 4 1 R .
One can verify very easily that the function G 3 has only one stagnation point.
( S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) = 1 + a 1 a 4 1 + a 4 , 1 1 + a 4 , a 2 a 4 1 + a 4 , a 3 a 4 1 + a 4 , 1 1 + a 4 , 1 1 + a 4 .
Furthermore, at ( S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) for V 2 ( S , V , E , I , H , R ) , the Hessian matrix is as follows:
B = 1 + a 1 a 4 S 2 ( 0 ) 0 0 0 0 0 0 1 V 2 ( 0 ) 0 0 0 0 0 0 a 2 a 4 E 2 ( 0 ) 0 0 0 0 0 0 a 3 a 4 I 2 ( 0 ) 0 0 0 0 0 0 1 H 2 ( 0 ) 0 0 0 0 0 0 1 R 2 ( 0 ) .
The Hessian matrix is obviously positive definite. As a result, the minimum value of G 3 ( S , V , E , I , H , R ) is G 3 ( S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) . We may assert that G 3 ( S , V , E , I , H , R ) has one and only one minimum value G 3 ( S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) inside R + 6 based on Equation (20) and the continuity of G 3 ( S , V , E , I , H , R ) . Then, as follows, we delineate a non-negative C 2 function G : R + 6 R + :
G ( S , V , E , I , H , R ) = G 3 ( S , V , E , I , H , R ) G 3 ( S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) .
Using I t o s formula and the suggested system, we arrive at
L ( G ) a 4 4 Π ( R 0 s ) 1 / 4 1 + a 1 α I N a 1 ω V S δ ( I + H ) + ϕ μ ( S + V + E + I + H + R ) δ ( I + H ) Π S + α I N ω V S + ( τ + μ ) + ξ 1 2 2 τ S V + ( μ + ω ) + μ + ξ 2 2 2 ρ I H + ( γ + δ + μ ) + ξ 5 2 2 γ H R + μ + ξ 6 2 2 .
Based on the above result, we have the following assertion:
L V a 4 a 5 + ( a 1 a 4 + 1 ) α I N ( a 1 a 4 + 1 ) ω V S δ ( a 4 + 1 ) ( I + H ) } + ϕ μ N Π S + ( τ + μ ) + ξ 1 2 2 τ S V + ( μ + ω ) + μ + ξ 2 2 2 ρ I H + ( γ + δ + μ ) + ξ 5 2 2 γ H R + μ + ξ 6 2 2 ,
where
a 5 = 4 ϕ ( R 0 s ) 1 / 4 1 > 0 .
Overall, for the solution to the model, we have the following set
D = { ϵ 1 < S < 1 ϵ 2 , ϵ 1 < V < 1 ϵ 2 , ϵ 1 < E < 1 ϵ 2 , ϵ 1 < I < 1 ϵ 2 , ϵ 1 < H < 1 ϵ 2 , ϵ 1 < R < 1 ϵ 2 } ,
where ϵ i > 0 are extremely small positive real values to be calculated later on. For the sake of simplicity, the whole set was partitioned into the following subsets:
D 1 = ( S , V , E , I , H , R ) R + 6 , 0 < S ϵ 1 , D 2 = ( S , V , E , I , H , R ) R + 6 , 0 < V ϵ 1 , S > ϵ 2 , D 3 = ( S , V , E , I , H , R ) R + 6 , 0 < E ϵ 1 , V > ϵ 2 , D 4 = ( S , V , E , I , H , R ) R + 6 , 0 < I ϵ 1 , E > ϵ 2 , D 5 = ( S , V , E , I , H , R ) R + 6 , 0 < H ϵ 1 , I > ϵ 2 , D 6 = ( S , V , E , I , H , R ) R + 6 , 0 < R ϵ 1 , I > ϵ 2 , D 7 = ( S , V , E , I , H , R ) R + 6 , S 1 ϵ 2 , D 8 = ( S , V , E , I , H , R ) R + 6 , V 1 ϵ 2 , D 9 = ( S , V , E , I , H , R ) R + 6 , E 1 ϵ 2 , D 10 = ( S , V , E , I , H , R ) R + 6 , I 1 ϵ 2 , D 11 = ( S , V , E , I , H , R ) R + 6 , H 1 ϵ 2 , D 12 = ( S , V , E , I , H , R ) R + 6 , R 1 ϵ 2 .
Then, we show that L G ( S , V , E , I , H , R ) < 0 on R + 6 D , which is the same is displaying the solution on the eight sub-regions.
Case 1.
If ( S , V , E , I , R ) D 1 , then, by Equation (24), we obtain
L G a 4 a 5 + ( a 1 a 4 + 1 ) α I N ( a 1 a 4 + 1 ) ω V S δ ( a 4 + 1 ) ( I + H ) } + ϕ μ N Π S + ( τ + μ ) + ξ 1 2 2 τ S V + ( μ + ω ) + μ + ξ 2 2 2 ρ I H + ( γ + δ + μ ) + ξ 5 2 2 γ H R + μ + ξ 6 2 2 , + ( a 1 a 4 + 1 ) α I N } + ϕ Π ϵ 1 + ( τ + μ ) + ξ 1 2 2 + ( μ + ω ) + μ + ξ 2 2 2 + ( γ + δ + μ ) + ξ 5 2 2 + μ + ξ 6 2 2 ,
For every ( S , V , E , I , H , R ) D 1 . , picking ϵ 1 > 0 , returns L G < 0 .
Just as in the proof above, we conclude that L G < 0 for ( S , V , E , I , H , R ) D i , i = 2 , 3 . . . 12 . .
As a result, we arrived to the point that there must be positive constant W > 0 in such a way
L G ( S , V , E , I , H , R ) < W < 0 for all ( S , V , E , I , H , R ) R + 6 D .
Thus,
d G ( S , V , E , I , H , R ) < W d t + [ ( a 4 + 1 ) S ( a 1 a 4 + 1 ) ξ 1 ] d W 1 ( t ) + [ ( a 4 + 1 ) V ξ 2 ] d W 2 ( t ) + [ ( a 4 + 1 ) E a 2 a 4 ξ 3 ] d W 3 ( t ) + [ ( a 4 + 1 ) I a 3 a 4 ξ 4 ] d W 4 ( t ) + [ ( a 4 + 1 ) H ξ 5 ] d W 5 ( t ) + [ ( a 4 + 1 ) R ξ 5 ] d W 6 ( t ) .
Consider ( S ( 0 ) , V ( 0 ) , E ( 0 ) , I ( 0 ) , H ( 0 ) , R ( 0 ) ) = ( ϰ 1 , ϰ 2 , ϰ 3 , ϰ 4 , ϰ 5 ) = ϰ R + 6 D , and τ ϰ is the amount of time it takes for a path to start from ϰ to achieve set D,
τ n = i n f { t : | X ( t ) | = n } & τ ( n ) ( t ) = min { τ ϰ , t , τ n } .
The next relation could be obtained if one integrates the above inequality from 0 to τ ( n ) ( t ) , considering the expectation and using Dynkin’s formula:
E G ( S ( τ ( n ) ( t ) ) , V ( τ ( n ) ( t ) ) , E ( τ ( n ) ( t ) ) , I ( τ ( n ) ( t ) ) , H ( τ ( n ) ( t ) ) , R ( τ ( n ) ( t ) ) ) G ( ϰ ) = E 0 τ ( n ) ( t ) L G ( S ( u ) , V ( u ) , E ( u ) , I ( u ) , H ( u ) , R ( u ) ) d u , E 0 τ ( n ) ( t ) W d u = W E τ ( n ) ( t ) .
Hence, G ( ϰ ) is non-negative; then,
E τ ( n ) ( t ) V ( ϰ ) W .
We have P { τ e = } = 1 as a result of the proof of Theorem (3). This also shows that model (3) is regular and, consequently, by letting n , t , almost surely we have τ ( n ) ( t ) τ ϰ .
Moreover, by utilizing the Fatou’s lemma, we have
E τ ( n ) ( t ) G ( ϰ ) W < ,
s u p ϰ K E τ ϰ < , where K is the compact set, i.e., a subset of R + 6 . This proves the condition (2) of Lemma (3) in an alternative approach.
Additionally, the diffusion matrix for the system (3) is
B = ξ 1 2 S 2 0 0 0 0 0 0 ξ 2 2 V 2 0 0 0 0 0 0 ξ 3 2 E 2 0 0 0 0 0 0 ξ 4 2 I 2 0 0 0 0 0 0 ξ 5 2 H 2 0 0 0 0 0 0 ξ 6 2 R 2 .
Taking M = min ( S , V , E , I , H , R ) D ¯ R + 6 { ξ 1 2 S 2 , ξ 2 2 V 2 , ξ 3 2 E 2 , ξ 4 2 I 2 , ξ 3 2 E 2 , ξ 5 2 H 2 , , ξ 6 2 R 2 } , we obtain
i , j = 1 6 a i j ( S , V , E , I , H , R ) ξ i ξ j = ξ 1 2 S 2 σ 2 + ξ 2 2 V 2 σ 2 2 + ξ 3 2 E 2 σ 2 + ξ 4 2 ξ 4 2 I 2 + ξ 5 2 σ 5 2 + ξ 6 2 σ 6 2 R 2 M | σ | 2 , ( S , V , E , I , H , R ) D ¯ ,
where ξ = ( σ 1 , σ 2 , σ 3 , σ 4 , σ 5 , σ 6 ) R + 6 . Similarly, this proves condition (1) of Lemma (3). Keeping in view the previous statements, the ergodicity of model (3) is insured by Lemma (3) and hence proving that the model has stationary distribution.

6. Parameter Estimation

Among the most important processes to carry through out the testing process is the utilization of real data (if available) to acquire findings for certain unidentified biological factors used in the epidemiology system under study. Real-world measles cases, as shown in Table 1, are used to test the proposed rubella model and to choose the best fitted parameters for numerous unknown biological parameters that emerge in the system. Considering the 2018 statistics of WHO, the natural death rate of a Pakistani individual is 1 / 66.5 since the life expectancy of a Pakistani is 66.5 . In addition, the entire size of the country is 207 , 862 , 518 , whereas the recruitment rate is calculated to be Π = 207 , 862 , 518 × μ 260 , 479 . Additionally, Memon et al. [13] indicates that the rate of measles vaccination is generally 97 percent effective, implying that the vaccines’ effectiveness τ equals 0.97 . Some values of the parameters are estimated and others were fitted against the data, and these were presented in Table 2. In Figure 2, the results via simulations for measles resurgence cases obtained by adapting the proposed model (3) with real data from the first ten months of the year 2019 are shown. As shown in Figure 2, it gives a rather good match, lending veracity to the predictions generated from the proposed measles model (3). We employed the relation 1 10 k = 1 10 ϰ k real ϰ k approximate ϰ k real 1.4685 × 10 1 to measure the associated relative error for fitting the model against the data.

7. Numerical Simulations and Discussion

It is of great concern to simulate each mathematical model against the real data and to verify theoretical results, and this step is very important when dealing with modeling biological problems. The researchers seek to simulate problem (3) employing classic numerical algorithms that rapidly converged. Almost all of the theoretical conclusions are quantitatively validated by using the well-known RK-4 approach.
S i + 1 = S i + ϕ α S i I i + ω V i ( τ + μ ) S i t + ξ 1 S i t ς 1 , i + ξ 1 2 2 S i ( ς 1 , i 2 1 ) t , V i + 1 = V i + τ S i ( μ + ω ) V i t + ξ 2 V i t ς 2 , i + ξ 2 2 2 V i ( ς 2 , i 2 1 ) t , E i + 1 = E i + α S i I i ( μ + β ) E i t + ξ 3 E i t ς 3 , i + ξ 3 2 2 E i ( ς 3 , i 2 1 ) t , I i + 1 = I i + β E i ( μ + δ + ρ ) I i t + ξ 4 I i t ς 4 , i + ξ 4 2 2 I i ( ς 4 , i 2 1 ) t , H i + 1 = H i + ρ I i ( γ + δ + μ ) H i t + ξ 4 H i t ς 5 , i + ξ 5 2 2 H i ( ς 5 , i 2 1 ) t , R i + 1 = R i + γ i μ R i t + ξ 6 R i t ς 5 , i + ξ 6 2 2 R i ( ς 6 , i 2 1 ) t .
where ς i , j ( i = 1 , 2 , 3 , 4 , 5 , 6 ) stands for the standard Gaussian stochastic variables, having distribution N ( 0 , 1 ) , and Δ t is the constant time-step. The terms ξ i > 0 , ( i = 1 , 2 , 3 , 4 , 5 , 6 ) reflect the intensities of the white noises.
To quantitatively validate the analytical conclusions, we require the parameters’ value being used in model (3). In Example 1 and 2, we established two sets of parameters’ values for this reason, and the population levels of each compartment at t = 0 were displayed. For every scenario, we simulated the model over the interval [ 0 , 80 ] .
We formulated Theorem 2 based on the stochastic theory of stability, which indicates that, under the condition of R s < 1 , the infection would continue to perish from the community regardless of the levels of the variables at t = 0 . Furthermore, the theory demonstrates that the infection will be eradicated from the community with unit probability. Figure 3 shows that the random curves reach the infection-free state after a limited period of time, confirming the analytical results.
Example 1.
The values of the parameter are assumed as: ϕ = 0.12 , τ = 0.002 , β = 0.008 , ω = 0.09 , μ = 0.005 , δ = 0.005 , α = 0.02 , γ = 0.004 and ρ = 0.05 , where the initial values of the state vector are S ( 0 ) = 50 , V ( 0 ) = 30 , E ( 0 ) = 10 , I ( 0 ) = 15 , H ( 0 ) = 25 , R ( 0 ) = 20 . Similarly, the intensities of the white noises are: ξ 1 = 0.55 , ξ 2 = 0.25 , ξ 3 = 0.25 , ξ 4 = 0.33 , ξ 5 = 0.55 , ξ 6 = 0.50 . Using those same model parameters, we estimated R s , which was found to be lower than 1. As a result, the assumption of Theorem 2 is met, and hence the components of the solution of the model adhere to the following relations:
lim t sup log E ( t ) t 0 , a . s .
lim t sup log I ( t ) t 0 , a . s .
Essentially, these relations explain the elimination of the virus from the community, as seen by Figure 3. As a consequence, the derived research findings on elimination are valid and may be relied on.
Likewise, Theorem 3 guarantees the disease’s prevalence in the population given permissible limits. By considering data from Example 2, we calculated the value of R 0 s , and it was found that it is greater than unity. Figure 4 depicts the numerical results based on the theorem’s assumptions. The figure implies that the infection persists inside the population, and that there is persistence of the solution of the proposed model (3). This verifies the conclusion of Theorem 3 in the case of deterministic model (1). These results further elaborate that, when the related threshold parameter of the perturbed system exceeds unity, the solution of the model (3) lies within the neighborhood of endemic equilibrium. Thus, policymakers must provide strong preventative measures against various variations in order to restrict the spreading of multiple strains throughout the community. Moreover, Theorem 3 ensures the existence of an ergodic stationary distribution for model (3), and it is confirmed by Figure 5.
Example 2.
In this case, the chosen parameter values are of the form: ϕ = 2.12 , β = 0.08 , ω = 0.07 , μ = 0.001 , δ = 0.005 , α = 0.2 , ρ = 0.004 and γ = 0.08 . Similarly, the initial population sizes of the state variables are S ( 0 ) = 50 , V ( 0 ) = 30 , E ( 0 ) = 10 , I ( 0 ) = 15 , H ( 0 ) = 25 , R ( 0 ) = 20 , whereas the intensities are given by ξ 1 = 0.50 , ξ 2 = 0.35 , ξ 3 = 0.70 , ξ 4 = 0.50 , ξ 5 = 0.41 , ξ 5 = 0.45 . Considering the measles data and hence the estimated and fitted parameters, we find that R 0 s exceeds unity. It is also explored that the model parameters taken in this example satisfy the premise of Theorem 3. We tested the model using this input, and the outcomes are displayed visually in Figure 4. The figure shows that the disease tends to stay inside the community, and the model exhibits a homogeneous stationary distribution in this situation.
From Theorem 3, there is a single stationary distribution of system (3). To numerically illustrate this statistical property, we present in Figure 4 and Figure 5, the trajectories and the associated joint density function for each class of the population. For a good visibility, we offer the 3 D and upper view of the aforementioned joint densities in Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10. This indicates that the infection is still present in the population over time. We talk here about persistence in the mean of the epidemic. In Figure 4, we illustrate the continuation of all groups of the studied population. We remark that the stochastic trajectories fluctuate around the deterministic solution with reasonable distances according to magnitude of the noises.

8. Conclusions and Future Research

In this study, we presented a detail analysis of a stochastic model that describes the spreading behavior of the measles virus while accounting for the white noises and the influence of immunizations. It is assumed that the randomness being used in the model is nonlinear, and that a person may lose the resistance after vaccination, implying that vaccination might create temporary protection against the disease. First of all, we formulated a deterministic model and then it was expanded to a stochastic model. It is shown that the stochastic model is both theoretically and practically viable by demonstrating that the model has a global solution which is positive and stochastically bounded. Next, we developed sufficient criteria for the disease’s elimination and permanence. Furthermore, the presence of a stationary distribution is examined by constructing a suitable Lyapunov function, wherein we noticed that the disease will persist for R 0 s > 1 and that the illness will vanish from the community for R 0 s < 1 . We simulated the model against the available data of measles in Pakistan during the first ten months of 2019, by using the conventional curve fitting methods, and the values of the parameters were calculated therein. These values of the parameters were used in simulating the model, and the theoretical findings of the research were evaluated and conclusions were made. Simulations of the study suggest that, in order to fully understand the dynamic behavior of the measles epidemic, time-delay must be included in such analyses, and that advancements in every vaccine campaign are unavoidable in order to stop or minimize the spreading of the disease.
We fit both stochastic and determinism models to known data on rubella in Pakistan during the first ten months of 2019, and the values of parameters were obtained using lsqcurvefit methods. These model parameters are used, and the threshold number, which is around 13, is determined. This demonstrates that measles is extremely harmful and might have a negative impact on this community if adequate control tactics are not implemented in time. It also suggests that, if no appropriate measures are implemented, the cases of the measles may increase in the coming years. To further reduce the transmission rate, health authorities and lawmakers must launch awareness campaigns, mass immunizations, particularly among youngsters, and, most crucially, care and treatment for people having chronic conditions. This was also discovered that α serves as the most critical indicator to the threshold parameter; thus, lowering the disease spreading co-efficient (by lowering the household and sexual interactions of chronic patients with vulnerable) for acutely infected individuals who become chronic is also an effective control to reduce the spread of measles.
In the next research studies, the authors intend to use fractional calculus and modify this and the related model while utilizing different definitions of fractional derivatives such as Riesz, Caputo, Atangana–Baleanu, Caputo–Fabrizio, and many others for capturing the real dynamics of such and related diseases.

Author Contributions

Methodology, B.G. Formal analysis, B.G. Investigation, B.G. Data curation, A.K. and A.D.; Writing—original draft, B.G. and A.D. Writing—review and editing, A.K. Supervision, A.D. Project administration, A.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research was sponsored by the Guangzhou Government Project under Grant No. 62216235, and the National Natural Science Foundation of China (Grant No. 12250410247).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

Authors are grateful to the handling editor and reviewers for their constructive comments and remarks which helped to improve the quality of the manuscript.

Conflicts of Interest

The authors declare that they have no competing interests.

References

  1. Mose, O.F.; Sigey, J.K.; Okello, J.A.; Okwoyo, J.M.; Kang’ethe, G.J. Mathematical modeling on the control of measles by vaccination: Case study of KISII County, Kenya. SIJ Trans. Comput. Sci. Eng. Appl. CSEA 2014, 2, 61–69. [Google Scholar]
  2. Garba, S.M.; Safi, M.A.; Usaini, S. Mathematical model for assessing the impact of vaccination and treatment on measles transmission dynamics. Math. Meth. Appl. Sci. 2017, 40, 6371–6388. [Google Scholar] [CrossRef]
  3. Roberts, M.G.; Tobias, M.I. Predicting and preventing measles epidemic in New Zealand: Application of mathematical model. Epidem. Infect. 2000, 124, 279–287. [Google Scholar] [CrossRef] [PubMed]
  4. World Health Organization. Measles, Preprint. 2018. Available online: https://www.who.int/news-room/fact-sheets/detail/measles (accessed on 27 December 2019).
  5. Perry, R.T.; Halsey, N.A. The clinical significance of measles: A review. J. Infect. Dis. 2004, 189, S4–S16. [Google Scholar]
  6. Ejima, K.; Omori, R.; Aihara, K.; Nishiura, H. Real-time investigation of measles epidemics with estimate of vaccine efficacy. Int. J. Biol. Sci. 2012, 8, 620–629. [Google Scholar] [CrossRef] [Green Version]
  7. Mossong, J.; Muller, C.P. Modelling measles re-emergence as a result of waning of immunity in vaccinated populations. Vaccine 2003, 21, 4597–4603. [Google Scholar] [CrossRef]
  8. Bolarin, G. On the dynamical analysis of a new model for measles infection. Int. J. Math. Trends Technol. 2014, 7, 144–154. [Google Scholar] [CrossRef]
  9. Taiwo, L.; Idris, S.; Abubakar, A.; Nguku, P.; Nsubuga, P.; Gidado, S. Factors affecting access to information on routine immunization among mothers of under 5 children in Kaduna state Nigeria, 2015. Pan Afr. Med. J. 2017, 27, 186. [Google Scholar] [CrossRef]
  10. Center for Disease Control. Available online: https://www.cdc.gov/vaccines/vpd/measles/index.html (accessed on 26 January 2021).
  11. World Health Organization. Eastern Mediterranean Vaccine Action Plan 2016–2020: A Framework for Implementation of the Global Vaccine Action Plan (No. WHO-EM/EPI/353/E); World Health Organization, Regional Office for the Eastern Mediterranean: Cairo, Egypt, 2019. [Google Scholar]
  12. Dawn. Curbing Measles. Available online: https://www.dawn.com/news/1520931 (accessed on 7 December 2019).
  13. Memon, Z.; Qureshi, S.; Memon, B.R. Mathematical analysis for a new nonlinear measles epidemiological system using real incidence data from Pakistan. Eur. Phys. J. Plus 2020, 135, 1–21. [Google Scholar] [CrossRef]
  14. Asamoah, J.K.K.; Nyabadza, F.; Jin, Z.; Bonyah, E.; Khan, M.A.; Li, M.Y.; Hayat, T. Backward bifurcation and sensitivity analysis for bacterial meningitis transmission dynamics with a nonlinear recovery rate. Chaos Solitons Fractals 2020, 140, 110237. [Google Scholar] [CrossRef]
  15. Khan, F.M.; Khan, Z.U.; Lv, Y.-P.; Yusuf, A.; Din, A. Investigating of fractional order dengue epidemic model with ABC operator. Results Phys. 2021, 24, 104075. [Google Scholar] [CrossRef]
  16. Li, X.-P.; Bayatti, H.A.; Din, A.; Zeb, A. A vigorous study of fractional order COVID-19 model via ABC derivatives. Results Phys. 2021, 29, 104737. [Google Scholar] [CrossRef]
  17. Din, A.; Li, Y. Lévy noise impact on a stochastic hepatitis B epidemic model under real statistical data and its fractal–fractional Atangana–Baleanu order model. Phys. Scr. 2021, 96, 124008. [Google Scholar] [CrossRef]
  18. Srivastava, H.M.; Saad, K.M. Numerical simulation of the fractal-fractional Ebola virus. Fractal Fract. 2020, 4, 49. [Google Scholar] [CrossRef]
  19. Liu, P.; Ikram, R.; Khan, A. The measles epidemic model assessment under real statistics: An application of stochastic optimal control theory. Comput. Methods Biomech. Biomed. Eng. 2022, 26, 138–159. [Google Scholar] [CrossRef]
  20. Liu, G.; Qi, H.; Chang, Z.; Meng, X. Asymptotic stability of a stochastic May mutualism system. Comput. Math. Appl. 2020, 79, 735–745. [Google Scholar] [CrossRef]
  21. Tran, K.Q.; Yin, G. Optimal harvesting strategies for stochastic ecosystems. IET Control. Theory Appl. 2017, 11, 2521–2530. [Google Scholar] [CrossRef]
  22. Sabbar, Y.; Din, A. Probabilistic analysis of a marine ecological system with intense variability. Mathematics 2022, 10, 2262. [Google Scholar] [CrossRef]
  23. Din, A. The stochastic bifurcation analysis and stochastic delayed optimal control for epidemic model with general incidence function. Chaos Interdiscip. J. Nonlinear Sci. 2021, 31, 123101. [Google Scholar] [CrossRef]
  24. Khan, A.; Sabbar, Y. Stochastic modeling of the Monkeypox 2022 epidemic with cross-infection hypothesis in a highly disturbed environment. Math. Biosci. Eng. 2022, 19, 13560–13581. [Google Scholar] [CrossRef]
  25. Din, A.; Khan, A.; Sabbar, Y. Long-Term Bifurcation and Stochastic Optimal Control of a Triple-Delayed Ebola Virus Model with Vaccination and Quarantine Strategies. Fractal Fract. 2022, 6, 578. [Google Scholar] [CrossRef]
  26. Din, A.; Ain, Q.T. Stochastic Optimal Control Analysis of a Mathematical Model: Theory and Application to Non-Singular Kernels. Fractal Fract. 2022, 6, 279. [Google Scholar] [CrossRef]
  27. Omame, A.; Abbas, M.; Din, A. Global asymptotic stability, extinction and ergodic stationary distribution in a stochastic model for dual variants of SARS-CoV-2. Math. Comput. Simul. 2022, 204, 302–336. [Google Scholar] [CrossRef] [PubMed]
  28. Liu, P.; Rahman, M.u.; Din, A. Fractal fractional based transmission dynamics of COVID-19 epidemic model. Comput. Methods Biomech. Biomed. Eng. 2022, 25, 1852–1869. [Google Scholar] [CrossRef]
  29. Olumuyiwa, J.P.; Ojo, M.M.; Viriyapong, R.; Oguntolu, F.A. Mathematical model of measles transmission dynamics using real data from Nigeria. J. Differ. Equations Appl. 2022, 28, 753–770. [Google Scholar]
  30. Jin, X.; Jia, J. Qualitative study of a stochastic SIRS epidemic model with information intervention. Phys. A Stat. Mech. Appl. 2020, 547, 123866. [Google Scholar] [CrossRef]
  31. Rajasekar, S.P.; Pitchaimani, M. Qualitative analysis of stochastically perturbed SIRS epidemic model with two viruses. Chaos Solitons Fractals 2019, 118, 207–221. [Google Scholar] [CrossRef]
  32. Bao, K.; Zhang, Q. Stationary distribution and extinction of a stochastic SIRS epidemic model with information intervention. Adv. Differ. Equations 2017, 2017, 352. [Google Scholar] [CrossRef] [Green Version]
  33. Zhao, Y.N.; Jiang, D.Q. The threshold of a stochastic SIS epidemic model with vaccination. Appl. Math. Comput. 2014, 243, 718–727. [Google Scholar] [CrossRef]
  34. Khasminskii, R. Stochastic Stability of Differential Equations; Springer Science and Business Media: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
Figure 1. Flow chart of the measles model (1) [29].
Figure 1. Flow chart of the measles model (1) [29].
Fractalfract 07 00130 g001
Figure 2. The figure shows the fitting of both ODE and SDE models against the real data by using values of the parameters shown in Table 2.
Figure 2. The figure shows the fitting of both ODE and SDE models against the real data by using values of the parameters shown in Table 2.
Fractalfract 07 00130 g002
Figure 3. Sample solutions of stochastic systems (3) with its associated deterministic counterpart.
Figure 3. Sample solutions of stochastic systems (3) with its associated deterministic counterpart.
Fractalfract 07 00130 g003
Figure 4. Sample solution profiles of the system (3) correspond to its deterministic counterpart.
Figure 4. Sample solution profiles of the system (3) correspond to its deterministic counterpart.
Fractalfract 07 00130 g004
Figure 5. Ergodic stationary distribution of (3).
Figure 5. Ergodic stationary distribution of (3).
Fractalfract 07 00130 g005
Figure 6. The right panel of the figure shows the joint two-dimensional densities at t = 3000 of individuals S , V , E and I of system (3) correspond to data from Example 2—Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of S , V , E and I collectively (in this case, α = 1.99 ).
Figure 6. The right panel of the figure shows the joint two-dimensional densities at t = 3000 of individuals S , V , E and I of system (3) correspond to data from Example 2—Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of S , V , E and I collectively (in this case, α = 1.99 ).
Fractalfract 07 00130 g006
Figure 7. The right panel of the figure shows the joint two-dimensional densities at t = 3000 of individuals S , H , R , V and E of system (3) correspond to data from Example 2—Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of S , H , R , V and E collectively (in this case, α = 1.99 ).
Figure 7. The right panel of the figure shows the joint two-dimensional densities at t = 3000 of individuals S , H , R , V and E of system (3) correspond to data from Example 2—Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of S , H , R , V and E collectively (in this case, α = 1.99 ).
Fractalfract 07 00130 g007
Figure 8. The right panel of the figure shows the joint two-dimensional densities at t = 3000 of individuals V , I , H , and R of system (3) correspond to data from Example 2—Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of V , I , H , and R collectively (in this case, α = 1.99 ).
Figure 8. The right panel of the figure shows the joint two-dimensional densities at t = 3000 of individuals V , I , H , and R of system (3) correspond to data from Example 2—Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of V , I , H , and R collectively (in this case, α = 1.99 ).
Fractalfract 07 00130 g008
Figure 9. The right panel of the figure shows the joint two-dimensional densities at t = 3000 of individuals E , I , H , and R of system (3) correspond to data from Example 2—Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of E , I , H , and R , collectively (in this case, α = 1.99 ).
Figure 9. The right panel of the figure shows the joint two-dimensional densities at t = 3000 of individuals E , I , H , and R of system (3) correspond to data from Example 2—Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of E , I , H , and R , collectively (in this case, α = 1.99 ).
Fractalfract 07 00130 g009
Figure 10. (First part) The right panel of the figure shows that the joint two-dimensional densities at t = 3000 of individuals I , H , and R of system (3) correspond to data from Example 2–Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of I , H , and R collectively (in this case, α = 1.99 ).
Figure 10. (First part) The right panel of the figure shows that the joint two-dimensional densities at t = 3000 of individuals I , H , and R of system (3) correspond to data from Example 2–Test 2 (2nd column), where different colors depict the density sizes. The panel on the left demonstrates the 3D graph of the two-dimensional densities of I , H , and R collectively (in this case, α = 1.99 ).
Fractalfract 07 00130 g010
Table 1. Real cases of the measles epidemic in Pakistan during January and October 2019.
Table 1. Real cases of the measles epidemic in Pakistan during January and October 2019.
Jan (01)Feb (02)Mar (03)Apr (04)May (05)June (06)July (07)Aug (08)Sep (09)Oct (10)
23725239739927616870282319
Table 2. Values of the parameters estimated and fitted against the real measles cases.
Table 2. Values of the parameters estimated and fitted against the real measles cases.
ParameterDescriptionSource
ϕ 260,479Estimated
α 1.253133 × 10 3 Estimated
ω 0.97Estimated
τ 1.60056 × 10 7 Fitted
μ 9.3408Fitted
β 9.2373 × 10 1 Fitted
δ 5.8306 × 10 1 Fitted
ρ 0Estimated
γ 5.8306 × 10 1 Fitted
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Guo, B.; Khan, A.; Din, A. Numerical Simulation of Nonlinear Stochastic Analysis for Measles Transmission: A Case Study of a Measles Epidemic in Pakistan. Fractal Fract. 2023, 7, 130. https://doi.org/10.3390/fractalfract7020130

AMA Style

Guo B, Khan A, Din A. Numerical Simulation of Nonlinear Stochastic Analysis for Measles Transmission: A Case Study of a Measles Epidemic in Pakistan. Fractal and Fractional. 2023; 7(2):130. https://doi.org/10.3390/fractalfract7020130

Chicago/Turabian Style

Guo, Bing, Asad Khan, and Anwarud Din. 2023. "Numerical Simulation of Nonlinear Stochastic Analysis for Measles Transmission: A Case Study of a Measles Epidemic in Pakistan" Fractal and Fractional 7, no. 2: 130. https://doi.org/10.3390/fractalfract7020130

APA Style

Guo, B., Khan, A., & Din, A. (2023). Numerical Simulation of Nonlinear Stochastic Analysis for Measles Transmission: A Case Study of a Measles Epidemic in Pakistan. Fractal and Fractional, 7(2), 130. https://doi.org/10.3390/fractalfract7020130

Article Metrics

Back to TopTop