# Inverse projection of axisymmetric orientation distributions

Philipp A. Kloza  James A. Elliott<EMAIL_ADDRESS>Department of Materials
Science & Metallurgy
University of Cambridge
27 Charles Babbage Road, Cambridge, United Kingdom

###### Abstract

We show that the projection of an axisymmetric three-dimensional orientation
distribution to two dimensions can be cast into an Abel transform. Based on
this correspondence, we derive an exact integral inverse, which allows for the
quantification of three-dimensional uniaxial alignment of rod-like units from
two-dimensional sliced images, thus providing an alternative to X-ray or
tomographic analysis. A matrix representation of the projection and its
inverse is derived, providing a direct relationship between two- and three-
dimensional order parameters for both polar and non-polar systems.

††preprint: APS/123-QED

##  I Introduction

The three-dimensional orientation distribution of synthetic and naturally
occurring fibrillar materials can have a major influence on their effective
macroscopic properties. For instance, aligning cellulose fibres can
substantially improve their stiffness, tensile strength and thermal stability
[1, 2]. Likewise, aligning carbon nanotube fibres has been shown to greatly
improve their mechanical properties as well as their electrical and thermal
conductivity [3, 4, 5, 6, 7], which was recently confirmed by a large meta-
analysis [8]. In biological systems, aligning collagen fibres also increases
their stiffness [9, 10] and is directly related to the mechanical strength of
tendons [11, 12] and bones [13]. Changes in collagen alignment have also been
identified as a predictor for breast [14, 15] and pancreatic [16] cancer, and
progressive organ diseases such as liver fibrosis [17].

Naturally, there is an interest in quantifying the alignment of the orientable
units in fibrillar structures. If the orientation distribution is axisymmetric
in three dimensions, i.e. if it can be described as a one-dimensional function
of the alignment angle with a specified axis, the alignment is often
quantified using mean values of Legendre polynomials, or order parameters,
$\langle P_{l}\rangle$ [18], where the value for $l=2$ is more commonly known
as the Hermans order parameter [19].

Traditionally, wide- or small-angle X-ray diffraction (XRD) is used to measure
the orientation distribution function (ODF) [20, 21, 22], with the first
theoretical considerations going back to the 1930s [23, 19]. Fibre orientation
can also be studied by direct three-dimensional imaging techniques such as
computed tomography [24]. However, it is more common to employ image analysis
methods based on optical or scanning electron micrograph (SEM) images using
the Fourier transform [25, 26, 27, 28], the Hough transform [29, 30], direct
fibre tracking [27, 31] or topological data analysis [32] as a fast and facile
alternative to a full three-dimensional analysis.

As image data is two-dimensional, the ability to extract true three-
dimensional data is limited, firstly because the observed structure may not be
axisymmetric in the first place, as is often the case with film materials
[25]. Secondly, even if the structure is axisymmetric, e.g. it is a CNT [33]
or cellulose [2] fibre or a liquid crystal [34], any image shows merely a
projection of the structure. Likewise, the three-dimensional ODF is projected
onto two dimensions [20]. In both cases, using the three-dimensional order
parameters to quantify alignment is inappropriate as the ODF is not truly
three-dimensional. Incorrectly treating the two-dimensional ODF as a three-
dimensional one and using it to calculate three-dimensional order parameters
will lead to values which deviate from other experimental techniques such as
XRD [20, 32]. Rather, the ODF should be treated as two-dimensional and mean
values of the Chebyshev polynomials $\langle T_{m}\rangle$ should be used as
order parameters instead [25, 28, 35]. Subsequently, rigorously extracting
information about the three-dimensional axisymmetric ODF from the two-
dimensional one obtained from image data would allow for the direct comparison
of the measurement to ODFs obtained from XRD experiments and is the main aim
of the method presented here.

In this work, we show that the projection applied to the three-dimensional ODF
$f(\theta_{\mathrm{3D}})$ to obtain the two-dimensional ODF
$g(\theta_{\mathrm{2D}})$ can be formulated as an integral transform using the
Leadbetter-Norris (LN) kernel [36], which can also be interpreted as an Abel
transform. Deutsch showed that this particular Abel transform can be inverted
to infer the three-dimensional ODF from XRD experiments [37]. We extend
Deutsch’s derivation of the exact inverse of the projection to account for
orientable structures and show that it can be used to retrieve the true three-
dimensional ODF from two-dimensional data. Additionally, we derive a matrix
representation of the transform and its inverse in the basis of Chebyshev and
Legendre polynomials which provides a direct and straightforward relationship
between two- and three-dimensional order parameters. The results of this work
can be readily implemented as an additional processing step after obtaining
the two-dimensional ODF from image data. The relationships and differences
between the orientation distribution functions and the structures they are
derived from are summarised in fig. 1, which justifies the need for a rigorous
inverse projection.

The methodology derived in this work is applicable to any projected, two-
dimensional ODF obtained of an axisymmetric three-dimensional structure and is
fully agnostic of the method used to obtain the two-dimensional ODF in the
first place. Hence, it may be readily implemented as an additional simple
processing step in a multitude of different software tools [28, 27, 38].

Figure 1: Illustration of projection methods derived in this work. A two-
dimensional image is derived from a three-dimensional structure. Two- and
three-dimensional ODFs can be computed from the structures and converted into
each other using the direct and inverse projection.

##  II Problem Definition

Suppose we have a three-dimensional unit vector $\hat{r}$ parametrised in
spherical coordinates:

$\hat{r}=\begin{pmatrix}\sin{\theta_{\mathrm{3D}}}\cos{\phi}\\\
\sin{\theta_{\mathrm{3D}}}\sin{\phi}\\\
\cos{\theta_{\mathrm{3D}}}\end{pmatrix}.$ (1)

Further suppose this vector is distributed according to a three-dimensional
ODF in spherical coordinates $f(\theta_{\mathrm{3D}})$. The probability
measure of the three-dimensional ODF is:

$\mathrm{d}\Omega_{\mathrm{3D}}\equiv\frac{1}{2\pi}\mathrm{d}\phi\,\mathrm{d}\theta_{\mathrm{3D}}\,\sin\theta_{\mathrm{3D}}f(\theta_{\mathrm{3D}}),$
(2)

where the factor of $1/2\pi$ is due to axisymmetry. Assuming the plane of
projection contains the axis of symmetry (the $x_{3}$-axis), we want to derive
the two-dimensional ODF $g(\theta_{\mathrm{2D}})$ which can be measured from
image data. The two-dimensional orientation angle $\theta_{\mathrm{2D}}$ is
defined by:

$\cos{\theta_{\mathrm{2D}}}=\frac{r_{3}}{\sqrt{r_{2}^{2}+r_{3}^{2}}}=\frac{\cos{\theta_{\mathrm{3D}}}}{\sqrt{\cos^{2}{\theta_{\mathrm{3D}}}+\sin^{2}{\theta_{\mathrm{3D}}}\sin^{2}{\phi}}},$
(3)

which also establishes all relations between $\phi$, $\theta_{\mathrm{2D}}$
and $\theta_{\mathrm{3D}}$ relevant in this work. Now, we may define the two-
dimensional probability measure via the two-dimensional ODF:

$\mathrm{d}\Omega_{\mathrm{2D}}\equiv\mathrm{d}\theta_{\mathrm{2D}}\,g(\theta_{\mathrm{2D}}).$
(4)

##  III Projection as an Abel Transform and its Inverse

We may transform from the three-dimensional probability measure to the two-
dimensional one by an appropriate change of variables from
$\theta_{\mathrm{3D}}$ to $\theta_{\mathrm{2D}}$ and by integrating out the
contribution of the azimuthal angle $\phi$. For simplicity, we will restrict
$\theta_{\mathrm{2D}}$ to the range $[0,\pi]$ and extend the 2D ODF to the
full range $[0,2\pi)$ by symmetry about $\phi=\pi$. We also note the we do not
consider symmetry of either ODF about
$\theta_{\mathrm{2D}},\theta_{\mathrm{3D}}=\pi/2$ as it not required for the
methodology shown in this work. This somewhat more general result may be
applied to measure the alignment of truly directional fibres, e.g. if they are
polar, or for the measurement of spin texture.

According to eq. (3), the differential of $\theta_{\mathrm{3D}}$ transforms
as:

$\mathrm{d}\theta_{\mathrm{3D}}=\mathrm{d}\theta_{\mathrm{2D}}\left|\frac{\partial\theta_{\mathrm{3D}}}{\partial\theta_{\mathrm{2D}}}\right|=\mathrm{d}\theta_{\mathrm{2D}}\frac{\sin{\phi}}{1-\cos^{2}{\theta_{\mathrm{2D}}}\cos^{2}{\phi}},$
(5)

so the two-dimensional ODF can be explicitly computed as follows:

$g(\theta_{\mathrm{2D}})=\frac{1}{\pi}\int_{0}^{\pi/2}\mathrm{d}\phi\,\left|\frac{\partial\theta_{\mathrm{3D}}}{\partial\theta_{\mathrm{2D}}}\right|\sin\theta_{\mathrm{3D}}f(\theta_{\mathrm{3D}}),$
(6)

where $\theta_{\mathrm{3D}}=\theta_{\mathrm{3D}}(\phi,\theta_{\mathrm{2D}})$
is implied. The above expression constitutes the projection of the three-
dimensional ODF to two dimensions and further accounts for the symmetry of the
contributions about $\phi=\pi/2$.

We may cast the integral in eq. (6) into a variant of the Abel transform by a
change of integration variable from $\phi$ to $t=\cos{\theta_{\mathrm{3D}}}$.
The change in integral measure is again computed by making use of eq. (3),
solving for $\phi$ and differentiating with respect to $t$. After simplifying
the resulting trigonometric functions, the integral is:

$g(\theta_{\mathrm{2D}})=\frac{1}{\pi|\cos{\theta_{\mathrm{2D}}}|}\int\displaylimits_{0}^{\cos{\theta_{\mathrm{2D}}}}\mathrm{d}t\,\frac{tf(\arccos{t})}{\sqrt{\cos^{2}{\theta_{\mathrm{2D}}}-t^{2}}}.$
(7)

This integral transform is an alternative formulation of a transform using the
Leadbetter-Norris (LN) kernel which was used to calculate XRD intensity
distributions of axisymmetric liquid crystal systems [36]. Ultimately, the LN
kernel proved to be incorrect for the use in XRD [39, 40], yet it is still
applicable in our context of projected ODFs. Hence, we may follow Deutsch’s
derivation of the inverse of this transform [37], where we also account for
the possibility of orientable structures. We may further substitute
$\tau=t^{2}$ and $\xi=\cos^{2}{\theta_{\mathrm{2D}}}$ in the integral in eq.
(7):

$g(\arccos(\pm\sqrt{\xi}))=\frac{1}{2\pi\sqrt{\xi}}\int_{0}^{\xi}\mathrm{d}\tau\,\frac{f(\arccos(\pm\sqrt{\tau}))}{\sqrt{\xi-\tau}}.$
(8)

In eq. 8, the $\pm$ symbol indicates two branches with equal signs.
Furthermore, the integral on the right-hand side of the equation is an Abel
transform, and we have the following inverse [41]:

$f(\arccos(\pm\sqrt{\tau}))=2\frac{\mathrm{d}}{\mathrm{d}\tau}\int\displaylimits_{0}^{\tau}\mathrm{d}\xi\,\frac{\sqrt{\xi}g(\arccos(\pm\sqrt{\xi}))}{\sqrt{\tau-\xi}}.$
(9)

Finally, we re-substitute $\tau=\cos^{2}{\theta_{\mathrm{3D}}}$ and
$\xi=\cos^{2}{\theta_{\mathrm{2D}}}$ to find the exact inverse projection:

$\displaystyle f(\theta_{\mathrm{3D}})=$
$\displaystyle\frac{2}{\sin{(2\theta_{\mathrm{3D}})}}\frac{\mathrm{d}}{\mathrm{d}\theta_{\mathrm{3D}}}$
(10)
$\displaystyle\times\int\displaylimits_{\pi/2}^{\theta_{\mathrm{3D}}}\mathrm{d}\theta_{\mathrm{2D}}\,\frac{|\cos{\theta_{\mathrm{2D}}}|\sin{(2\theta_{\mathrm{2D}})}g(\theta_{\mathrm{2D}})}{\sqrt{\cos^{2}{\theta_{\mathrm{3D}}}-\cos^{2}{\theta_{\mathrm{2D}}}}}.$

While this is insightful on its own, numerical evaluation of expression (10)
can be difficult due to the singularities both in the integral and
multiplicative factor. We will therefore derive representations of the
projection and its inverse in terms of moment expansions which allow us to
avoid using the integral transforms. The derivation of the projection in the
ODF moment bases follows in the next sections.

##  IV Legendre Expansion

In order to compute the relationship between the moments of the two- and
three-dimensional ODFs, we expand the three-dimensional ODF in terms of
Legendre polynomials $P_{l}(\cos{\theta}_{\mathrm{3D}})$:

$f(\theta_{\mathrm{3D}})=\sum_{l=0}^{\infty}\frac{2l+1}{2}\langle P_{l}\rangle
P_{l}(\cos{\theta}_{\mathrm{3D}}),$ (11)

where the coefficients are a set of three-dimensional order parameters given
by:

$\langle
P_{l}\rangle=\int_{0}^{\pi}\mathrm{d}\theta_{\mathrm{3D}}\,\sin{\theta_{\mathrm{3D}}}f(\theta_{\mathrm{3D}})P_{l}(\cos{\theta}_{\mathrm{3D}}).$
(12)

By expanding the Legendre polynomials $P_{l}(\cos{\theta_{\mathrm{3D}}})$
using the Rodrigues formula [42], the two-dimensional ODF can be written as:

$\displaystyle g(\theta_{\mathrm{2D}})=$
$\displaystyle\sum_{l=0}^{\infty}\frac{2l+1}{2\pi}\frac{\langle
P_{l}\rangle}{2^{l}}$ $\displaystyle\times\sum_{k=0}^{\left\lfloor
l/2\right\rfloor}(-1)^{k}\binom{l}{k}\binom{2l-2k}{l}I_{l-2k}(\theta_{\mathrm{2D}}),$

where $I_{n}(\theta_{\mathrm{2D}})$ denotes the following integral:

$I_{n}(\theta_{\mathrm{2D}})=\int_{0}^{\pi/2}\mathrm{d}\phi\,\left|\frac{\partial\theta_{\mathrm{3D}}}{\partial\theta_{\mathrm{2D}}}\right|\sin\theta_{\mathrm{3D}}\cos^{n}{\theta_{\mathrm{3D}}}.$
(14)

By solving eq. (3) for $\theta_{\mathrm{3D}}$, the integral is:

$\displaystyle I_{n}(\theta_{\mathrm{2D}})=$
$\displaystyle\int_{0}^{\pi/2}\mathrm{d}\phi\,\frac{\sin^{n+1}{\phi}\sin{\theta_{\mathrm{2D}}}\cos^{n}{\theta_{\mathrm{2D}}}}{\left(1-\cos^{2}{\theta_{\mathrm{2D}}}\cos^{2}{\phi}\right)^{(n+3)/2}}$
$\displaystyle=$
$\displaystyle\frac{\sqrt{\pi}}{2}\frac{\Gamma\left(\frac{n+2}{2}\right)}{\Gamma\left(\frac{n+3}{2}\right)}\cos^{n}\theta_{\mathrm{2D}},$
(15)

where $\Gamma$ is the Euler gamma function.

##  V Chebyshev Expansion

The next step in finding a matrix representation of the projection is to
expand the two-dimensional ODF in an appropriate basis. A natural choice is
the Fourier series due to periodicity, or equivalently an expansion in terms
of Chebyshev polynomials $T_{m}$:

$g(\theta_{\mathrm{2D}})=\frac{\langle
T_{0}\rangle}{2\pi}+\sum_{m=1}^{\infty}\frac{\langle
T_{m}\rangle}{\pi}T_{m}(\cos\theta_{\mathrm{2D}}),$ (16)

where the equivalence to the Fourier series is given due to
$T_{m}(\cos\theta_{\mathrm{2D}})=\cos(m\theta_{\mathrm{2D}})$. By symmetry of
the ODF, the function is even and we only have cosine terms in the expansion.
The coefficients are given by:

$\langle
T_{m}\rangle=\int_{0}^{2\pi}\mathrm{d}\theta_{\mathrm{2D}}\,T_{m}(\cos\theta_{\mathrm{2D}})g(\theta_{\mathrm{2D}}).$
(17)

Due to the linearity of the Abel transform, there exists a matrix
representation of the projection $\mathcal{P}_{ml}$ between the coefficients
$\langle T_{m}\rangle$ and $\langle P_{l}\rangle$ of the respective two- and
three-dimensional ODFs:

$\langle T_{m}\rangle=\sum_{l=0}^{\infty}\mathcal{P}_{ml}\langle
P_{l}\rangle.$ (18)

The matrix element $\mathcal{P}_{ml}$ is obtained by calculating the
coefficients $\langle T_{m}\rangle$ of the expansion given in eq. (IV) and
extracting the contribution due to $\langle P_{l}\rangle$:

$\displaystyle\mathcal{P}_{ml}=$
$\displaystyle\frac{2l+1}{2\pi}\frac{1}{2^{l}}\sum_{k=0}^{\left\lfloor
l/2\right\rfloor}(-1)^{k}\binom{l}{k}\binom{2l-2k}{l}$
$\displaystyle\times\int_{0}^{2\pi}\mathrm{d}\theta_{\mathrm{2D}}\,\cos(m\theta_{\mathrm{2D}})I_{l-2k}(\theta_{\mathrm{2D}}).$

The integral in the above equation can be solved using the trigonometric power
formula [43] and the orthogonality of trigonometric functions. We may also
derive the coefficients of the inverse matrix using the same approach as above
by using the Chebyshev expansion of the two-dimensional ODF from eq. (IV) in
the explicit inverse transform in eq. (10) and then calculating the
contributions of $\langle T_{m}\rangle$ to $\langle P_{l}\rangle$. A detailed
derivation can be found in the Supplemental Material [44]. We find for the
matrix coefficients of the projection and its inverse:

$\displaystyle\mathcal{P}_{ml}=$
$\displaystyle\begin{dcases}\frac{\sqrt{\pi}}{2}\frac{2l+1}{4^{l}}\sum_{k=0}^{(l-m)/2}(-4)^{k}\binom{l}{k}\binom{2l-2k}{l}\binom{l-2k}{\frac{l-m-2k}{2}}\frac{\Gamma\left(\frac{l-2k+2}{2}\right)}{\Gamma\left(\frac{l-2k+3}{2}\right)},&l\geq
m\text{ and }l-m\text{ even,}\\\ 0,&\text{ otherwise},\end{dcases}$ (20)
$\displaystyle\mathcal{P}_{lm}^{-1}=$ $\displaystyle\begin{dcases}1,&l=m=0,\\\
\frac{m}{2\sqrt{\pi}}\sum_{k=0}^{\lfloor
m/2\rfloor}(-1)^{k}2^{m-2k}(2+m-2k)\frac{(m-k-1)!}{k!(m-2k)!}\frac{\Gamma\left(\frac{m-2k+3}{2}\right)}{\Gamma\left(\frac{m-2k+4}{2}\right)}L_{l,m-2k},&m>0,\\\
0,&\text{ otherwise},\end{dcases}$ (21)

where an additional matrix element $L_{mn}$ is given by:

$L_{mn}=\frac{1}{2^{m}}\sum_{k=0}^{\lfloor
m/2\rfloor}(-1)^{k}\binom{m}{k}\binom{2m-2k}{m}\frac{1+(-1)^{m+n}}{1+m+n-2k}.$
(22)

Using the matrix representation of the projection and its inverse, the first
even moments of the ODFs can be expressed as:

$\displaystyle\langle T_{2}\rangle=$ $\displaystyle\frac{5}{4}\langle
P_{2}\rangle-\frac{3}{8}\langle P_{4}\rangle+\frac{13}{64}\langle
P_{6}\rangle-\cdots,$ (23) $\displaystyle\langle P_{2}\rangle=$
$\displaystyle\frac{4}{5}\langle T_{2}\rangle+\frac{8}{35}\langle
T_{4}\rangle-\frac{4}{105}\langle T_{6}\rangle+\cdots.$ (24)

In particular, eq. (24) establishes a simple algebraic relation between the
Hermans order parameter and two-dimensional order parameters.

##  VI Numerical Example

Figure 2: ODFs in two and three dimensions, both computed directly (solid) and inferred from the ODF in the respectively other dimension (dashed). Treating the two-dimensional ODF as a three-dimensional one yields the dashed-dotted line for reference. | 2D Order Parameters | 3D Order Parameters
---|---|---
$n$ | $\langle T_{n}\rangle$ | $\mathcal{P}_{nm}\langle P_{m}\rangle$ | $\langle P_{n}\rangle$ | $\mathcal{P}_{nm}^{-1}\langle T_{m}\rangle$ | $\langle\bcancel{P}_{n}\rangle$
0 | 1. | 1. | 1. | 1. | 1.
1 | 0.245151 | 0.244988 | 0.233294 | 0.233466 | 0.243665
2 | 0.519902 | 0.520364 | 0.473539 | 0.473091 | 0.636382
3 | 0.161026 | 0.160886 | 0.147580 | 0.147576 | 0.190669
4 | 0.269882 | 0.270380 | 0.241801 | 0.241137 | 0.445402
5 | 0.095817 | 0.096274 | 0.087147 | 0.086754 | 0.146383
Table 1:  Order parameters obtained directly from two- and three-dimensional
ODFs and by inferring from respectively other dimension as well as by treating
the two-dimensional ODF as a three-dimensional one. Einstein summation
convention is used for brevity. The inverse projection from two to three
dimensions is highlighted in grey as the central result of this Letter.

To illustrate how the conversion from a two-dimensional to a three-dimensional
ODF and vice-versa works, we have generated a random set of $10^{6}$ vectors
according to a normal distribution $\mathcal{N}(0,1)$ for the $x_{1,2}$
components and $\mathcal{N}(1,3)$ for the $x_{3}$ component. The vectors were
normalised and their projection onto the $x_{2}$-$x_{3}$ plane was computed.
Then, the corresponding ODFs $f(\theta_{\mathrm{3D}})$ and
$g(\theta_{\mathrm{2D}})$ were computed by statistical binning. For the
projection, we used 21 Chebyshev and Legendre terms in the series expansions
of the ODFs combined with the matrix representation of both transforms. fig. 2
shows that the lines of the true ODFs and the respective ODFs obtained by
projection and inverse projection coincide visually, showing the validity of
the inverse projection derived in this work. Additionally, the figure also
shows a false ODF $\bcancel{f}(\theta_{\mathrm{3D}})$ which is obtained by
naively treating the two-dimensional ODF as three-dimensional with the measure
given in eq. (5). It clearly deviates from any other ODF, especially close to
the edges of the domain, showing that this treatment is indeed incorrect.
Table 1 shows that the order parameters in two and three dimensions are
correctly calculated by using the appropriate (inverse) projection matrix to
within $0.5\%$ of the values obtained directly from the true ODFs. It further
illustrates how using $\bcancel{f}(\theta_{\mathrm{3D}})$ to calculate the
moments $\langle\bcancel{P}_{l}\rangle$ leads to systematically incorrect
order parameters compared to the true three-dimensional ODF. The false ODF
substantially overestimates alignment (even order parameters) and polarity
(odd order parameters), further emphasising that a correct treatment of the
inverse projection is crucial when analysing two-dimensional data.

Figure 3: Convergence study showing absolute error of two- and three-
dimensional order parameters vs number of terms included in projection series
expansion.

fig. 3 shows that the three-dimensional order parameters converge faster
compared to the two-dimensional ones based on the series expansion approach
given in eqs. (20-21). For the three-dimensional order parameters, two terms
in the series expansion are sufficient to get within $1\%$ of the converged
value, supporting the fast convergence of the expansion formulation of the
inverse projection.

##  VII Conclusions

We have demonstrated that the projection of a three-dimensional ODF to a two-
dimensional ODF can be formulated as an Abel transform. Hence, we derived an
exact integral transform which recovers the true three-dimensional ODF from
two-dimensional data. By using a series expansion approach in Chebyshev and
Legendre polynomials for the two- and three-dimensional ODFs respectively, we
calculated matrix representations of the direct and inverse projection.
Numerical simulations show that our method correctly reconstructs both the
three-dimensional order parameters and ODF from two-dimensional data, while
naively treating the two-dimensional ODF as a three-dimensional one
overestimated the order parameters by up to 85%. The simplicity of the matrix
representation allows for the straightforward implementation of the inverse
projection in existing image analysis tools. By establishing a mathematically
rigorous footing for the reconstruction of three-dimensional alignment
information from two-dimensional optical or SEM images, our inverse projection
method provides a facile alternative to more involved experimental methods
such as XRD or tomography.

## References

  * Wang and Chen [2014] Y. Wang and L. Chen, Cellulose nanowhiskers and fiber alignment greatly improve mechanical properties of electrospun prolamin protein fibers, ACS Applied Materials & Interfaces 6, 1709 (2014), pMID: 24387200, https://doi.org/10.1021/am404624z .
  * Li  _et al._ [2021] K. Li, C. M. Clarkson, L. Wang, Y. Liu, M. Lamm, Z. Pang, Y. Zhou, J. Qian, M. Tajvidi, D. J. Gardner, H. Tekinalp, L. Hu, T. Li, A. J. Ragauskas, J. P. Youngblood, and S. Ozcan, Alignment of cellulose nanofibers: Harnessing nanoscale properties to macroscale benefits, ACS Nano 15, 3646 (2021), pMID: 33599500, https://doi.org/10.1021/acsnano.0c07613 .
  * Fischer  _et al._ [2003] J. Fischer, W. Zhou, J. Vavro, M. Llaguno, C. Guthy, R. Haggenmueller, M. Casavant, D. Walters, and R. Smalley, Magnetically aligned single wall carbon nanotube films: Preferred orientation and anisotropic transport properties, Journal of Applied Physics 93 (2003).
  * He  _et al._ [2016] X. He, W. Gao, L. Xie, B. Li, Q. Zhang, S. Lei, J. Robinson, E. Haroz, S. Doorn, W. Weipeng, R. Vajtai, P. Ajayan, W. Adams, R. Hauge, and J. Kono, Wafer-scale monodomain films of spontaneously aligned single-walled carbon nanotubes, Nature Nanotechnology 11 (2016).
  * Yamaguchi  _et al._ [2019] S. Yamaguchi, I. Tsunekawa, N. Komatsu, W. Gao, T. Shiga, T. Kodama, J. Kono, and J. Shiomi, One-directional thermal transport in densely aligned single-wall carbon nanotube films, Applied Physics Letters 115, 223104 (2019), https://doi.org/10.1063/1.5127209 .
  * Alemán  _et al._ [2015] B. Alemán, V. Reguero, B. Mas, and J. J. Vilatela, Strong carbon nanotube fibers by drawing inspiration from polymer fiber spinning, ACS Nano 9, 7392 (2015), pMID: 26082976, https://doi.org/10.1021/acsnano.5b02408 .
  * Issman  _et al._ [2022] L. Issman, P. A. Kloza, J. Terrones Portas, B. Collins, A. Pendashteh, M. Pick, J. J. Vilatela, J. A. Elliott, and A. Boies, Highly oriented direct-spun carbon nanotube textiles aligned by in situ radio-frequency fields, ACS Nano 16, 9583 (2022), pMID: 35638849, https://doi.org/10.1021/acsnano.2c02875 .
  * Bulmer  _et al._ [2021] J. Bulmer, A. Kaniyoor, and J. Elliott, A meta‐analysis of conductive and strong carbon nanotube materials, Advanced Materials 33, 2008432 (2021).
  * Lanir [1983] Y. Lanir, Constitutive equations for fibrous connective tissues, Journal of Biomechanics 16, 1 (1983).
  * Chen and Cheng [1996] C.-H. Chen and C.-H. Cheng, Effective elastic moduli of misoriented short-fiber composites, International Journal of Solids and Structures 33, 2519 (1996).
  * Lynch  _et al._ [2003] H. A. Lynch, W. Johannessen, J. P. Wu, A. Jawa, and D. M. Elliott, Effect of Fiber Orientation and Strain Rate on the Nonlinear Uniaxial Tensile Material Properties of Tendon , Journal of Biomechanical Engineering 125, 726 (2003), https://doi.org/10.1115/1.1614819 .
  * Lake  _et al._ [2009] S. P. Lake, K. S. Miller, D. M. Elliott, and L. J. Soslowsky, Effect of fiber distribution and realignment on the nonlinear and inhomogeneous mechanical properties of human supraspinatus tendon under longitudinal tensile loading, Journal of Orthopaedic Research 27, 1596 (2009), https://doi.org/10.1002/jor.20938 .
  * Martin and Ishida [1989] R. Martin and J. Ishida, The relative effects of collagen fiber orientation, porosity, density, and mineralization on bone strength, Journal of Biomechanics 22, 419 (1989).
  * Conklin  _et al._ [2011] M. W. Conklin, J. C. Eickhoff, K. M. Riching, C. A. Pehlke, K. W. Eliceiri, P. P. Provenzano, A. Friedl, and P. J. Keely, Aligned collagen is a prognostic signature for survival in human breast carcinoma, The American Journal of Pathology 178, 1221 (2011).
  * Riching  _et al._ [2014] K. Riching, B. L. Cox, M. Salick, C. Pehlke, A. Riching, S. M. Ponik, B. Bass, W. Crone, Y. Jiang, A. M. Weaver, K. Eliceiri, and P. Keely, 3d collagen alignment limits protrusions to enhance breast cancer cell persistence, Biophysical Journal 107, 2546 (2014).
  * Drifka  _et al._ [2016] C. R. Drifka, A. G. Loeffler, K. Mathewson, A. Keikhosravi, J. C. Eickhoff, Y. Liu, S. M. Weber, W. J. Kao, and K. W. Eliceiri, Highly aligned stromal collagen is a negative prognostic factor following pancreatic ductal adenocarcinoma resection, Oncotarget 7, 76197 (2016).
  * Bataller and Brenner [2005] R. Bataller and D. A. Brenner, Liver fibrosis, J Clin Invest 115, 209 (2005).
  * van Gurp [1995] M. van Gurp, The use of rotation matrices in the mathematical description of molecular orientations in polymers, Colloid and Polymer Science 273, 607 (1995).
  * Hermans  _et al._ [1946] J. J. Hermans, P. H. Hermans, D. Vermaas, and A. Weidinger, Quantitative evaluation of orientation in cellulose fibres from the x-ray fibre diagram, Recueil des Travaux Chimiques des Pays-Bas 65, 427 (1946), https://doi.org/10.1002/recl.19460650605 .
  * Vainio  _et al._ [2014] U. Vainio, T. I. W. Schnoor, S. Koyiloth Vayalil, K. Schulte, M. Müller, and E. T. Lilleodden, Orientation distribution of vertically aligned multiwalled carbon nanotubes, The Journal of Physical Chemistry C 118, 9507 (2014), https://doi.org/10.1021/jp501060s .
  * Goh  _et al._ [2005] K. L. Goh, J. Hiller, J. L. Haston, D. F. Holmes, K. E. Kadler, A. Murdoch, J. R. Meakin, and T. J. Wess, Analysis of collagen fibril diameter distribution in connective tissues using small-angle x-ray scattering, Biochimica et Biophysica Acta (BBA) - General Subjects 1722, 183 (2005).
  * Crawshaw  _et al._ [2002] J. Crawshaw, W. Bras, G. R. Mant, and R. E. Cameron, Simultaneous saxs and waxs investigations of changes in native cellulose fiber microstructure on swelling in aqueous sodium hydroxide, Journal of Applied Polymer Science 83, 1209 (2002), https://onlinelibrary.wiley.com/doi/pdf/10.1002/app.2287 .
  * Kratky [1933] O. Kratky, Zum deformationsmechanismus der faserstoffe, i, Kolloid-Zeitschrift 64, 213 (1933).
  * Tan  _et al._ [2006] J. C. Tan, J. A. Elliott, and T. W. Clyne, Analysis of tomography images of bonded fibre networks to measure distributions of fibre segment length and fibre orientation, Advanced Engineering Materials 8, 495 (2006), https://doi.org/10.1002/adem.200600033 .
  * Shaffer  _et al._ [1998] M. Shaffer, X. Fan, and A. Windle, Dispersion and packing of carbon nanotubes, Carbon 36, 1603 (1998).
  * Marquez [2006] J. P. Marquez, Fourier analysis and automated measurement of cell and fiber angular orientation distributions, International Journal of Solids and Structures 43, 6413 (2006).
  * Usov and Mezzenga [2015] I. Usov and R. Mezzenga, Fiberapp: An open-source software for tracking and analyzing polymers, filaments, biomacromolecules, and fibrous objects, Macromolecules 48, 1269 (2015), https://doi.org/10.1021/ma502264c .
  * Kaniyoor  _et al._ [2021] A. Kaniyoor, T. S. Gspann, J. E. Mizen, and J. A. Elliott, Quantifying alignment in carbon nanotube yarns and similar two-dimensional anisotropic systems, Journal of Applied Polymer Science 138, 50939 (2021), https://doi.org/10.1002/app.50939 .
  * Zhou and Zheng [2008] Y. Zhou and Y.-P. Zheng, Estimation of muscle fiber orientation in ultrasound images using revoting hough transform (rvht), Ultrasound in medicine & biology 34, 1474 (2008).
  * Pourdeyhimi and Kim [2002] B. Pourdeyhimi and H. Kim, Measuring fiber orientation in nonwovens: The hough transform, Textile Research Journal 72, 803 (2002).
  * Bredfeldt  _et al._ [2014] J. S. Bredfeldt, Y. Liu, C. A. Pehlke, M. W. Conklin, J. M. Szulczewski, D. R. Inman, P. J. Keely, R. D. Nowak, T. R. Mackie, and K. W. Eliceiri, Computational segmentation of collagen fibers from second-harmonic generation images of breast cancer, Journal of Biomedical Optics 19, 1 (2014).
  * Dong  _et al._ [2022] L. Dong, H. Hang, J. G. Park, W. Mio, and R. Liang, Detecting carbon nanotube orientation with topological analysis of scanning electron micrographs, Nanomaterials 12, 1251 (2022).
  * Behabtu  _et al._ [2013] N. Behabtu, C. C. Young, D. E. Tsentalovich, O. Kleinerman, X. Wang, A. W. Ma, E. A. Bengio, R. F. ter Waarbeek, J. J. de Jong, R. E. Hoogerwerf, _et al._ , Strong, light, multifunctional fibers of carbon nanotubes with ultrahigh conductivity, science 339, 182 (2013).
  * Bates and Luckhurst [2003] M. A. Bates and G. R. Luckhurst, X-ray scattering patterns of model liquid crystals from computer simulation: Calculation and analysis, The Journal of Chemical Physics 118, 6605 (2003), https://doi.org/10.1063/1.1557525 .
  * Zamora-Ledezma  _et al._ [2008] C. Zamora-Ledezma, C. Blanc, M. Maugey, C. Zakri, P. Poulin, and E. Anglaret, Anisotropic thin films of single-wall carbon nanotubes from aligned lyotropic nematic suspensions, Nano Letters 8, 4103 (2008), pMID: 18956923, https://doi.org/10.1021/nl801525x .
  * Leadbetter and Norris [1979] A. Leadbetter and E. Norris, Distribution functions in three liquid crystals from x-ray diffraction measurements, Molecular Physics 38, 669 (1979), https://doi.org/10.1080/00268977900101961 .
  * Deutsch [1991] M. Deutsch, Orientational order determination in liquid crystals by x-ray diffraction, Phys. Rev. A 44, 8264 (1991).
  * Boudaoud  _et al._ [2014] A. Boudaoud, A. Burian, D. Borowska-Wykret, M. Uyttewaal, R. Wrzalik, D. Kwiatkowska, and O. Hamant, Fibriltool, an imagej plug-in to quantify fibrillar structures in raw microscopy images, Nature protocols 9, 457 (2014).
  * Ruland and Burger [2006] W. Ruland and C. Burger, Evaluation of equatorial orientation distributions, Journal of Applied Crystallography 39, 889 (2006), https://doi.org/10.1107/S0021889806038957 .
  * Agra-Kooijman  _et al._ [2018] D. M. Agra-Kooijman, M. R. Fisch, and S. Kumar, The integrals determining orientational order in liquid crystals by x-ray diffraction revisited, Liquid Crystals 45, 680 (2018), https://doi.org/10.1080/02678292.2017.1372526 .
  * Bracewell and Bracewell [2000] R. Bracewell and R. Bracewell,  _The Fourier Transform and Its Applications_, Electrical engineering series (McGraw Hill, 2000).
  * Abramowitz  _et al._ [1988] M. Abramowitz, I. A. Stegun, and R. H. Romer,  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ (American Association of Physics Teachers, 1988).
  * Zwillinger [2018] D. Zwillinger,  _CRC Standard Mathematical Tables and Formulas_, Advances in Applied Mathematics (CRC Press, 2018).
  * [44] See Supplemental Materials at [URL] for a detailed derivation of the inverse projection matrix.


#  Insights on Stochastic Dynamics for Transmission of Monkeypox: Biological
and Probabilistic Behaviour111

Ghaus ur Rahman<EMAIL_ADDRESS>Olena Tymoshenko<EMAIL_ADDRESS>Giulia Di Nunno<EMAIL_ADDRESS>“Department of Mathematics and
Statistics, University of Swat, Pakistan”  “Department of Mathematics,
University of Oslo, Norway”  “Department of Mathematical Analysis and
Probability Theory, NTUU Igor Sikorsky Kyiv Polytechnic Institute, Ukraine”

###### Abstract

The transmission of monkeypox is studied using a stochastic model taking into
account the biological aspects, the contact mechanisms and the demographic
factors together with the intrinsic uncertainties. Our results provide insight
into the interaction between stochasticity and biological elements in the
dynamics of monkeypox transmission. The rigorous mathematical analysis
determines threshold parameters for disease persistence. For the proposed
model, the existence of a unique global almost sure non-negative solution is
proven. Conditions leading to disease extinction are established. Asymptotic
properties of the model are investigated such as the speed of transmission.

###### keywords:

Monkeypox Virus; Stochastic Model; Extinction; Long-Time Behaviour

##  1 Introduction

The modeling of disease transmission heavily relies on stochastic differential
equations (SDEs), which offer a more realistic framework for comprehending the
inherent uncertainty involved in the innate variety within a population. The
authors of [1] provided a comprehensive framework for understanding the
modeling of infectious diseases, blending both theoretical and practical
aspects of epidemiology, while [2] offers a deep dive into the principles of
epidemiological modeling and control strategies. Ecological models have been
studied using time series analysis [3] and tools of stochastic calculus [4, 5,
6], also see the references therein.

Monkeypox is an extremely contagious zoonotic disease that results from
infection by the monkeypox virus, a member of the Orthopoxvirus genus within
the Poxviridae family [7, 8]. Its impact is largely felt in rural regions
situated close to tropical rainforests in central and West Africa. However,
there is also an increasing number of reported cases in UK and in Europe.
Rodents are responsible for the spread and transmission of the zoonotic viral
disease monkeypox. Humans can either get monkeypox directly or indirectly
through rodents (mostly rats) and primates wildlife is major source of this
disease.

Furthermore, it is commonly transmitted between individuals through various
modes, such as respiratory droplets, contact with contaminated bodily fluids
or objects, as well as through the skin lesions of infected individuals [9].
After the successful eradication of smallpox, the incidence of the orthopox
virus known as monkeypox has surged, making it the most commonly encountered
variant [10]. The comprehension of the ecological dynamics of rodents is
important for controlling the propagation of the virus by the decreasing
rodent-human interaction [11, 12]. Among various methods in modeling, the
compartmental approach provides an invaluable framework for describing and
forecasting the dynamics of infectious diseases, like monkeypox. General
theory of the compartmental approach can be found in e.g., [13, 14, 15, 16],
and also references therein.

In the literature, we find some papers modelling monkeyfox transmission, with
deterministic perspective. As illustration, [17, 18, 19, 20] explore the
matter explaining the integer order formulation, as well as the their
extensions to the fractional order system by employing $\Psi$-Hilfer
fractional derivative and with Mittag-Leffler kernel.

We also use compartmental modeling. In comparison to the existing literature
on dynamics of monkeypox, our paper contributes to explore complex aspects of
the diseases including the uncertainty in the model. Indeed, an epidemic’s
behaviour is unpredictable due to external environmental disturbances that
impact its propagation. We propose a compartmental model where we have two
populations. The human population is composed of four distinct categories:
susceptible individuals $S_{h}$, infected individuals $I_{h}$, isolated
individuals $Q_{h}$, and recovered individuals $R_{h}$. On the other hand, the
rodent population can be divided into two groups: susceptible rodents $S_{r}$
and infected rodents $I_{r}$.

Figure 1: The schematic diagram of Monkey Pox virus transition in the
deterministic case

Starting from the deterministic model (see Figure 1), to be able to account
for the uncertainty, we add stochastic noises modelled by Brownian motions.
These act as a perturbation of the system mostly due to the effects of the
uncertainty of death and birth rates as well as the imprecise count of
individuals. Indeed, we consider the following stochastic compartment model
for monkeypox: $(S_{h}~{},I_{h},~{}Q_{h},~{}R_{h},)$ related to human
population and $(S_{r},~{}I_{r})$ for rodents:

$\displaystyle dS_{h}(t)$ $\displaystyle=$
$\displaystyle\left(\theta_{h}-(1-p)\cdot\frac{\eta_{1}I_{r}(t)+\eta_{2}I_{h}(t)}{N_{h}(t)}\cdot
S_{h}(t)-\mu_{h}S_{h}(t)\right)dt$ (3)
$\displaystyle-(1-p)\cdot\frac{\sigma_{1}I_{r}(t)dB_{1}(t)+\sigma_{2}I_{h}(t)dB_{2}(t)}{N_{h}(t)}\cdot
S_{h}(t)+\sigma_{3}S_{h}(t)dB_{3}(t);$ $\displaystyle dI_{h}(t)$
$\displaystyle=$
$\displaystyle\left((1-p)\cdot\frac{\eta_{1}I_{r}(t)+\eta_{2}I_{h}(t)}{N_{h}(t)}\cdot
S_{h}(t)-(\mu_{h}+\delta_{h}+\zeta)I_{h}(t)\right)dt$ (7)
$\displaystyle+(1-p)\cdot\frac{\sigma_{2}S_{h}(t)}{N_{h}(t)}\cdot
I_{h}(t)dB_{2}(t)+\sigma_{4}I_{h}(t)dB_{4}(t);$ $\displaystyle dQ_{h}(t)$
$\displaystyle=$ $\displaystyle\left(\zeta
I_{h}(t)-(\mu_{h}+\gamma_{h}+(1-\theta)\delta_{h})Q_{h}(t)\right)dt+\sigma_{5}Q_{h}(t)dB_{5}(t);$
(9) $\displaystyle dR_{h}(t)$ $\displaystyle=$
$\displaystyle\left(\gamma_{h}Q_{h}(t)-\mu_{h}R_{h}(t)\right)dt+\sigma_{6}R_{h}(t)dB_{6}(t);$
(11) $\displaystyle dS_{r}(t)$ $\displaystyle=$
$\displaystyle\left(\theta_{r}-\frac{\eta_{3}S_{r}(t)I_{r}(t)}{N_{r}(t)}-\mu_{r}S_{r}(t)\right)dt+\sigma_{7}S_{r}(t)dB_{7}(t);$
(13) $\displaystyle dI_{r}(t)$ $\displaystyle=$
$\displaystyle\left(\frac{\eta_{3}S_{r}(t)I_{r}(t)}{N_{r}(t)}-(\mu_{r}+\delta_{r})I_{r}(t)\right)dt+\sigma_{8}I_{r}(t)dB_{8}(t),$
(15)

where
$(S_{h}(0),I_{h}(0),Q_{h}(0),R_{h}(0),S_{r}(0)),I_{r}(0)\in\mathbb{R}^{6}_{+}$
is initial value, $B_{i}(t)$, $i=1,...,8,$ are mutually independent standard
Brownian motions defined on a filtered probabilistic space
$(\Omega,\mathbb{F},\\{{\cal F}_{t}\\}_{t\geq 0},\mathrm{P})$ such that
${\cal{F}}_{0}$ contains the null sets, and $\sigma_{i},i=1,2,...,8,$ are
their intensities (or volatilises). We also denote the involved human
population by $N_{h}:=S_{h}+I_{h}+Q_{h}+R_{h}$, and the rodent population by
$N_{r}:=S_{r}+I_{r}.$ Note that in this paper, we only consider the case when

$N_{r}(t)\leq k(t)N_{h}(t),$ (16)

where $k(t),~{}t\geq 0,$ is some continuous positive function with
$\bar{k}=\underset{t\geq 0}{\sup}\,k(t).$ Note in particular that the rate of
contacts, $\eta_{1}$ and $\eta_{2}$, are the basic parameters that are
responsible for monkeypox disease’s transmission. The susceptible human class
$S_{h}$ is generated by a constant recruitment rate $\theta_{h}$ and reduced
by acquiring infection after interaction with infectious humans and rodents.
The force of infection is $\eta_{1}I_{r}+\eta_{2}I_{h}$, where $\eta_{1}$ and
$\eta_{2}$ are the effective contact rates (capable of transmitting infection)
corresponding to infectious rodents and humans, respectively. The mortality
rate of each class within the human population decreases over time due to
natural deaths occurring at a rate of $\mu_{h}$. Meanwhile, the death rate
caused by monkeypox in the infected and isolated human compartments is
indicated by $\delta_{h}$. The infectious human population is isolated and
joins the $Q_{h}$ class at the transmission rate $\zeta$. The recovery rate of
the isolated/quarantine human population is denoted by $\gamma_{h}$. The value
$p$, with $0\leq p\leq 1$, measures the effectiveness of enlightenment
campaign, and $\theta$, with $0\leq\theta\leq 1$ is the effectiveness of
quarantine and treatment. It is assumed that the death in $Q_{h}$ due to
disease is influenced by the effectiveness of treatment, hence it is
$(1-\theta)\delta_{h}$ .

In the rodent population, $\theta_{r}$ shows the recruitment rate while the
natural death rate is $\mu_{r}$. The infection rate of the virus from an
infected rodent to another rodent is denoted as $\eta_{3}I_{r}$, with
$\eta_{3}$ representing the transmission rate. The population of each class of
rodents decreases naturally due to a death rate denoted as $\mu_{r}$, while
the death rate caused by monkeypox infection in the infected rodent group is
denoted as $\delta_{r}$. Since wild rodents do not generically have access to
treatment, we assume that they do not recover from the disease. A summary of
the model elements is presented in Table 1.

Model’s Elements | Description of Elements
---|---
$N_{h}$ | Humans population
$S_{h}$ | Humans who are susceptible
$I_{h}$ | Humans who are infected
$Q_{h}$ | Infected Humans who are quarantine
$R_{h}$ | Humans who are Recovered
$N_{r}$ | Rodents population
$S_{r}$ | Rodents who are susceptible
$I_{r}$ | Rodents who are infected
$\theta_{h}$ | Recruitment Rate of Humans
$p$ | Effectiveness Public Enlightenment Campaign
$\eta_{1}$ | Rate of Contact Between Rodents and Humans
$\eta_{2}$ | Rate of Contact Between Humans and Humans
$\eta_{3}$ | Rate of Contact Between Rodents and Rodents
$\mu_{h}$ | The rate at which humans naturally pass away
$\delta_{h}$ | Mortality rate caused by diseases in humans
$\zeta$ | Progression Rate from Infected to Quarantine
$\gamma_{h}$ | Recovery Rate of Humans
$\theta$ | Effectiveness of Quarantine and Treatment
$\theta_{r}$ | Rodents growth rate from birth to adult stage
$\mu_{r}$ | The rate at which rodents naturally pass away
$\delta_{r}$ | Mortality rate caused by diseases in rodents
Table 1: Descriptions of the elements

One described the stochastic model for monkeypox transmission, we proceed with
a complete analysis of the existence and positivity of the solution to the
stochastic dynamic system (3) in Section 2. We study stochastic ultimate
boundedness and permanence of the solution in Section 3. Conditions for the
disease extinction are detailed in Section 4, where we also provide
illustration of the obtained results with numerical examples realized with
Python. Section 5 presents an extended model in which the parameters are time
varying (non-autonomous system). With this, we study the long-time asymptotic
behaviour of the model and give results on the speed of the transmission.

##  2 Global Existence and Positivity of the Solution

In studying dynamic behavior, it is crucial to establish whether a solution
exists globally, i.e., there is no explosion in a finite time. In addition, in
the context of epidemic population dynamics, it is also important to assess
whether the solutions are non-negative.

A necessary condition for a stochastic differential equation to have a unique
global solution for any given initial value, is that the coefficients have the
linear growth, and are locally Lipschitz. (e.g. see [21], [22], [23]).
However, although the coefficients of system (3) do not satisfy linear growth
condition due to the bilinear incidence, they are locally Lipschitz
continuous. Nonetheless, the following result demonstrates that the solution
of system (3) is exists globally, it is unique, and non-negative.

######  Theorem 2.1.

The system represented by equation (3) has a unique solution,

$\mathcal{X}(t)=(S_{h}(t),I_{h}(t),Q_{h}(t),R_{h}(t),S_{r}(t),I_{r}(t)),$

on $t\geq 0$ given an initial condition of
$\mathcal{X}(0)\in\mathbb{R}^{6}_{+}$, and this solution remains in the space
$\mathbb{R}^{6}_{+}$ with a probability of 1, namely,
$\mathcal{X}(t)\in\mathbb{R}^{6}_{+}$ for all $t\geq 0$ almost surely (a.s.).

###### Proof.

Since the coefficients of the equation are locally Lipschitz continuous for
any given initial value $\mathcal{X}(0)\in\mathbb{R}^{6}_{+}$, by the theory
of stochastic differential equations there is a unique local solution
$\mathcal{X}(t)$ on $t\in[0,\tau_{e})$, where $\tau_{e}$ is the explosion time
(see [22]). To prove that this solution is global, we have to show that that
$\tau_{e}=\infty$ a.s.

This can be done by setting a sufficiently large value for $m_{0}\geq 0$, such
that $\mathcal{X}(0)\in\left[\frac{1}{m_{0}},m_{0}\right]^{\times 6}$. Then we
define a stopping time for each integer $m\geq m_{0}$:

$\displaystyle\tau_{m}=$ $\displaystyle\inf$
$\displaystyle\left\\{t\in[0,\tau_{e}):S_{h}(t)\notin\left(\frac{1}{m},m\right)\lor
I_{h}(t)\notin\left(\frac{1}{m},m\right)\lor
Q_{h}(t)\notin\left(\frac{1}{m},m\right)\lor\right.$
$\displaystyle\left.R_{h}(t)\notin\left(\frac{1}{m},m\right)\lor
S_{r}(t)\notin\left(\frac{1}{m},m\right)\lor
I_{r}(t)\notin\left(\frac{1}{m},m\right)\right\\}.$

By definition, $\tau_{m}$ is increasing as $m\to\infty$, and we set
$\tau_{\infty}=\underset{m\to\infty}{\lim}\tau_{m}$, then

$\tau_{\infty}\leq\tau_{e}\,\,\textnormal{ a.s.}$

To complete the proof, we have to show that $\tau_{\infty}=\infty$ a.s. If
this statement is false, then we can find a pair of constants $T>0$ and
$\varepsilon\in(0,1)$ such that $P\\{\tau_{\infty}\leq T\\}>\varepsilon$.
Hence, there is an integer $m_{1}\geq m_{0}$ such that

$\displaystyle\mathrm{P}\\{\tau_{m}\leq T\\}\geq\varepsilon\,\,\textnormal{for
all}\,\,\,m\geq m_{1}.$ (17)

With this we consider the following function
$\mathcal{V}:\mathbb{R}^{6}_{+}\to\mathbb{R}_{+}$:

$\displaystyle\mathcal{V}(S_{h},I_{h},Q_{h},R_{h},S_{r},I_{r})$
$\displaystyle=$ $\displaystyle(S_{h}-1-\log S_{h})+(I_{h}-1-\log I_{h})$
$\displaystyle+(Q_{h}-1-\log Q_{h})+(R_{h}-1-\log R_{h})$
$\displaystyle+(S_{r}-1-\log S_{r})+(I_{r}-1-\log I_{r}).$

This function is non-negative, since $x-1-\log x\geq 0$ for all $x>0$.

By Itô’s formula, it holds that

$\displaystyle d\mathcal{V}(\mathcal{X}(t))=L\mathcal{V}(\mathcal{X}(t))dt$
(18) (19)
$\displaystyle+\frac{(1-p)\cdot(\sigma_{1}I_{r}(t)S_{h}(t)dB_{1}(t)+\sigma_{2}I_{h}(t)S_{h}(t)dB_{2}(t))}{N_{h}(t)}\cdot\left(\frac{1}{S_{h}(t)}-\frac{1}{I_{h}(t)}\right)$
(20) (21)
$\displaystyle+(S_{h}(t)-1)\sigma_{3}dB_{3}(t)+(I_{h}(t)-1)\sigma_{4}dB_{4}(t)+(Q_{h}(t)-1)\sigma_{5}dB_{5}(t)$
(22) (23)
$\displaystyle+(R_{h}(t)-1)\sigma_{6}dB_{6}(t)+(S_{r}(t)-1)\sigma_{7}dB_{7}(t)+(I_{r}(t)-1)\sigma_{8}dB_{8}(t),$
(24)

where $L\mathcal{V}:\mathbb{R}^{6}_{+}\to\mathbb{R}_{+}$ is defined by

$\displaystyle
L\mathcal{V}(\mathcal{X}(t))=\theta_{h}+\theta_{r}+4\mu_{h}+\delta_{h}+\zeta+2\mu_{r}+\delta_{r}+\gamma_{h}+(1-\theta)\delta_{h}$
$\displaystyle+(1-p)\cdot\frac{(\eta_{1}I_{r}(t)+\eta_{2}I_{h}(t))}{N_{h}(t)}+\frac{\eta_{3}I_{r}(t)}{N_{r}(t)}-\mu_{h}S_{h}(t)-(\mu_{h}+\delta_{h})I_{h}(t)$
$\displaystyle-(\mu_{h}+(1-\theta)\delta_{h})Q_{h}(t)-\mu_{h}R_{h}(t)-\mu_{r}S_{r}(t)-(\mu_{r}+\delta_{r})I_{r}(t)-\frac{\theta_{h}}{S_{h}(t)}$
$\displaystyle-(1-p)\cdot\frac{(\eta_{1}I_{r}(t)+\eta_{2}I_{h}(t))S_{h}(t)}{N_{h}(t)I_{h}(t)}-\frac{\zeta
I_{h}(t)}{Q_{h}(t)}-\frac{\gamma_{h}Q_{h}(t)}{R_{h}(t)}-\frac{\theta_{r}}{S_{r}(t)}-\frac{\eta_{3}S_{r}(t)}{N_{r}(t)}$
$\displaystyle+(1-p)^{2}\cdot\frac{\sigma^{2}_{1}I^{2}_{r}(t)+\sigma^{2}_{2}I^{2}_{h}(t)}{2N^{2}_{h}(t)}+(1-p)^{2}\cdot\frac{\sigma^{2}_{2}S^{2}_{h}(t)}{2N^{2}_{h}(t)}$
$\displaystyle+\frac{\sigma^{2}_{3}}{2}+\frac{\sigma^{2}_{4}}{2}+\frac{\sigma^{2}_{5}}{2}+\frac{\sigma^{2}_{6}}{2}+\frac{\sigma^{2}_{7}}{2}+\frac{\sigma^{2}_{8}}{2}.$

Following assumption (16) we have that $\dfrac{I_{r}}{N_{h}}\leq\bar{k}$ and
$S_{h}<N_{h}$, $I_{r}<N_{r}$, $S_{r}<N_{r}$, $I_{h}<N_{h}$, then we can see
that

$\displaystyle L\mathcal{V}(\mathcal{X}(t))$ $\displaystyle\leq K,$ (25)

with the constant $K$ given by

$\displaystyle K$ $\displaystyle=$
$\displaystyle\theta_{h}+\theta_{r}+4\mu_{h}+\delta_{h}+\zeta+2\mu_{r}+\delta_{r}+\gamma_{h}+(1-\theta)\delta_{h}+(1-p)\cdot(\eta_{1}\bar{k}+\eta_{2})$
$\displaystyle+\eta_{3}+(1-p)^{2}\cdot\frac{\sigma^{2}_{1}}{2}\cdot{\bar{k}}^{2}+(1-p)^{2}\sigma^{2}_{2}+\frac{\sigma^{2}_{3}}{2}+\frac{\sigma^{2}_{4}}{2}+\frac{\sigma^{2}_{5}}{2}+\frac{\sigma^{2}_{6}}{2}+\frac{\sigma^{2}_{7}}{2}+\frac{\sigma^{2}_{8}}{2}.$

Integrating both sides of equation (18) from 0 to $\tau_{m}\wedge
t=\min\\{\tau_{m},t\\}$, where $t\in[0,T]$, and by (25) we obtain the
following relation

$\displaystyle\int^{\tau_{m}\wedge
t}_{0}d\mathcal{V}\bigg{(}\mathcal{X}(r)\bigg{)}$ $\displaystyle\leq$
$\displaystyle\int^{\tau_{m}\wedge t}_{0}Kdr-(1-p)\int^{\tau_{m}\wedge
t}_{0}\frac{(S_{h}(r)-1)\sigma_{1}I_{r}(r)}{N_{h}(r)}dB_{1}(r)$
$\displaystyle+(1-p)\int^{\tau_{m}\wedge
t}_{0}\frac{(I_{h}(r)-S_{h}(r))\sigma_{2}}{N_{h}(r)}dB_{2}(r)$
$\displaystyle+\int^{\tau_{m}\wedge
t}_{0}(S_{h}(r)-1)\sigma_{3}dB_{3}(r)+\int^{\tau_{m}\wedge
t}_{0}(I_{h}(r)-1)\sigma_{4}dB_{4}(r)$ $\displaystyle+\int^{\tau_{m}\wedge
t}_{0}(Q_{h}(r)-1)\sigma_{5}dB_{5}(r)+\int^{\tau_{m}\wedge
t}_{0}(R_{h}(r)-1)\sigma_{6}dB_{6}(r)$ $\displaystyle+\int^{\tau_{m}\wedge
t}_{0}(S_{r}(r)-1)\sigma_{7}dB_{7}(r)+\int^{\tau_{m}\wedge
t}_{0}(I_{r}(r)-1)\sigma_{8}dB_{8}(r),$

since $\mathcal{X}(\tau_{m}\wedge t)\in\mathbb{R}^{6}_{+}$, $t\in[0,T]$.

The rest of the proof follows the arguments in Mao et. al. [25]. Taking the
expectation on both sides above, we have

$\displaystyle\mathrm{E}\mathcal{V}\left(\mathcal{X}(\tau_{m}\wedge T)\right)$
$\displaystyle\leq$
$\displaystyle\mathcal{V}(\mathcal{X}(0))+\mathrm{E}\int^{\tau_{m}\wedge
T}_{0}Kdr$ (28)
$\displaystyle=\mathcal{V}(\mathcal{X}(0))+K\mathrm{E}(\tau_{m}\wedge
T)\leq\mathcal{V}(\mathcal{X}(0))+KT.$

Set $A_{m}=\\{\tau_{m}\leq t\\}$ for $m\geq m_{1}$. By (17), we have
$P(A_{m})\geq\varepsilon$. Note that for every $\omega\in A_{m}$, there is at
least one term among these five $S_{h}(\tau_{m}\wedge T)$,
$I_{h}(\tau_{m}\wedge T)$, $Q_{h}(\tau_{m}\wedge T)$, $R_{h}(\tau_{m}\wedge
T)$, $S_{r}(\tau_{m}\wedge T)$ and $I_{r}(\tau_{m}\wedge T)$ which is equal to
either $m$ or $\dfrac{1}{m}$. Then

$\mathcal{V}(\mathcal{X}(\tau_{m}\wedge T))\geq\left(\left(m-1-\log
m\right)\wedge\left(\frac{1}{m}-1-\log\frac{1}{m}\right)\right).$

Therefore, from (17) and (28), we can easily obtain

$\displaystyle\mathcal{V}(\mathcal{X}(0))+KT$ $\displaystyle\geq$
$\displaystyle\mathrm{E}\left(\mathbb{I}_{A_{m}(\omega)}\mathcal{V}(\mathcal{X}(\tau_{m}))\right)$
$\displaystyle\geq\varepsilon\left(\left(m-1-\log
m\right)\wedge\left(\frac{1}{m}-1-\log\frac{1}{m}\right)\right),$

where $\mathbb{I}_{A_{m}(\omega)}$ is the indicator function of $A_{m}$.
Letting $m$ increase to infinity, we find a contradiction:
$\infty\leq\mathcal{V}(\mathcal{X}(0))+KT=\infty.$ So it must be
$\tau_{\infty}=\infty$ a.s. The non-negativity of the solution is guaranteed
as a product of the argument above thanks to the definition of
$(\tau_{m})_{m\geq m_{0}}$ and $\tau_{\infty}=\infty$ a.s. ∎

The result above shows that the proposed model is well-posed and provides
plausible future scenarios for the evolution of the disease depending on the
initial circumstances (i.e., initial conditions), representing the population
conditions at the start of the epidemic, and the set of parameters described
in Table 1.

Particularly, it is important that the trajectories keep positive since a
negative solution would not have meaningful interpretation.

##  3 Stochastically Ultimately Bounded and Permanent

After establishing the well-posedness of the system (3), we proceed detailing
other properties. Now, we present the definition of stochastic ultimate
boundedness and permanent solution of the system, see [24].

######  Definition 3.1.

The solutions $\mathcal{X}(t)$, $t\geq 0$, of model (3) are said to be
stochastically ultimately bounded, if for any $\epsilon\in(0,1)$, there is a
positive constant $\kappa=\kappa(\epsilon)$, such that for any initial value
$\mathcal{X}(0)\in\mathbb{R}^{6}_{+}$, the solution to (3) satisfies

$\underset{t\to\infty}{\lim\sup}\,\mathrm{P}\left\\{|\mathcal{X}(t)|>\kappa\right\\}<\epsilon.$

######  Definition 3.2.

The solutions $\mathcal{X}(t)$, $t\geq 0$ of model (3) are said to be
stochastically permanent if for any $\epsilon\in(0,1)$, there exists a pair of
positive constants $\kappa=\kappa(\epsilon)$ and $\chi=\chi(\epsilon)$, such
that for any initial value $\mathcal{X}(0)\in\mathbb{R}^{6}_{+}$, the solution
to (3) has the following properties

$\underset{t\to\infty}{\lim\inf}\,\mathrm{P}\left\\{|\mathcal{X}(t)|\geq\kappa\right\\}\geq
1-\epsilon;$
$\underset{t\to\infty}{\lim\inf}\,\mathrm{P}\left\\{|\mathcal{X}(t)|\leq\chi\right\\}\geq
1-\epsilon.$

######  Proposition 3.3.

For any initial value $\mathcal{X}(0)\in\mathbb{R}^{6}_{+}$, let model (3)
satisfy

$\mathcal{K}(t)=N_{h}(t)+N_{r}(t)>0\,\,\,\text{for all }\,\,\,t\geq 0.$ (29)

Then the solution $\mathcal{X}(t)$, $t\geq 0$ is stochastically ultimately
bounded and permanent.

From a biological point of view, assumption (29) is the relevant case, since
$\underset{t\geq 0}{\inf}\,\mathcal{K}(t)=0$ would also mean that the
populations are extinct.

###### Proof.

Define $\mathcal{V}(\mathcal{X}(t))=\mathcal{K}(t)+\frac{1}{\mathcal{K}(t)}$.
Applying the Itô formula, we get

$\displaystyle d\mathcal{V}(\mathcal{X}(t))=L\mathcal{V}(\mathcal{X}(t))dt$
$\displaystyle+\left(1-\frac{1}{\mathcal{K}^{2}(t)}\right)\bigg{[}\frac{(p-1)\sigma_{1}I_{r}(t)S_{h}(t)}{N_{h}(t)}dB_{1}(t)+\sigma_{3}S_{h}(t)dB_{3}(t)+\sigma_{4}I_{h}(t)dB_{4}(t)$
$\displaystyle+\sigma_{5}Q_{h}(t)dB_{5}(t)+\sigma_{6}R_{h}(t)dB_{6}(t)+\sigma_{7}S_{r}(t)dB_{7}(t)+\sigma_{8}I_{r}(t)dB_{8}(t)\bigg{]},$

where

$\displaystyle
L\mathcal{V}(\mathcal{X}(t))=\left(1-\frac{1}{\mathcal{K}^{2}(t)}\right)\bigg{(}\theta_{h}-\mu_{h}(S_{h}(t)+I_{h}(t)+Q_{h}(t)+R_{h}(t))$
$\displaystyle-\delta_{h}I_{h}(t)-(1-\theta)\delta_{h}Q_{h}(t)+\theta_{r}-\mu_{r}(S_{r}(t)+I_{r}(t))-\delta_{r}I_{r}(t)\bigg{)}$
$\displaystyle+\frac{\sigma^{2}_{1}(1-p)^{2}I^{2}_{r}(t)S^{2}_{h}(t)}{N^{2}_{h}\mathcal{K}^{3}(t)}+\frac{2\sigma^{2}_{2}(1-p)^{2}I^{2}_{h}(t)S^{2}_{h}(t)}{N^{2}_{h}\mathcal{K}^{3}(t)}$
$\displaystyle+\frac{\sigma^{2}_{3}S^{2}_{h}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{4}I^{2}_{h}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{5}Q^{2}_{h}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{6}R^{2}_{h}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{7}S^{2}_{r}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{8}I^{2}_{r}(t)}{\mathcal{K}^{3}(t)}.$

We can rewrite the function $L\mathcal{V}(\mathcal{X}(t)),\,t\geq 0,$ as

$\displaystyle
L\mathcal{V}(\mathcal{X}(t))=\theta_{h}-\delta_{h}I_{h}(t)+\theta_{r}-(1-\theta)\delta_{h}Q_{h}(t)+\frac{(1-\theta)\delta_{h}Q_{h}(t)}{\mathcal{K}^{2}(t)}$
$\displaystyle-
H_{1}\left(\mathcal{K}(t)+\frac{1}{\mathcal{K}(t)}\right)+\frac{2H_{1}}{\mathcal{K}(t)}+\frac{\delta_{r}I_{r}(t)}{\mathcal{K}^{2}(t)}$
$\displaystyle+\frac{\sigma^{2}_{1}(1-p)^{2}I^{2}_{r}(t)S^{2}_{h}(t)}{N^{2}_{h}\mathcal{K}^{3}(t)}+\frac{2\sigma^{2}_{2}(1-p)^{2}I^{2}_{h}(t)S^{2}_{h}(t)}{N^{2}_{h}\mathcal{K}^{3}(t)}$
$\displaystyle+\frac{\sigma^{2}_{3}S^{2}_{h}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{4}I^{2}_{h}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{5}Q^{2}_{h}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{6}R^{2}_{h}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{7}S^{2}_{r}(t)}{\mathcal{K}^{3}(t)}+\frac{\sigma^{2}_{8}I^{2}_{r}(t)}{\mathcal{K}^{3}(t)},$

for any $t$, with $H_{1}=\max\\{\mu_{h},\mu_{r}\\}$. Since (29) holds, and
$\dfrac{S_{h}}{N_{h}}<1$, we find the upper estimate:

$\displaystyle
L\mathcal{V}(\mathcal{X}(t))\leq\mathcal{F}-H_{1}\left(\mathcal{K}(t)+\frac{1}{\mathcal{K}(t)}\right),$

which is

$L\mathcal{V}(\mathcal{X}(t))\leq\mathcal{F}-H_{1}\mathcal{V}(t),$ (30)

where
$\mathcal{F}=\theta_{h}+\theta_{r}+\dfrac{\delta_{h}+\delta_{r}+2H_{1}}{\underset{t\geq
0}{\inf}\,\mathcal{K}(t)}$.

Let $\mathcal{V}_{1}(\mathcal{X}(t))=e^{H_{1}t}\mathcal{V}(\mathcal{X}(t))$,
using Itô formula, by (30) we have

$\displaystyle d(e^{H_{1}t}\mathcal{V}(\mathcal{X}(t))$ $\displaystyle\leq$
$\displaystyle
H_{1}e^{H_{1}t}\mathcal{V}(\mathcal{X}(t)))dt+e^{H_{1}t}(\mathcal{F}-H_{1}\mathcal{V}(\mathcal{X}(t)))dt$
$\displaystyle+e^{H_{1}t}\left(1-\frac{1}{\mathcal{K}^{2}(t)}\right)\bigg{[}\frac{(p-1)\sigma_{1}I_{r}(t)S_{h}(t)}{N_{h}(t)}dB_{1}(t)$
$\displaystyle+\sigma_{3}S_{h}(t)dB_{3}(t)+\sigma_{4}I_{h}(t)dB_{4}(t)+\sigma_{5}Q_{h}(t)dB_{5}(t)$
$\displaystyle+\sigma_{6}R_{h}(t)dB_{6}(t)+\sigma_{7}S_{r}(t)dB_{7}(t)+\sigma_{8}I_{r}(t)dB_{8}(t)\bigg{]}$
$\displaystyle=e^{H_{1}t}\mathcal{F}dt+e^{H_{1}t}\left(1-\frac{1}{\mathcal{K}^{2}(t)}\right)\bigg{[}\frac{(p-1)\sigma_{1}I_{r}(t)S_{h}(t)}{N_{h}(t)}dB_{1}(t)$
$\displaystyle+\sigma_{3}S_{h}(t)dB_{3}(t)+\sigma_{4}I_{h}(t)dB_{4}(t)+\sigma_{5}Q_{h}(t)dB_{5}(t)$
$\displaystyle+\sigma_{6}R_{h}(t)dB_{6}(t)+\sigma_{7}S_{r}(t)dB_{7}(t)+\sigma_{8}I_{r}(t)dB_{8}(t)\bigg{]}.$

Then

$\displaystyle\mathrm{E}(\mathcal{V}(\mathcal{X}(t)))$ $\displaystyle=$
$\displaystyle\mathrm{E}\left(e^{-H_{1}t}\mathcal{V}(\mathcal{X}(0))+\frac{\mathcal{F}}{H_{1}}\right)\leq\frac{\mathcal{F}}{H_{1}}+\mathcal{V}(\mathcal{X}(0))=\alpha.$

Let’s fix $\epsilon>0$, then exist $\kappa>0$ such that
$\dfrac{\alpha}{\kappa}<\epsilon$. By Markov’s inequality, for a sufficiently
large constant $\kappa$ we have

$P\left\\{\mathcal{K}(t)+\frac{1}{\mathcal{K}(t)}>\kappa\right\\}\leq\frac{1}{\kappa}\mathrm{E}\left(\mathcal{K}(t)+\frac{1}{\mathcal{K}(t)}\right)\leq\frac{\alpha}{\kappa}<\epsilon.$

Then

$\lim_{t\to\infty}\sup\mathrm{P}\left\\{\mathcal{K}(t)+\frac{1}{\mathcal{K}(t)}>\kappa\right\\}\leq\epsilon.$

That is

$\lim_{t\to\infty}\sup\mathrm{P}\\{\mathcal{K}(t)>\kappa\\}<\epsilon.$ (31)

Observe that, for $x\in\mathbb{R}^{d}_{+}$, we have
$|x|\leq\sqrt{d}\sum_{i=1}^{d}x_{i}$. Then by (31), we get

$\lim_{t\to\infty}\sup\mathrm{P}\\{|\mathcal{X}(t)|>\kappa\sqrt{6}\\}\leq\lim_{t\to\infty}\sup\mathrm{P}\\{\mathcal{K}(t)>\kappa\\}\leq\epsilon.$

According to the Definition 3.1, it means that the solutions $\mathcal{X}(t)$,
$t\geq 0$, of model (3) are stochastically ultimately bounded. Similarly can
be easily obtained the stochastically permanent property for solutions
according to the Definitions 3.2 (see [24]). This completes the proof. ∎

Stochastically ultimately boundedness indicates that the model’s depiction of
population dynamics won’t show unbounded expansion or decline in the actual
situation of monkeypox disease transmission. The system’s variables, which
represent various population segments, have certain upper bounds that they
will not ultimately surpass. According to the concept of stochastic eventually
boundedness, there is a strong likelihood that the solutions of the stochastic
differential equations driving the monkeypox model will remain within certain
bounds throughout time.

##  4 Study of Extinction of the Disease

In this section, we determine the factors that could lead to the extinction of
the disease.

######  Theorem 4.1.

Let $\mathcal{X}(t)$, $t\geq 0$, be a solution of system (3) with initial
value $\mathcal{X}(0)\in\mathbb{R}^{6}_{+}$. Assume that

$\mathcal{R}_{0}=\frac{(1-p)(\eta_{1}+\eta_{2})+\eta_{3}}{\min\\{\mu_{h},\mu_{r}\\}+\min\\{\delta_{h},\delta_{r}\\}}<1,$

then the number of infected individuals of system (3) goes to extinction
almost surely, i.e.,

$\underset{t\to\infty}{\lim}(I_{h}(t)+Q_{h}(t)+I_{r}(t))=0\,\,\,\textnormal{a.s.}$
(32)

###### Proof.

The Itô’s formula applied to $\mathcal{U}=ln(I_{h}+Q_{h}+I_{r})$, where
$\mathcal{U}:\mathbb{R}^{3}_{+}\to\mathbb{R}_{+}$, leads to

$\displaystyle d\mathcal{U}(t)$ $\displaystyle=$ $\displaystyle
L\mathcal{U}(t)dt+\frac{(1-p)\sigma_{2}I_{h}(t)S_{h}(t)dB_{2}(t)}{N_{h}(t)(I_{h}(t)+Q_{h}(t)+I_{r}(t))}$
(35)
$\displaystyle+\frac{\sigma_{4}I_{h}(t)dB_{4}(t)+\sigma_{5}Q_{h}(t)dB_{5}(t)+\sigma_{8}I_{r}(t)dB_{8}(t)}{I_{h}(t)+Q_{h}(t)+I_{r}(t)},$

where

$\displaystyle
L\mathcal{U}(t)=\frac{(1-p)}{N_{h}(t)(I_{h}(t)+Q_{h}(t)+I_{r}(t))}\cdot\Big{(}\eta_{1}I_{r}(t)S_{h}(t)+\eta_{2}I_{h}(t)S_{h}(t)\Big{)}$
(36) (37)
$\displaystyle+\frac{\eta_{3}S_{r}(t)I_{r}(t)}{N_{r}(t)(I_{h}(t)+Q_{h}(t)+I_{r}(t))}$
(38) (39)
$\displaystyle-\frac{(\mu_{h}+\delta_{h})I_{h}(t)+(\mu_{h}+\gamma_{h}+(1-\theta)\delta_{h})Q_{h}(t)+(\mu_{r}+\delta_{r})I_{r}(t)}{I_{h}(t)+Q_{h}(t)+I_{r}(t)}$
(40) (41)
$\displaystyle-\frac{(1-p)^{2}S^{2}_{h}(t)\sigma^{2}_{2}I^{2}_{h}(t)}{2N^{2}_{h}(t)(I_{h}(t)+Q_{h}(t)+I_{r}(t))^{2}}-\frac{\sigma^{2}_{4}I^{2}_{h}(t)+\sigma^{2}_{5}Q^{2}_{h}(t)+\sigma^{2}_{8}I^{2}_{r}(t)}{2(I_{h}(t)+Q_{h}(t)+I_{r}(t))^{2}}.$
(42)

Since $S_{h}\leq N_{h}$, $Q_{h}\leq N_{h}$, from (36), we get

$\displaystyle L\mathcal{U}(t)$ $\displaystyle\leq$
$\displaystyle(1-p)(\eta_{1}+\eta_{2})+\eta_{3}-(\min\\{\mu_{h},\mu_{r}\\}+\min\\{\delta_{h},\delta_{r}\\})$

Integrating on both side equation (35), and substituting (36), we obtain

$\displaystyle\mathcal{U}(t)-\mathcal{U}(0)\leq(\min\\{\mu_{h},\mu_{r}\\}+\min\\{\delta_{h},\delta_{r}\\})(\mathcal{R}_{0}-1)\,t$
$\displaystyle+(1-p)\sigma_{2}B_{2}(t)+\sigma_{4}B_{4}(t)+\sigma_{5}B_{5}(t)+\sigma_{8}B_{8}(t).$

By the strong low of large number (see [23]), for $i=2,4,5,8,$ we have

$\underset{t\to\infty}{\lim}\frac{B_{i}(t)}{t}=0\,\,\,\textnormal{a.s.}$ (43)

Since $\mathcal{R}_{0}<1$ (by assumption), $L(\mathcal{U})(t)\leq 0$, and by
(43) we obtain

$\underset{t\to\infty}{\lim\sup}\frac{\mathcal{U}(t)}{t}\leq(\min\\{\mu_{h},\mu_{r}\\}+\min\\{\delta_{h},\delta_{r}\\})(\mathcal{R}_{0}-1)<0\,\,\,\textnormal{a.s.}$

leading to (32), which completes the proof. ∎

To conform the analytical results above, we provide some numerical
illustrations. First we demonstrate the case for $R_{0}<1$.

######  Example 4.2.

Consider the model (3) with parameters value:

  * 1.

$(S_{h}(0),Q_{h}(0),I_{h}(0),R_{h}(0),S_{r}(0),I_{r}(0))=(90,50,60,70,80,30);$

  * 2.

$(\sigma_{1},~{}\sigma_{1},~{}\sigma_{3},~{}\sigma_{4},~{}\sigma_{5},~{}\sigma_{6},~{}\sigma_{7},~{}\sigma_{8})=(0.05,0.04,0.01,0.05,0.04,0.01,0.05,0.04).$

We also consider the values from Table 2 that will result in producing,
$R_{0}<1$, where the conditions of Theorem 4.1 are satisfied.

Model’s Elements | Particular Values of Elements
---|---
$\theta$ | 0.043
$\theta_{h}$ | 10
$p$ | 0.041
$\eta_{1}$ | 0.009
$\eta_{2}$ | 0.002
$\eta_{3}$ | 0.0027
$\mu_{h}$ | 0.05
$\delta_{h}$ | 0.003
$\zeta$ | 0.5
$\gamma_{h}$ | 0.2
$\theta_{r}$ | 10
$\mu_{r}$ | 0.02
$\delta_{r}$ | 0.004
Table 2: Value of the elements

The graphs are obtained while using python software for

$R_{0}=\frac{(1-0.041)(0.009+0.002)+0.0027}{0.02+0.003}=0.567<1.$

The obtained curves show the dynamics of the monkeypox disease transmission
between humans and rodents population by comparing the deterministic and
stochastic versions of the model. Also the uncertainty is quantified through
the means and standard deviation of $S_{h}$, $I_{h}$, $Q_{h}$, $R_{h}$,
$S_{r}$ and $I_{r}$. Standard deviation curves represent the variability of
each compartment over multiple simulations, indicating the uncertainty
associated with each curve.

Figure 2: The curves represent the behaviour of susceptible(a), infected(b),
isolated(c), recovered(d) humans populations and of different rodents classes
(e)-(f) as time evolve. Figure 3: Paths of the different classes of humans
(a)-(d) and rodents (e)-(f) along with mean of each family of curves. Means
are shown with dark curves surrounded by family of related trajectories
Figure 4: Curves show the variability of each compartment humans (a)-(d) and
rodents (e)-(f) over multiple simulations, indicating the uncertainty
associated with each curve. Figure 5: The graphs show the distribution of
values for each compartment humans (a)-(d) and rodents (e)-(f): The histograms
for Susceptible Humans(a) and Susceptible Rodents (e) indicate a high number
of individuals in these compartments, suggesting a large pool of potential
hosts for the disease. The histograms for Infected Humans (b) and Infected
Rodents (f) show a smaller number of individuals, indicating the disease
prevalence. The histograms for Isolated Humans (c) and Recovered Humans (d)
indicate the effectiveness of isolation and recovery measures.

Next we discuss the case for $R_{0}>1$.

######  Example 4.3.

To illustrate this we consider the formulated model (3) with the same initial
values and volatility $\sigma_{i},i=1,...,8$ defined in Example 4.2. We also
consider the values from Table 3 that will results in producing,

$R_{0}=\frac{(1-0.041)(0.009+0.002)+0.027}{0.02+0.003}=1.6326>1.$

In particular, we remark that in Example 4.2 the value $\eta_{3}=0.0027$,
while here it is $\eta_{3}=0.027.$

Model’s Elements | Particular Values of Elements
---|---
$\theta$ | 0.043
$\theta_{h}$ | 10
$p$ | 0.041
$\eta_{1}$ | 0.009
$\eta_{2}$ | 0.002
$\eta_{3}$ | 0.027
$\mu_{h}$ | 0.05
$\delta_{h}$ | 0.003
$\zeta$ | 0.5
$\gamma_{h}$ | 0.2
$\theta_{r}$ | 10
$\mu_{r}$ | 0.02
$\delta_{r}$ | 0.004
Table 3: Descriptions values of the elements

All trajectories after long time (800 days) with $R_{0}>1$ have different
nature as compared to the case when, $R_{0}<1$. Indeed, for $R_{0}<1$, it is
shown in Example 4.2 that all trajectories reach the endemic equilibrium point
after 40 days.

Figure 6: For $R_{0}>1$: the curves represent the behavior of susceptible(a),
infected(b), isolated(c), recovered(d) humans populations and of different
rodents classes (e)-(f) as time evolve for 800 days. Figure 7:  For $R_{0}>1$:
the curves without noise and mean of cures with noise for human population of
susceptible(a), infected(b), isolated(c), recovered(d) humans populations and
of different rodents classes (e)-(f) as time evolve for 800 days.  Figure 8:
For $R_{0}>1$, paths of the different classes of humans (a)-(d) and rodents
(e)-(f) along with mean of each family of curves. Means of each family of
paths are shown with dark curves surrounded by family of related trajectories.
After a long time 800 days the infected population still fluctuate.

Figure 9: For $R_{0}>1$, curves show the variability of each compartment
humans (a)-(d) and rodents (e)-(f) over multiple simulations, indicating the
uncertainty associated with each curve.

##  5 Speed of growth of time-depending stochastic model for monkey-pox

Several generalized SIR (Susceptible, Infected, and Recovered) models have
been studed in terms of their long-term behaviour. For instance, the SIR model
with perturbed disease transmission coefficients of [26] has the rates of
convergence determined for $t\to\infty$. It is shown that the rate is not too
far from exponential in the sense that the convergence speed has a polynomial
form of any degree. In [27], [28], the long-time behaviour of densities of the
distributions of the stochastic SIR model are studied.

In this section, we focus on the influence of the various system parameters on
the asymptotic behaviour of the non-autonomous stochastic model, which
corresponds system (3). For this study, the parameters of the model are
allowed to depend on time, to incorporate the seasonal variation. Hence the
model is of the form

$\displaystyle dS_{h}(t)$
$\displaystyle=\left(\theta_{h}(t)-(1-p)\left(\frac{\eta_{1}(t)I_{r}(t)+\eta_{2}(t)I_{h}(t)}{N_{h}(t)}\right)S_{h}(t)-\mu_{h}(t)S_{h}(t)\right)dt$
$\displaystyle-(1-p)\left(\frac{\sigma_{1}(t)I_{r}(t)dB_{1}(t)+\sigma_{2}(t)I_{h}(t)dB_{2}(t)}{N_{h}(t)}\right)S_{h}(t)+\sigma_{3}(t)S_{h}(t)dB_{3}(t),$
$\displaystyle dI_{h}(t)$
$\displaystyle=\left((1-p)\left(\frac{\eta_{1}(t)I_{r}(t)+\eta_{2}(t)I_{h}(t)}{N_{h}(t)}\right)S_{h}(t)-(\mu_{h}(t)+\delta_{h}(t)+\zeta(t))I_{h}(t)\right)dt$
$\displaystyle+(1-p)\frac{\sigma_{2}(t)S_{h}(t)I_{h}(t)}{N_{h}(t)}dB_{2}(t)+\sigma_{4}(t)I_{h}(t)dB_{4}(t),$
$\displaystyle dQ_{h}(t)$
$\displaystyle=\left(\zeta(t)I_{h}(t)-(\mu_{h}(t)+\gamma_{h}(t)+(1-\theta(t))\delta_{h}(t))Q_{h}(t)\right)dt+\sigma_{5}(t)Q_{h}(t)dB_{5}(t),$
$\displaystyle dR_{h}(t)$
$\displaystyle=\left(\gamma_{h}(t)Q_{h}(t)-\mu_{h}(t)R_{h}(t)\right)dt+\sigma_{6}(t)R_{h}(t)dB_{6}(t)$
$\displaystyle dS_{r}(t)$
$\displaystyle=\left(\theta_{r}(t)-\frac{\eta_{3}(t)S_{r}(t)I_{r}(t)}{N_{r}(t)}-\mu_{r}(t)S_{r}(t)\right)dt+\sigma_{7}(t)S_{r}(t)dB_{7},$
$\displaystyle dI_{r}(t)$
$\displaystyle=\left(\frac{\eta_{3}(t)S_{r}(t)I_{r}(t)}{N_{r}(t)}-(\mu_{r}(t)+\delta_{r}(t))I_{r}(t)\right)dt+\sigma_{8}(t)I_{r}(t)dB_{8}(t).$
(44)

Our goal is to determine which parameters have a significant effect on the
growth rate of the solution of this non-autonomous system. We follow similar
arguments to [29, 30] to achieve the results.

We stress that determining the rate of growth is essential for establishing
the severity of the outbreak. In this context, we consider the model (3) with
time-dependent parameters leading to (5). This makes it possible to take into
account, e.g., seasonal changes. We emphasize that the rate of growth provides
important insights into the dynamic nature of the disease, indicating whether
it is accelerating or decelerating.

We can still denote the solution of the time-depending stochastic epidemic
model (5), as

$\mathcal{X}(t)=(S_{h}(t),I_{h}(t),Q_{h}(t),R_{h}(t),S_{r}(t),I_{r}(t)),$

given an initial condition of $\mathcal{X}(0)\in\mathbb{R}^{6}_{+}$ on $t\geq
0$. The following theorem establish the well-posedness of the system (5).

######  Theorem 5.1.

Let the coefficients of the stochastic system (5) be continuous functions with
upper bounds. Then for every initial condition
$\mathcal{X}(0)\in\mathbb{R}^{6}_{+}$, the stochastic system (5) admit a
unique, a.s. non-negative solution.

###### Proof.

The proof is based on similar arguments as in the proof of Theorem 2.1. ∎

To achieve our goals, we first study the conditions under which the terms
involving differentials of Brownian motion vanishes as $t\to\infty$. To do
this, let us consider a function $F:\mathbb{R}^{6}_{+}\to\mathbb{R}_{+}$.
Assume that $F:=F\left(\sum_{i=1}^{6}x_{i}\right)$ is an increasing twice
continuously differentiable function:

$1^{\circ}$.
$F^{\prime}=F^{\prime}_{x_{1}}=F^{\prime}_{x_{2}}=F^{\prime}_{x_{3}}=F^{\prime}_{x_{4}}=F^{\prime}_{x_{5}}=F^{\prime}_{x_{6}}$;

$2^{\circ}$.
$F^{\prime\prime}=F^{\prime\prime}_{x_{1}x_{1}}=F^{\prime\prime}_{x_{2}x_{2}}=F^{\prime\prime}_{x_{3}x_{3}}=F^{\prime\prime}_{x_{4}x_{4}}=F^{\prime\prime}_{x_{5}x_{5}}=F^{\prime\prime}_{x_{6}x_{6}}$,

with $F^{\prime}_{x_{k}}=\frac{\partial
F\left(\sum_{i=1}^{6}x_{i}\right)}{\partial x_{k}}$ and
$F^{\prime\prime}_{x_{k}x_{k}}=\frac{\partial^{2}F\left(\sum_{i=1}^{6}x_{i}\right)}{\partial{x_{k}}^{2}}$
for $k=1,...,6$.

Let $\sigma_{j},$ $j=3,...,8$ be continuous positive functions. Denote

$\mathcal{I}(t)=\int_{0}^{t}F^{\prime}(\mathcal{K}(s))\cdot\left[\sigma_{3}(s)S_{h}(s)dB_{3}(s)+\sigma_{4}(s)I_{h}(s)dB_{4}(s)+\sigma_{5}(s)Q_{h}(s)dB_{5}(s)\right.$
$\left.+\sigma_{6}(s)R_{h}(s)dB_{6}(s)+\sigma_{7}(s)S_{r}(s)dB_{7}(s))+\sigma_{8}(s)I_{r}(s)dB_{8}(s)\right],$

where $\mathcal{K}$ such as in (29), and we get the following result.

######  Theorem 5.2.

Let $\left\\{B_{j}(t):t\geq 0\right\\},j=3,4,5,6,7,8,$ be independent Brownian
processes, $\sigma_{j},$ $j=3,...,8$ be continuous positive functions, and for
$i=1,...,6,$ function $x_{i}F^{\prime}\left(\sum_{i=1}^{6}x_{i}\right)$,
$x=(x_{1},...,x_{6})\in\mathbb{R}^{6},$ be bounded. Assume

$\sum_{n=1}^{\infty}\frac{1}{2^{2n}}\int_{0}^{2^{n+1}}\Xi(s)\,ds<\infty,$ (45)

where $\Xi(t)=max\\{\sigma_{j}(t),\,t\geq 0,\,j=3,4,5,6,7,8\\}$. Then

$\underset{t\to\infty}{\lim}\frac{\mathcal{I}(t)}{t}=0\,\,\,\textit{a.s.}$
(46)

###### Proof.

Let us fix $\epsilon>0$. Using the Doob’s inequality

$\displaystyle\mathrm{P}\left\\{\underset{k\geq
n}{\sup}\frac{|\mathcal{I}(2^{k})|}{2^{k}}>\epsilon\right\\}\leq\mathrm{P}\left\\{\underset{t\geq
2^{n}}{\sup}\frac{|\mathcal{I}(t)|}{t}>\epsilon\right\\}\leq\sum_{k=n}^{\infty}\mathrm{P}\left\\{\underset{2^{k}\leq
t\leq 2^{k+1}}{\sup}\frac{|\mathcal{I}(t)|}{t}>\epsilon\right\\}$ (47) (48)
$\displaystyle\leq\sum_{k=n}^{\infty}\mathrm{P}\left\\{\underset{2^{k}\leq
t\leq
2^{k+1}}{\sup}\left(\frac{1}{t}\Big{|}\int_{0}^{t}F^{\prime}(\mathcal{K}(s))\sigma_{3}(s)S_{h}(s)dB_{3}(s)\Big{|}\right)>\frac{\epsilon}{6}\right\\}$
(49) (50)
$\displaystyle+\sum_{k=n}^{\infty}\mathrm{P}\left\\{\underset{2^{k}\leq t\leq
2^{k+1}}{\sup}\left(\frac{1}{t}\Big{|}\int_{0}^{t}F^{\prime}(\mathcal{K}(s))\sigma_{4}(s)I_{h}(s)dB_{4}(s)\Big{|}\right)>\frac{\epsilon}{6}\right\\}$
(51) (52)
$\displaystyle+\sum_{k=n}^{\infty}\mathrm{P}\left\\{\underset{2^{k}\leq t\leq
2^{k+1}}{\sup}\left(\frac{1}{t}\Big{|}\int_{0}^{t}F^{\prime}(\mathcal{K}(s))\sigma_{5}(s)Q_{h}(s)dB_{5}(s)\Big{|}\right)>\frac{\epsilon}{6}\right\\}$
(53) (54)
$\displaystyle+\sum_{k=n}^{\infty}\mathrm{P}\left\\{\underset{2^{k}\leq t\leq
2^{k+1}}{\sup}\left(\frac{1}{t}\Big{|}\int_{0}^{t}F^{\prime}(\mathcal{K}(s))\sigma_{6}(s)R_{h}(s)dB_{6}(s)\Big{|}\right)>\frac{\epsilon}{6}\right\\}$
(55) (56)
$\displaystyle+\sum_{k=n}^{\infty}\mathrm{P}\left\\{\underset{2^{k}\leq t\leq
2^{k+1}}{\sup}\left(\frac{1}{t}\Big{|}\int_{0}^{t}F^{\prime}(\mathcal{K}(s))\sigma_{7}(s)S_{r}(s)dB_{7}(s)ds\Big{|}\right)>\frac{\epsilon}{6}\right\\}$
(57) (58)
$\displaystyle+\sum_{k=n}^{\infty}\mathrm{P}\left\\{\underset{2^{k}\leq t\leq
2^{k+1}}{\sup}\left(\frac{1}{t}\Big{|}\int_{0}^{t}F^{\prime}(\mathcal{K}(s))\sigma_{8}(s)I_{r}(s)dB_{8}(s)\Big{|}\right)>\frac{\epsilon}{6}\right\\}$
(59) (60)
$\displaystyle\leq\frac{144}{\epsilon^{2}}\sum_{k=n}^{\infty}\frac{1}{2^{2k}}\int_{0}^{2^{k+1}}\mathrm{E}\Big{(}F^{\prime}(\mathcal{K}(s))^{2}\cdot(\sigma^{2}_{3}(s)S^{2}_{h}(s)+\sigma^{2}_{4}(s)I^{2}_{h}(s)+\sigma^{2}_{5}(s)Q^{2}_{h}(s)$
(61) (62)
$\displaystyle+\sigma^{2}_{6}(s)R^{2}_{h}(s)+\sigma^{2}_{7}S^{2}_{r}(s)+\sigma^{2}_{8}I^{2}_{r}(s)\Big{)}ds\leq$
(63) (64)
$\displaystyle\leq\frac{144}{\epsilon^{2}}\sum_{k=n}^{\infty}\frac{1}{2^{2k}}\int_{0}^{2^{k+1}}\Xi(s)\cdot\mathrm{E}\Big{(}F^{\prime}(\mathcal{K}(s))^{2}\cdot(S^{2}_{h}(s)+I^{2}_{h}(s)+Q^{2}_{h}(s)$
(65) (66)
$\displaystyle+R^{2}_{h}(s)+S^{2}_{r}(s)+I^{2}_{r}(s)\Big{)}ds\leq\frac{144}{\epsilon^{2}}\sum_{k=n}^{\infty}\frac{1}{2^{2k}}\int_{0}^{2^{k+1}}\Xi(s)\cdot
6\tilde{C}^{2}_{1}ds$ (67) (68)
$\displaystyle=\frac{864\cdot\tilde{C}^{2}_{1}}{\epsilon^{2}}\sum_{k=n}^{\infty}\frac{1}{2^{2k}}\int_{0}^{2^{k+1}}\Xi(s)ds,$
(69)

where
$\tilde{C}_{1}=\underset{x_{i}}{\sup}\left\\{x_{i}F^{\prime}\left(\sum_{i=1}^{6}x_{i}\right)\right\\}$.
Thus

$\mathrm{P}\left\\{\underset{k\geq
n}{\sup}\frac{\Big{|}\mathcal{I}(2^{k})\Big{|}}{2^{k}}>\epsilon\right\\}\leq\frac{864\cdot\tilde{C}^{2}_{1}}{\epsilon^{2}}\sum_{k=n}^{\infty}\frac{1}{2^{2k}}\int_{0}^{2^{k+1}}\Xi(s)ds$

for any $n\geq 0$ and $\epsilon>0$. By assumption, we have that

$\displaystyle\sum_{k=n}^{\infty}\frac{1}{2^{2k}}\int_{0}^{2^{k+1}}\Xi(s)ds\to\,\,0,\,\,\textnormal{as}\,\,n\to\infty.$

Since $\epsilon>0$ is arbitrary, we have the relations

$\mathrm{P}\left\\{\underset{t\to\infty}{\lim\sup}\frac{\Big{|}\mathcal{I}(t)\Big{|}}{t}>0\right\\}=0\Rightarrow\mathrm{P}\left\\{\underset{t\to\infty}{\lim\sup}\frac{\Big{|}\mathcal{I}(t)\Big{|}}{t}=0\right\\}=1.$

Therefore,

$\underset{t\to\infty}{\lim}\frac{\Big{|}\mathcal{I}(t)\Big{|}}{t}=0\,\,\,\,\textnormal{a.s.}$

i.e. relation (46) holds, and Theorem 5.2 is proved. ∎

######  Remark 5.3.

For the monkeypox model (3) wich has constant volatilities
$\sigma_{i},i=3,4,5,6,7,8,$ the conditions of the Theorem 5.2 are fulfilled
automatically. Indeed, if
$m_{\sigma}=max\\{\sigma_{3},\sigma_{4},\sigma_{5},\sigma_{6},\sigma_{7},\sigma_{8}\\},$
then

$\sum_{n=1}^{\infty}\frac{1}{2^{2n}}\int_{0}^{2^{n+1}}\Xi(s)\,ds=\sum_{n=1}^{\infty}\frac{1}{2^{2n}}\int_{0}^{2^{n+1}}m_{\sigma}\,ds=m_{\sigma}\sum_{n=1}^{\infty}\frac{1}{2^{n-1}}<\infty.$

Hence, we obtain (46) from the properties of Brownian motion.

Now let’s move to the study of the extinction of the disease. Consider the
operator $L$ which acts on the Lyapunov function
$\mathcal{V}:\mathbb{R}^{6}_{+}\to\mathbb{R}_{+}$ given by
$\mathcal{V}(\mathcal{X}(t))=F(\mathcal{K}(t))$ have the form

$\displaystyle
L\mathcal{V}(\mathcal{X}(t))=F^{\prime}(\mathcal{K}(t))\bigg{(}\theta_{h}(t)-\mu_{h}(t)(S_{h}(t)+I_{h}(t)+Q_{h}(t)+R_{h}(t))$
$\displaystyle-\delta_{h}(t)I_{h}(t)-(1-\theta(t))\delta_{h}(t)Q_{h}(t)+\theta_{r}(t)-\mu_{r}(t)(S_{r}(t)+I_{r}(t))-\delta_{r}(t)I_{r}(t)\bigg{)}$
$\displaystyle+F^{\prime\prime}(\mathcal{K}(t))\bigg{(}\sigma^{2}_{3}(t)S^{2}_{h}(t)+\sigma^{2}_{4}(t)I^{2}_{h}(t)+\sigma^{2}_{5}(t)Q^{2}_{h}(t)+\sigma^{2}_{6}(t)R^{2}_{h}(t)$
$\displaystyle+\sigma^{2}_{7}(t)S^{2}_{r}(t)+\sigma^{2}_{8}(t)I^{2}_{h}(t)\bigg{)}$
$\displaystyle+\frac{(1-p)^{2}\sigma^{2}_{1}(t)I^{2}_{r}(t)S^{2}_{h}(t)F^{\prime\prime}(\mathcal{K}(t))}{N^{2}_{h}(t)}+\frac{2(1-p)^{2}\sigma^{2}_{2}(t)I^{2}_{h}(t)S^{2}_{h}(t)F^{\prime\prime}(\mathcal{K}(t))}{N^{2}_{h}(t)}.$

The following theorem describes the parameters that mostly affect the rate of
the disease spread.

######  Theorem 5.4.

Assume that all conditions Theorem 5.2 hold and suppose

I.
$\underset{t\to\infty}{\lim}\frac{1}{t}\int_{0}^{t}\Big{|}\mathcal{V}^{\prime}(\mathcal{X}(s))\Big{|}(\theta_{h}(s)+\theta_{r}(s))ds=a\in[0;+\infty)$
a.s.;

II.
$\underset{t\to\infty}{\lim}\frac{1}{t}\int_{0}^{t}H_{1}(s)ds=b\in(0;+\infty)$
with $H_{1}(t)=\mu_{h}(t)\vee\mu_{r}(t)$;

III.
$\underset{t\to\infty}{\lim}\frac{1}{t}\int_{0}^{t}\left((2-\theta(s))\delta_{h}(s)+\delta_{r}(s)\right)ds=c\in[0;+\infty).$

IV.
$\lim\underset{t\to\infty}{\lim}\frac{1}{t}\int_{0}^{t}(\sigma^{2}_{1}(s)+2\sigma^{2}_{2}(s))ds=d\in[0;+\infty)$;
Then

$\underset{t\to\infty}{\lim}\frac{\mathcal{V}(\mathcal{X}(t))}{t}\leq
a+3\tilde{C}_{1}\cdot(2b+c)+(1-p)^{2}\cdot\tilde{C}_{2}\cdot d\,\,\
\textit{a.s.},$

where $\tilde{C}_{1}=\sup\left\\{\mathcal{V}^{\prime}x_{1}\right\\}$,
$\tilde{C}_{2}=\sup\\{\mathcal{V}^{\prime\prime}x^{2}_{1}\\}$.

###### Proof.

Applying Itô formula to $\mathcal{V},$ we obtain

$\displaystyle d(\mathcal{V}(\mathcal{X}(t)))$ $\displaystyle=$ $\displaystyle
L\mathcal{V}(\mathcal{X}(t))dt+\mathcal{V}^{\prime}(\mathcal{X}(t))\bigg{(}\frac{(p-1)\sigma_{1}I_{r}(t)S_{h}(t)}{N_{h}(t)}dB_{1}(t)$
$\displaystyle+\sigma_{3}(t)S_{h}(t)dB_{3}(t)+\sigma_{4}(t)I_{h}(t)dB_{4}(t)+\sigma_{5}(t)Q_{h}(t)dB_{5}(t)$
$\displaystyle+\sigma_{6}(t)R_{h}(t)dB_{6}(t)+\sigma_{7}(t)S_{r}(t)dB_{7}(t))+\sigma_{8}(t)I_{r}(t)dB_{8}(t)\bigg{)},$

or in integral form

$\mathcal{V}(\mathcal{X}(t))=\mathcal{V}(\mathcal{X}(0))+\int_{0}^{t}L\mathcal{V}(\mathcal{X}(s))ds+(p-1)\sigma_{1}\int_{0}^{t}\frac{\mathcal{V}^{\prime}(\mathcal{X}(s))I_{r}(s)S_{h}(s)}{N_{h}(s)}dB_{1}(s)+\mathcal{I}(t).$

By Theorem 5.2, we get

$\underset{t\to\infty}{\lim}\frac{\mathcal{I}(t)}{t}=0\,\,\,\textit{a.s.}$

Since $S_{h}<N_{h}$, and $x_{i}\mathcal{V}^{\prime}$, $i=1,...,6,$ is bounded,
then the strong low of large number gives

$\underset{t\to\infty}{\lim}\frac{(p-1)\sigma_{1}}{t}\int_{0}^{t}\frac{\mathcal{V}^{\prime}(\mathcal{X}(s))I_{r}(s)S_{h}(s)}{N_{h}(s)}dB_{1}(s)=(p-1)\sigma_{1}\tilde{C}_{1}\underset{t\to\infty}{\lim}\frac{B_{1}(t)}{t}=0.$

While (45) holds, we have

$\displaystyle\mathcal{I}_{1}(t)=\int_{0}^{t}\mathcal{V}^{\prime\prime}(\mathcal{X}(s))(\sigma^{2}_{3}(s)S^{2}_{h}(s)+\sigma^{2}_{4}(s)I^{2}_{h}(s)+\sigma^{2}_{5}(s)Q^{2}_{h}(s)$
$\displaystyle+\sigma^{2}_{6}(s)R^{2}_{h}(s)+\sigma^{2}_{7}(s)S^{2}_{r}(s)+\sigma^{2}_{8}(s)I^{2}_{r}(s))ds$
$\displaystyle\leq\int_{0}^{t}\Xi(s)\cdot
F^{\prime\prime}(\mathcal{K}(s))(S^{2}_{h}(s)+I^{2}_{h}(s)+Q^{2}_{h}(s)+R^{2}_{h}(s)+S^{2}_{r}(s)+I^{2}_{r}(s))ds$
$\displaystyle\leq 6\tilde{C}_{2}\int_{0}^{t}\Xi(s)ds.$

Therefore
$\displaystyle\underset{t\to\infty}{\lim}\frac{\mathcal{I}_{1}(t)}{t}=0\,\,\textnormal{a.s.}$
Since
$\dfrac{S^{2}_{h}(s)\mathcal{V}^{\prime\prime}(\mathcal{X}(s))}{N^{2}_{h}(s)}\leq\mathcal{V}^{\prime\prime}(\mathcal{X}(s))$.
It is easy to see, that

$\displaystyle\mathcal{I}_{2}(t)=(1-p)^{2}\int_{0}^{t}(\sigma^{2}_{1}(s)I^{2}_{r}(s)+2\sigma^{2}_{2}(s)I^{2}_{h}(s))\frac{S^{2}_{h}(s)\mathcal{V}^{\prime\prime}(\mathcal{X}(s))}{N^{2}_{h}(s)}ds$
$\displaystyle\leq(1-p)^{2}\cdot\tilde{C}_{2}\int_{0}^{t}(\sigma^{2}_{1}(s)+2\sigma^{2}_{2}(s))ds.$

Next estimation,

$\displaystyle\mathcal{I}_{3}(t)=\Big{|}\int_{0}^{t}\mathcal{V}^{\prime}(\mathcal{X}(s))\bigg{(}\theta_{h}(s)+\theta_{r}(s)-\mu_{h}(s)(S_{h}(s)+I_{h}(s)+Q_{h}(s)+R_{h}(s))$
$\displaystyle-\mu_{r}(s)(S_{r}(s)+I_{r}(s))-\delta_{h}(s)I_{h}(s)-(1-\theta(s))\delta_{h}(s)Q_{h}(s)-\delta_{r}(s)I_{r}(s)\bigg{)}ds\Big{|}$
$\displaystyle\leq\int_{0}^{t}\Big{|}\mathcal{V}^{\prime}(\mathcal{X}(s))\Big{|}(\theta_{h}(s)+\theta_{r}(s))ds+6\tilde{C}_{1}\int_{0}^{t}H_{1}(s)ds$
$\displaystyle+3\tilde{C}_{1}\int_{0}^{t}\left(\delta_{h}(s)+(1-\theta(s))\delta_{h}(s)+\delta_{r}(s)\right)ds.$

Here $H_{1}(t)=\mu_{h}(t)\vee\mu_{r}(t)$. Finally,

$\displaystyle\underset{t\to\infty}{\lim}\frac{\mathcal{V}(\mathcal{X}(t))}{t}$
$\displaystyle=$
$\displaystyle\underset{t\to\infty}{\lim}\frac{\mathcal{V}(\mathcal{X}(0))}{t}+\underset{t\to\infty}{\lim}\frac{1}{t}\int_{0}^{t}L\mathcal{V}(\mathcal{X}(s))ds$
$\displaystyle+\underset{t\to\infty}{\lim}\frac{(p-1)\sigma_{1}}{t}\int_{0}^{t}\frac{\mathcal{V}^{\prime}(\mathcal{X}(s))I_{r}(s)S_{h}(s)}{N_{h}(s)}dB_{1}(s)+\underset{t\to\infty}{\lim}\frac{\mathcal{I}(t)}{t}$
$\displaystyle=\underset{t\to\infty}{\lim}\frac{1}{t}\int_{0}^{t}\Big{|}\mathcal{V}^{\prime}(\mathcal{X}(s))\Big{|}(\theta_{h}(s)+\theta_{r}(s))ds+\underset{t\to\infty}{\lim}\frac{6\tilde{C}_{1}}{t}\int_{0}^{t}H_{1}(s)ds$
$\displaystyle+(1-p)^{2}\tilde{C}_{2}\cdot\underset{t\to\infty}{\lim}\frac{1}{t}\int_{0}^{t}(\sigma^{2}_{1}(s)+2\sigma^{2}_{2}(s))ds$
$\displaystyle+3\tilde{C}_{1}\cdot\underset{t\to\infty}{\lim}\frac{1}{t}\int_{0}^{t}\left(\delta_{h}(s)+(1-\theta(s))\delta_{h}(s)+\delta_{r}(s)\right)ds$
$\displaystyle\leq a+6\tilde{C}_{1}\cdot b+(1-p)^{2}\cdot\tilde{C}_{2}\cdot
d+3\tilde{C}_{1}\cdot c.$

Thus, Theorem 5.4 is proved. ∎

######  Remark 5.5.

The verification of assumption I. is quite tedious, as it depends on the
stochastic processes. Therefore, it is advisable to use the following
sufficient conditions:
$\underset{x}{\sup}\Big{|}\mathcal{V}^{\prime}(x)\Big{|}<\infty$ and
$\underset{t\to\infty}{\lim}\frac{1}{t}\int_{0}^{t}(\theta_{h}(s)+\theta_{r}(s))ds=A\in[0;+\infty)$
a.s.

######  Corollary 5.6.

Let $I_{h},Q_{h},I_{r}$ satisfy the system (5). Assume

$\sum_{n=1}^{\infty}\frac{1}{2^{2n}}\int_{0}^{2^{n+1}}\tilde{\Xi}(s)ds<\infty,$

where $\tilde{\Xi}(t)=max\\{\sigma_{i}(t),t\geq 0,i=2,4,5,6,\\}$. If

$\underset{t\to\infty}{\lim}\sup\frac{1}{t}\int_{0}^{t}\left((1-p)(\eta_{1}(s)+\eta_{2}(s))+\eta_{3}(s)\right)ds<\underset{t\to\infty}{\lim}\sup\frac{1}{t}\int_{0}^{t}(M_{1}(s)+M_{2}(s))ds,$
(70)

for $M_{1}(t)=\min\\{\mu_{h}(t),\mu_{r}(t)\\}$ and
$M_{2}(t)=\min\\{\delta_{h}(t),\delta_{r}(t)\\}$, $t\geq 0$, then the classes
$I_{h},Q_{h},I_{r}$ will die out exponentially with probability one, i.e.,

$\displaystyle\underset{t\to\infty}{\lim\sup}\frac{\ln(I_{h}(t)+Q_{h}(t)+I_{r}(t))}{t}$
$\displaystyle<$ $\displaystyle 0\,\,\textnormal{a.s.}$

###### Proof.

The proof follows from the Itô formula for $\ln(I_{h}+Q_{h}+I_{r})$ and
Theorem 5.2. ∎

######  Remark 5.7.

The result of Corollary 5.6 generalizes the result of Theorem 4.1 for the
monkeypox model (3) where parameters are in the form of constants. It should
be noted that assumption (70) corresponds to $\mathcal{R}_{0}<1$ in Theorem
4.1 in cases of constant.

The rate of expansion shows how rapidly the monkeypox is spreading. The growth
rate over time is a good indicator of the severity and speed of the epidemic.
A faster rate of increase indicates that the population’s health is declining
more quickly, putting more strain on the healthcare systems.

##  6 Conclusion and Final Comments

In the present research work, we have introduced and analyzed a stochastic
monkeypox epidemic model with cross-infection between animals and humans. Our
results discuss the well-posedness of the model to secure that the simulations
and prediction will be meaningfull. Specifically, we have seen global
existence, positivity and uniqueness of solution. Subsequently, we have
expounded results for stochastic ultimate boundedness and permanence and
results related to the long-run behavior via exploring asymptotic properties
of the model in a time-depended variation.

By anticipating and preparing for possible outbreaks, public health
professionals may utilize these forecasts to lessen the burden of the disease
with early and focused actions. The compartment modeling technique can
incorporate compartments for immune people to better understand the fading and
development of immunity over time and the immunity of the population. It may
also be used to find possible amplification hosts, such rodents in the case of
monkeypox, and reservoirs. This information is critical for understanding the
ecology of the disease, directing monitoring efforts, and putting control
measures in place. Randomness in the spread of disease given by the uncertain
environment is taken into account. This facilitates the creation of
probabilistic forecasts, increases prediction reliability, and directs
decision-making in the face of uncertainty.

## Statements and Declarations

Funding. The present research is carried out within the frame and support of
the ToppForsk project nr. 274410 of the Research Council of Norway with the
title STORM: Stochastics for Time-Space Risk Models, and the project ”Dynamics
for Transmission of Monkeypox: Biological and Probabilistic Behaviour”
supported by the Higher Education Commission of Pakistan, and supported by the
MSCA4Ukraine grant (AvH ID:1233636), which is funded by the European
Commission.

Acknowledgements. The first author would like to thank the University of Swat,
KP Pakistan, for providing the opportunity to complete nine months postdoc
with stay at the University of Oslo, Norway, and CIMPA, Nice, France to have
provided the grounds to meet and form the present scientific collaboration.

## References

  * [1] M.J. Keeling, P. Rohani, Modeling Infectious Diseases in Humans and Animals. _Princeton University Press_ , (2008).
  * [2] R.M. Anderson, R.M. May, Infectious Diseases of Humans: Dynamics and Control._Oxford University Press_ , (1992).
  * [3] O.N. Bjornstad, B.T. Grenfell, Noisy Clockwork: Time series analysis of population fluctuations in animals. _Science_ , 293, 5539 (2001): 638–643, DOI: 10.1126/science.1062226
  * [4] A.J. Black, A.J. McKane, Stochastic formulation of ecological models and their applications. _Trends in Ecology and Evolution_ , 27(6) (2012), 337–345.
  * [5] R.P. Agarwal, Q. Badshah, G. ur Rahman, S. Islam, Optimal control & dynamical aspects of a stochastic pine wilt disease model. _Journal of the Franklin Institute_ , 356 (2019),33914025.
  * [6] G. ur Rahman, Q. Badshah, Ravi P. Agarwal, S. Islam, Ergodicity & dynamical aspects of a stochastic childhood disease model. _Mathematics and Computers in Simulation_ ,182 (2021) 738–764.
  * [7] K.N. Durski, et. al.: _Emergence of monkeypox-west and central Africa_ , 1970-2017, Morbidity and mortality weekly report, 67, (2018).
  * [8] Z. Jezek, M. Szczeniowski, K.M. Paluku, M. Mutombo, B. Grab, Human monkeypox: confusion with chickenpox. _Acta Tropica_ , 45, (1988).
  * [9] E. Alakunle, U. Moens, G. Nchinda, M. Okeke, Monkeypox Virus in Nigeria: Infection Biology, Epidemiology, and Evolution. _Viruses_ , 12, (2020).
  * [10] A. Kantele, K. Chickering, O. Vapalahti, A. Rimoin, Emerging diseases-the monkeypox epidemic in the Democratic Republic of the Congo. _Clinical Microbiology and Infection_ , 22, (2016).
  * [11] J.O. Wolff and P.W. Sherman, Rodent Societies: An Ecological and Evolutionary Perspective. _University of Chicago Press_ , (2007).
  * [12] B.G. Meerburg1, G.R. Singleton, A. Kijlstra, Rodent-borne diseases and their risks for public health. _Critical Reviews in Microbiology_ , 35(3), (2009), 221–270.
  * [13] Castillo-Chavez, S. Blower, P. Driessche, D. Kirschner, A. Yakubu, Mathematical Approaches for Emerging and Re-Emerging Infectious Diseases:An Introduction. _Springer_ , Berlin, 126, (2002).
  * [14] J.P. LaSalle, The Stability of Dynamical Systems. _SIAM, Philadelphia_ , 25, (1976), 176.
  * [15] R. Khasminskii, Stochastic Stability of Differential Equations. _Springer_ , 2012.
  * [16] F.C. Klebaner, Introduction to Stochastic Calculus with Applications. _Imperial College Press_ , (2005).
  * [17] B. Al-Hdaibat, M.Z. Meetei, I. Ahmad, N. Althobaiti, M. A. Safi, M.A. Khan, M.B. Riaz, The dynamics of monkeypox disease under $\Psi$-Hilfer fractional derivative: Application to real data. _Results in Physics_ , 55 (2023) 107127.
  * [18] M.A. Khan, M.Z. Meetei, K. Shah, T. Abdeljawad, M.Y. Alshahrani, Modeling the monkeypox infection using the Mittag-Leffler kernel. _Open Physics_ , 21(1), (2023), 20230111.
  * [19] C.P. Bhunu, S. Mushayabasa, Modelling the transmission dynamics of pox-like infections. _IAENG International Journal of Applied Mathematics_ , 41(2), 09 (2011).
  * [20] S. Usman, I.I. Adamu, Modeling the transmission dynamics of the Monkeypox virus infection with treatment and vaccination Interventions. _Journal of Applied Mathematics and Physics_ , 5, (2017), 2335–2353.
  * [21] X. Mao, Stochastic Differential Equations and Applications. _Horwood Publishing Chichester, UK_ , (1997).
  * [22] L. Arnold Stochastic, Differential Equations: Theory and Applications. _Wiley, New York_ , (1972).
  * [23] I. I. Gihman, A. V. Skorohod, Stochastic Differential Equations. _Springer-Verlag, Berlin-Heidelberg-New York_ , (1972).
  * [24] Y. Cai, X. Wang, W. Wang, M. Zhao, Stochastic dynamics of a SIRS epidemic model with ratio-dependent incidence rate. _Abstract and Applied Analysis_ , (2013) 415–425.
  * [25] X. Mao, G. Marion, and E. Renshaw, Environmental Brownian noise suppresses explosions in population dynamics. _Stochastic Processes and their Applications_ , 97(1)(2002), 95–110.
  * [26] N.T. Dieu, D H. Nguyen, N.H. Du and G. Yin, Classification of Asymptotic Behaviour in a Stochastic SIR Model. _SIAM Journal on Applied Dynamical System_ , 15 (2016), 1062–1084.
  * [27] D. Jiang, C. Ji, N. Shi, J. Yu, The long time behaviour of DI SIR epidemic model with stochastic perturbation. _Journal of Mathematical Analysis and Applications_ , 372 (2010), 162-180.
  * [28] Y. G. Lin and D. Q. Jiang, Long-time behaviour of a perturbed SIR model by white noise._Discrete Continuous Dynamical System Series B_ , 18 (2013), 1873–1887.
  * [29] V. V. Buldygin and O.A. Tymoshenko, On the exact order of growth of solutions of stochastic differential equations with time-dependent coefficients. _Theory of Stochastic Processes_ , 16 (2010), 12-22.
  * [30] O.I. Klesov, O.A. Tymoshenko, Almost sure asymptotic properties of solutions of a class of non-homogeneous stochastic differential equations. _Modern Mathematics and Mechanics. Fundamentals Problems and Challenges_ , Switzerland: Springer, (2019) 97-114.