Stationary distribution of an SIR epidemic model with three correlated Brownian motions and general Lévy measure

1.   Introduction

Mathematical epidemiology represents a vital and growing field of science in which analytical concepts, advanced approaches, and dynamical systems are applied to a wide variety of issues in ecological and real-life sciences ^[1]. Many epidemiological frameworks and processes are inherently probabilistic because the environmental perturbations influence the expansion of the infection and make it very complicated to prophesy the possible pandemic scenarios ^[2,3,4]. Thus, deterministic explanations, although capable of making eminently useful predictions and speculations, are not empirical and not sufficiently applicable ^{[5,6,7,8,9,10]}. There is, therefore, a pressing need to generate and develop a refined framework that takes into account stochasticity interventions, mainly when examining the attitude of a highly infectious disease such as Coronavirus and Monkeypox ^{[11,12,13,14,15,16]}.

To correctly model a large amount of noise with some discontinuities, Poisson processes are used, which are renowned for their ability to simulate the spread of infection in a small population size, in the case of medical and financial crises, or when implementing certain drugs and non-pharmaceutical interventions such as quarantine, vaccination and global immunization ^[17]. Recently, the authors of ^[18] developed a new analytical framework to survey the long run of a SIR model driven by independent jumps noises and a finite Lévy measure ^[19]. However, the Lévy noises associated with finite measures are characterized by their slight tails which are not able to model intense phenomena. The huge variations can produce an unforeseen increase in the density of individuals depending on the illness characteristics ^[20,21,22]. In the present research, we expose an alternative scope that takes into account a general Lévy measure (finite or infinite) and fully correlated Lévy noises ^[23]. For the purpose of comparison, we keep the same epidemic model studied in ^[18], but this time with the following representation:

where the positive constants $\mathfrak{A}$, $\mathfrak{B}$, $\mathfrak{C}$, $\mathfrak{D}$ and $\mathfrak{D}_{d}$ are respectively the flow into the susceptible class $S$, the prevalence rate, the permanent cure ratio (the transfer from the contaminated class $I$ to the permanently recovered compartment $R_p$), the demise rate and the disease-related mortality rate. Here and elsewhere, $S(t_{-})$, $I(t_{-})$ and $R_p(t_{-})$ are respectively the left limits of the functions $S(t)$, $I(t)$ and $R_p(t)$. The vector ${\mathbb A} = ({\mathbb A}_1, {\mathbb A}_2, {\mathbb A}_3)$ indicates the three-dimensional Lévy process with its associated Lévy-Khintchine formula:

where ${\mathbb E}$ is the mathematical expectation, $k = (k_1, k_2, k_3)\in{\mathbb R}^3$ and $\mathcal{K} = (\mathfrak{k}_{\ell, j})_{1\leq \ell, j\leq 3}$ represents a specific positive definite matrix. The jumps amplitude $\mathbf{e}_\ell:{\mathcal H}\subset{\mathbb R}^3_+\to {\mathbb R}$ ($\ell = 1, 2, 3$) are continuous functions. $\mathbf{v}_\text{Lévy}$ is the general Lévy measure (finite or infinite) such that

In accordance with the theory exhibited in ^[20] and ^[24], the process ${\mathbb A}$ takes the following form:

We consider a filtered probability space $(\Omega_{{\mathbb P}}, \mathcal{E}, \{\mathcal{E}_t\}_{t\geq 0 }, {\mathbb P})$ such that $\{\mathcal{E}_t\}_{t\geq 0}$ is filtration that verifies these hypotheses: right continuous, increasing and $\mathcal{E}_0$ includes all ${\mathbb P}$-null sets. In (1.2), ${\mathbb W}^{\mathfrak{k}} = \big({\mathbb W}^{\mathfrak{k}}_1, {\mathbb W}^{\mathfrak{k}}_2, {\mathbb W}^{\mathfrak{k}}_3\big)$ is referring to a Gaussian process with its co-variance matrix $\mathcal{K}$. $\mathcal{N}_{\mathbf{v}_\text{Lévy}}$ is a Poisson measure which is independent of ${\mathbb W}^{\mathfrak{k}}$ (natural hypothesis). ${\mathbb C}$ is the compensator process and $\mathbf{v}_\text{Lévy}$ is its corresponding Lévy measure, where ${\mathbb C}(t, {\text{d}} u) = \mathcal{N}_{\mathbf{v}_\text{Lévy}}(t, {\text{d}} u)-t\mathbf{v}_\text{Lévy}({\text{d}} u)$. The co-variances of ${\mathbb A}$ are as follows:

In ^[24], Privault and Wang obtained sufficient criteria for the disease vanishing and its insistence in the case of SIR model with the representation of (1.2). However, the existence of a unique stationary distribution has not been investigated due to some technical difficulties. It must be mentioned that the stationarity is an important statistical property for random processes. In this survey, we address this issue using a novel procedure different from that exhibited in ^[18].

To illustrate the importance of our work, we numerically treat an example of a general Lévy measure using the tempered stable processes. This class of Lévy distributions is widely applied in the case of an infinite Lévy measure ^[24]. According to ^[20], we consider the following tempered $\alpha$-stable Lévy measure:

where $\mathfrak{k}_{\alpha}(\cdot)$ is an infinite Lévy measure defined on ${\mathcal H}$ such that

Here, the measure $\mathfrak{k}_{\alpha}({\text{d}} r)$ can be decomposed as follows:

for all positive quantities $n_{-}, n_{+}, \varsigma_{-}, \varsigma_{+} > 0$, where $\delta_d$ indicates the Dirac mass measure at $d\in{\mathbb R}^3$. From (1.3), the infinite measure $\mathbf{v}_\text{Lévy}$ can be expressed as follows:

The remaining of this study is structured as follows: in Section 2, we prove the ergodic characteristic of the disturbed model (1.1) with the representation (1.2). In Section 3, we treat a numerical example to belay and emphasize the proposed approach. Furthermore, we explore the effect of tempered $\alpha$-stable quadratic Lévy noises on the long run attitude of the infection. In Section 4, we present the main conclusions of our article.

2.   Theoretical result: the existence of a unique stationary distribution

Before exhibiting the pivotal outcome of this work, we firstly present the hypothetical setting, some useful lemmas and some conventions to lighten its statement. For simplicity, it will be practical to use the following notations throughout the rest of the paper:

$\bullet$ $\mathbf{N_1} = { } \int_{{\mathcal H}}\mathbf{e}_\ell^2(u)\mathbf{v}_\text{Lévy}({\text{d}} u), \; \; \; \ell = 1, 2, 3$.

$\bullet$ $\mathbf{N_2} = { }\int_{{\mathcal H}}\Big(\ln\big(\mathbf{e}_\ell(u)+1\big)\Big)^2\mathbf{v}_\text{Lévy}({\text{d}} u), \; \; \; \ell = 1, 2, 3$.

$\bullet$ $\mathbf{N_3} = { }\int_{{\mathcal H}}\Big(\mathbf{e}_\ell(u)-\ln\big(\mathbf{e}_\ell(u)+1\big)\Big)\mathbf{v}_\text{Lévy}({\text{d}} u), \; \; \; \ell = 1, 2, 3$.

$\bullet$ $\mathbf{N_4} = { }\max\big\{\mathbf{e}_1(u), \mathbf{e}_2(u), \mathbf{e}_3(u)\big\}$.

$\bullet$ $\mathbf{N_5} = { } \max \big\{1, 2^{2p-3}\big\}(2p^2-p)\int_{{\mathcal H}}\Big\{\mathbf{N_4}^2+\mathbf{N_4}^{2p}\Big\}\mathbf{v}_\text{Lévy}({\text{d}} u)$.

$\bullet$ $\mathbf{N_6} = { } \mathfrak{D}-0.5(2p-1)\underset{i = 1, 2, 3}{\max}\sum\limits_{j = 1}^{3}|\mathfrak{k}_{i, j}|-\frac{0.5\mathbf{N_5}}{p}$.

According to some analytic and mathematical reasons, we presume that

$\bullet$ A1: $\mathbf{e}_i(u) > -1$ and the quantities $\mathbf{N_1}$, $\mathbf{N_2}$ and $\mathbf{N_3}$ are finite.

$\bullet$ A2: $\mathbf{N_6} > 0$ for some $p > 1$.

The next lemma emphasizes the mathematical well-posedness of system (1.1) with the perturbation (1.2).

Lemma 2.1 (Well-posedness of the model, ^[25]). Under the assumption A1, the probabilistic system (1.1) is well posed.

To proceed, we will establish some estimates of the total number of population $\mathbf{N}_{ \text{total}}(t) = S(t)+I(t)+R_p(t)$.

Lemma 2.2 (Moments estimates of $\mathbf{N}_{ \text{total}}(t)$). Let assumptions $\text{A1}$ and A2 hold. Then for any $p\in(1, \infty)$ such that $\mathbf{N_6} > 0$, we have

1) ${ }{\mathbb E}\Big\{(1+\mathbf{N}_{ \text{total}}(t))^{2p}\Big\}\leq \frac{2p\mathfrak{H}}{c}+\Big\{1+\mathbf{N}_{ \text{total}}(0)\Big\}^{2p}e^{-ct}$,

2) ${ }\underset{t\to \infty}{\lim \sup} \; \frac{1}{t}\int^t_0{\mathbb E}\Big\{(1+\mathbf{N}_{ \text{total}}(s))^{2p}\Big\}{\text{d}} s\leq \frac{2p\mathfrak{H}}{c}$,

where $c\in(0, 2p\mathbf{N_6})$ and ${ }\mathfrak{H} = \underset{z\in{\mathbb R}_+}{\sup}\bigg\{z^{2p-2}\bigg(-\Big(\mathbf{N_6}-\frac{c}{2p}\Big)z^2+\Big(\mathfrak{A}-\mathfrak{D}+\frac{c}{p}\Big)z+\mathfrak{A}+\frac{c}{2p}\bigg)\bigg\}+1$.

Proof. Set $\mathbf{N}_{ \text{total}, \textbf{e}}(t, u) = \textbf{e}_1(u)S(t)+\textbf{e}_2(u)I(t)+\textbf{e}_3(u)R_p(t)$ for all $u\in {\mathcal H}$. Making use of Itô's lemma to the function $\mathcal{V}(x) = (1+x)^{2p}$, $p > 1$, we get

where $\mathcal{L}\mathcal{V}(\mathbf{N}_{ \text{total}})$ is defined as follows:

For simplicity, we define

For all $u\in{\mathcal H}$ and $t > 0$, there is $0 < \varepsilon < 1$ such that

where $\mathfrak{c}_{2p} = \max \big(1, 2^{2p-3}\big)(2p^2-p)$. Then, we get

where $\mathbf{N_6}$ is positive due to $\mathbf{A_2}$. Now, we choose any $c\in(0, 2p\mathbf{N_6})$, then

Then, by taking the integration and the mathematical expectation of (2.2), we have

Therefore, we get

Obviously, we obtain

This completes the proof.

Theorem 2.1. Assume that $\mathbf{A1}$ and $\mathbf{A2}$ hold. If

then, a unique stationary distribution exists for the probabilistic model (1.1).

Proof. Consider the following initial value problem with the presentation (1.2):

Let $f(t) = -\frac{\mathfrak{B}}{\mathfrak{D}}\big(-S(t)+\Xi(t)\big)+\ln I(t)$ and $\mathfrak{X}(t) = S(s)I(s)$. The application of Itô's lemma implies that

An integration of (2.4) gives

According to Lemma 2.2 of ^[24] and $\mathbf{A1}$, we obtain

Now, we consider $\mathbf{E}_1 = \big\{(t, \omega)\in{\mathbb R}_+\times\Omega_{{\mathbb P}}| \quad S(t, \omega)\geq \chi, \quad \mbox{and}, \quad I(t, \omega)\geq \chi\big\}$, $\mathbf{E}_2 = \big\{(t, \omega)\in{\mathbb R}_+\times\Omega_{{\mathbb P}}| \quad S(t, \omega)\leq \chi\big\}$, and $\mathbf{E}_3 = \big\{(t, \omega)\in{\mathbb R}_+\times\Omega_{{\mathbb P}}| \quad I(t, \omega)\leq \chi\big\}$, where $\chi > 0$ will be chosen later. Then

We can choose ${ }\chi \leq \frac{\mathfrak{D}^2}{4\mathfrak{B}^2\mathfrak{A}}(\mathfrak{D}+\mathfrak{D}_{d}+\mathfrak{C})\big(\mathcal{T}_{\circ, s}-1\big)$, therefore

Let $p\in(1, \infty)$ such that $\mathbf{N_6} > 0$ and $q$ is given by $\frac{1}{q} = 1-\frac{1}{p}$. By using Young inequality, we obtain

where $\zeta > 0$ verifies $\zeta^p\leq\frac{\mathfrak{D}c\mathfrak{B}^{-(p+1)}}{8 \mathfrak{H}}(\mathfrak{D}+\mathfrak{D}_{d}+\mathfrak{C})\big(\mathcal{T}_{\circ, s}-1\big)$. From (2.5) and Remark 2, we have

Consider

where $v > \chi > 0$ will be defined in the next. Now, by employing Markov's inequality, we get

We choose ${ }\frac{1}{v}\leq \frac{\mathfrak{D}^2 q\zeta^q}{8\mathfrak{B} \mathfrak{A}}(\mathfrak{D}+\mathfrak{D}_{d}+\mathfrak{C})\big(\mathcal{T}_{\circ, s}-1\big)$, then we obtain

Via (2.6), we conclude that

In conclusion, we found the compact subset ${\bf D}\subset{\mathbb R}^3_+$ that verifies

where $x_0 = (S(0), I(0), R_p(0))$. Identical to the demonstration of (Lemma 3.2., ^[26]), we establish that the unique solution of the model (1.1) has the Feller property. According to the mutually limited possibilities lemma ^[27], a unique ergodic stable distribution exists for (1.1) with the representation (1.2).

3.   Numerical application: SIR system with tempered Lévy noise and infinite measure

This part of the manuscript is dedicated to the presentation of numerical examples and to the verification of the results of Theorem 2.1. Using computer simulations, we exhibit some plots for the paths and their corresponding histograms, from which the complex dynamic behaviors of the probabilistic system (1.1) can be easily explored. Additionally, we select some simulated parameter values to support our theoretical framework. Henceforth, the units adopted for time and number of individuals are one day and one million population.

In accordance with the study presented in ^[20], we consider a compensated tempered Poisson process of the following form:

with the well defined infinite Lévy measure (1.4). We suppose that the vector ${\mathbb W}^{\mathfrak{k}}$ is given as follows: ${\mathbb W}^{\mathfrak{k}}_1 = \mathfrak{k}_{1, 1}\mathcal{W}_a$, ${\mathbb W}^{\mathfrak{k}}_2 = \mathfrak{k}_{2, 1}\mathcal{W}_a+\mathfrak{k}_{2, 2}\mathcal{W}_b$, ${\mathbb W}^{\mathfrak{k}}_3 = \mathfrak{k}_{3, 1}\mathcal{W}_a+\mathfrak{k}_{3, 2}\mathcal{W}_b+\mathfrak{k}_{3, 3}\mathcal{W}_c$, where $\mathcal{W}_a$, $\mathcal{W}_b$ and $\mathcal{W}_c$ are three independent Brownian motions.

For simplicity, we choose $\mathbf{e}_\ell(u) = e_\ell u$ $(\ell = 1, 2, 3)$, where $e_h > 0$, and we deal with the following probabilistic model:

Remark 3.1. Indeed, we mention that the assumption A1 is naturally verified and the quantity $\mathbf{N5}$ is finite when $p > \alpha$.

By applying the algorithm proposed in ^[20], we will exhibit some numerical illustrations in the case of the unique-sided tempered process ${\mathbb Y}(t)$ with $n_{-} = 0$. So, we select $\alpha = 0.7$, $n_{+} = 2.8$, $\varsigma_{-} = \varsigma_{+} = 1.2$, $\mathfrak{k}_{1, 1} = 0.2$, $\mathfrak{k}_{2, 1} = 0.16$, $\mathfrak{k}_{3, 1} = 0.15$, $\mathfrak{k}_{2, 2} = 0.12$, $\mathfrak{k}_{3, 2} = 0.12$, $\mathfrak{k}_{3, 3} = 0.1$, $e_1 = 0.2$, $e_2 = 0.8$ and $e_3 = 0.5$.

For the stochastic system (3.2), the biological parameters are taken as follows: $\mathfrak{A} = 8, \; \; \mathfrak{B} = 5.1, \; \; \mathfrak{D} = 5.3\; \; \mathfrak{D}_{d} = 0.5, \; \mathfrak{C} = 0.7.$ Then, A2 holds and $\mathcal{T}_{\circ, s} = 1.0948 > 1$. In line with Theorem 2.1, we infer the existence of a unique stationary distribution for each class which is exactly illustrated in Figure 1. In the same Figure, we show the permanence of all trajectories. For a complete overview of the marginal densities of the solution, we plot the two-dimensional empirical distribution in Figure 2.

Figure 1.  First column presents the frequency histogram fitting curves at time $t = 1000$ and the associated density functions, respectively. Second column presents the trajectories of the stochastic solution.

DownLoad: Full-Size Img PowerPoint

Figure 2.  The $3D$ graph of the joint two-dimensional densities. Different colors represent different sizes of the individuals concentration.

DownLoad: Full-Size Img PowerPoint

Now, we select $\mathfrak{A} = 8, \; \; \mathfrak{B} = 4.8, \; \; \mathfrak{D} = 5.3\; \; \mathfrak{D}_{d} = 0.5, \; \mathfrak{C} = 0.7.$ Based on the results demonstrated in ^[24], we conclude that the illness will almost certainly go away because we have $\mathcal{T}_{\circ, s} = 0.9799 < 1$. To probe the effect of Lévy jumps in this case, we compare the solution of (3.2) with deterministic solutions and free jumps. From Figure 3, we observe that the noise of the jumps leads to the extinction of the infection whereas the non-stochastic solution and the paths without leaps both persist. This means that jumps have a positive effect on the pandemic situation, and noise with infinite measure can change the long-term behavior of the dynamical system.

Figure 3.  Computer simulation of the trajectories of system (3.2), the deterministic solution and the stochastic trajectories with only white noise.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

By considering correlated noise items and an infinite Lévy measure, in this research, we have exhibited an analytical and numerical framework to explore the dynamics of our perturbed epidemic model. Explicitly, we have studied the ergodic property by employing the lemma of mutually limited possibilities and some analytical tools. Ergodicity indicates that the epidemic will prevail and persist in long-term evolution. In the numerical simulations part, we have ensured the accuracy of our threshold. Further, we explored the impact of noises and heavy tails on the infection dynamics. In particular, we proved that jumps have a negative influence on the long-term behavior of the illness in the sense that they lead to complete extinction.

Generally, we pointed out that this research upgrades many previous papers to the case of general Lévy tails. In addition, it presents new insights into understanding illness spread with complex real-world hypotheses. In other words, the method described in this paper opens up many possibilities for future research.

Some fascinating topics deserve more attention. For example, we can consider our model with fractal-fractional differentiation ^{[28,29,30,31]}. This framework is an attractive branch of applied mathematics that deals with derivatives and integrals of non-integer order. Due to its amazing features, it is preferred for describing and simulating real-world problems in various fields such as biological mechanisms, material science, hydrological modeling, and economic phenomena ^{[32,33,34,35,36]}. We will address this idea in our future work.

Acknowledgments

This research was supported by Taif University Researchers Supporting Project number (TURSP2020/144), Taif University, Taif, Saudi Arabia.

Conflict of interest

The authors state that there is no conflict of interest.