# Single-species population models with
age structure and psychological effect in a polluted environment

Jiawei Wang  Ruiyang Zhou  Fengying Wei<EMAIL_ADDRESS>
###### Abstract

This paper considers a single-population model with age structure and
psychological effects in a polluted environment. We divide the single
population into two stages of larval and adult structure. The model uses
Logistic input, and the larvae are converted into adult bodies by constant
ratio. We only consider adulthood. The role of psychological effects makes the
contact between adult and environmental toxins a functional form, while the
contact between larvae and environmental toxins is linear.

For the deterministic model embodied as a nonlinear time-varying system, we
discuss the asymptotic stability of the system by Lyapunov one-time
approximation theory, and give a sufficient condition for stability to be
established.

Considering that the contact rate between biological and environmental toxins
in nature is not always constant, we make the contact rate interfere with
white noise, and then modify the contact rate into a stochastic process, thus
establishing a corresponding random single-population model. According to Itô
formula and Lyapunov in the function method, we first prove the existence of
globally unique positive solutions for stochastic models under arbitrary
initial conditions, and then give sufficient conditions for weak average long-
term survival and random long-term survival for single populations in the
expected sense.

###### keywords:

polluted environment , stochastic single-species population models , age
structure , psychological effect , stability , persistence

††journal: Mathematical Biosciences

###### Proof.

[a] organization=School of Mathematics and Statistics, Fuzhou University,
city=Fuzhou, postcode=350116, state=Fujian, country=China

[b] organization=Institute of Fundamental and Frontier Sciences, University of
Electronic Science and Technology of China, city=Chengdu, postcode=611731,
state=Sichuan, country=China

[c] organization=Key Laboratory of Operations Research and Control of
Universities in Fujian, Fuzhou University, city=Fuzhou, postcode=350116,
state=Fujian, country=China

##  1 Introduction

With the development of society, the problem of environmental pollution caused
by human activities and industrial production has become increasingly serious,
environmental pollution has always been the most threatening social and
ecological problem. All kinds of pollutants discharged into the environment
have a great impact on the normal survival of all kinds of organisms in the
environment. All kinds of organisms, including human beings, are threatened by
all kinds of poisons in the environment. At present, many scholars have
established mathematical models to describe it.

Since Hallam et al. [1, 2, 3, 4, 5], scholars have paid more and more
attention to this problem. In 1983, Hallam et al. [1] first proposed a
classical deterministic mathematical model for the survival of a single
population in a polluted environment. Other scholars continue to study the
environmental pollution model [6, 7, 8, 9, 10, 11]. Srinivasu [12] considered
the impact of biological death on the concentration of environmental toxins,
studied the problem of population’s long-term survival and extinction, and
considered to ensure the population’s long-term survival by increasing the
amount of environmental toxins removed.

The growth of species is basically accompanied by a process of development,
i.e. from adulthood to old age, from immature to mature, from young individual
to adult individual. Each stage of the growth process often shows different
characteristics, such as the young individual does not have the ability of
reproduction and predation. Due to the limitation of intelligence development
and survival experience, the young individual often does not have the ability
of self avoiding danger. Compared with other stages of population, the ability
of survival and resources competition of young individuals is weak. Young
individuals are often easy to die, and unable to complete large-scale
migration in space.

The adult population has the ability of reproduction and predation, a high
level of intellectual development. It is also superior in survival experience
and can actively avoid danger, thus has a strong survival ability and the
ability to compete with other populations for limited survival resources. Some
of them experience that one side is injured or even killed in order to compete
for spouse. There are significant differences in the behavior of the same
stage. In addition, there is always a relationship of mutual transformation
between organisms in different stages, which will have an impact on the
extinction and long-term survival of organisms. Therefore, it is of practical
significance to study the population model with stage structure, which has
also attracted the attention of many mathematicians and biologists.

The related work of stage structure population model can be traced back to the
three-stage structure self feeding model proposed by Landahl and Hansen [13],
in the 1970s, and the two-stage structure stochastic model proposed by
Tognetti [14]. Later, many scholars also proposed different stage structure
models [15]. But until 1990, Aiello and Freedman [16] proposed the single
group delay model with two-stage structure. The related research only then
truly ushers in the upsurge [17, 18, 19, 20, 21, 22, 23].

For vertebrates with a certain level of intelligence, when life safety is
threatened by diseases, predators, environmental pollution and other external
threats, there are often a series of psychological effects such as fear.
Compared with other creatures with a low level of intelligence, psychological
effects affect the daily behavior of the whole population to a certain extent.

In recent years, many scholars have tried to make the mathematical models to
depict the influence of psychological effects. Many of them study human groups
directly [24, 25, 26, 27, 28, 29, 30], and some of them study other creatures.
In 2018, Kumar [31], studied the prey predator model of Allee effect induced
by fear. In 2019, Lan et al. [32], considered the long-term behavior of a
single population model with psychological effects under the pulse input of
environmental toxins. Assuming that the internal growth rate of a single
population is interfered by white noise, a corresponding stochastic model was
established.

In this paper, we mainly consider vertebrates with sensory organs and highly
differentiated nervous system. Compared with organisms without this structure
in nature, sensory organs can interpret information in dynamic environment.
For vertebrates, such as fish and birds, we can detect and process information
in the physical world. When the environment is polluted, these organisms
either bear the toxin concentration of surrounding environment or escapes from
the living area when the psychological effect plays a role. Combined with the
age structure, we naturally think that the intelligence level of young
individuals is lower than that of adult individuals. Therefore, when
considering the psychological effect, we should consider the two groups
separately. Based on the above characteristics, we will study the single
population with age structure and psychological effect in the polluted
environment dynamic behavior.

In the first part of this paper, based on the previous studies, a single
population model with age structure and psychological effects is established.
The model is a nonlinear time-varying system. According to Lyapunov’s first
approximation theory, we study the related properties of the original system
by studying the stability of the zero solution of the approximate system.

In the second part of this paper, we consider that the contact rate between
organisms and environmental poisons is disturbed by white noise, and then we
establish a stochastic model. At the same time, we give the proof of the weak
average persistence and stochastic persistence of the single population under
certain conditions.

##  2 Stability analysis of deterministic model

###  2.1 Deterministic model

A single population model with psychological effects in polluted environment
is proposed [29]

$\dot{x}(t)=x(t)(b-d)-cx^{2}(t)-\alpha x(t)C_{o}(t)-\lambda
x(t)g(C_{e},\beta),t\in[0,\infty)$ (1)

The variables of the model are $x=x(t)$, which means the number of population
at time $t$, $C_{o}=C_{o}(t)$, the degree of individual toxins at time $t$,
$C_{e}=C_{e}(t)$, the concentration of toxins in living environment at time
$t$, $b,d,c$ represent the natural birth rate, mortality rate and biological
intra-specific restriction factors when the toxin is zero. $\alpha$ represent
the response intensity of organism growth to the toxin. $\lambda$ represent
the contact rate between organism and environmental toxicant. $\beta$ refers
to the inhibition factors or psychological effects of the population in the
polluted environment. To some extent, $\beta$ describes the sensitivity of the
population in the polluted environment. $\alpha x(t)C_{o}(t)$ describes the
removal amount of endotoxin in the organism at time $t$. The toxic effect of
the polluted environment is represented by $g(C_{e},\beta)$ [26]

$g(C_{e},\beta)=\frac{C_{e}}{1+\beta C_{e}^{2}}$ (2)

also known as contact rate or psychological effect. Age structure is
introduced, and the specific model is as follows

$\begin{array}[]{lll}\dot{x}(t)=by(t)-d_{1}x(t)-\gamma
x(t)-c_{1}x^{2}(t)-\alpha_{1}x(t)C_{o}(t)-\lambda_{1}x(t)g_{1}(C_{e}(t),\beta),\\\\[8.53581pt]
\dot{y}(t)=\gamma
x(t)-d_{2}y(t)-c_{2}y^{2}(t)-\alpha_{2}y(t)C_{o}(t)-\lambda_{2}y(t)g_{2}(C_{e}(t),\beta),\end{array}$
(3)

$x,y$ refers to the number of young population and adult population, both of
which are related to time $t$. In this paper, we stipulate that the values of
$x$ and $y$ are all in the range of $[1,\infty)$. $\gamma$ represents the
transformation rate from larva to adult, and the meaning of other symbols is
consistent with the original. The toxin action in the polluted environment at
time $t$ is indicated by $g_{i}(C_{e},\beta),i=1,2$. Because there are
differences in intelligence level and survival experience between larva and
adult, we choose the following functions to describe the effect of
environmental toxin respectively

$g_{1}(C_{e},\beta)=C_{e},g_{2}(C_{e},\beta)=\frac{C_{e}^{p}}{1+\beta
C_{e}^{q}},1\leqslant p\leqslant q,p,q\in\mathbb{N}_{+}$ (4)

It is assumed that the change rule of endotoxin concentration meets the
following requirements

$\dot{C_{o}}(t)=kC_{e}(t)-gC_{o}(t)-mC_{o}(t)-bC_{o}(t)$ (5)

Among them, $kC_{e}(t)$ is the ratio of biological absorption of environmental
toxins at the time of $t$. $gC_{o}(t)$ is the ratio of biological elimination
of toxins at the time of $t$. $mC_{o}(t)$ is the ratio of biological
purification of toxins at the time of $t$. $bC_{o}(t)$ is the ratio of toxins
lost by new individuals at the time of $t$. We always assume that the
environmental capacity is large enough. We also assume that the toxins in new
individuals should be very small, and the accumulation of toxins from birth to
adulthood is very small, so the effect of toxin concentration on the density
of a single population can be ignored. The change rate of toxin concentration
in the environment is as follows

$\dot{C_{e}}(t)=u_{e}(t)-hC_{e}(t),$ (6)

$hC_{e}(t)$ is the ratio of toxin reduction caused by environmental self
purification at time $t$. $hC_{e}(t)$ is the concentration of exogenous input
toxin entering the environment at time $t$, and it is assumed to be a bounded
non-negative differentiable function of time $t$.

A single population model with age structure and psychological effects in the
polluted environment is established:

$\left\\{\begin{array}[]{ll}\dot{x}(t)=by(t)-d_{1}x(t)-\gamma
x(t)-c_{1}x^{2}(t)-\alpha_{1}x(t)C_{o}(t)-\lambda_{1}x(t)g_{1}(C_{e}(t),\beta),\\\\[8.53581pt]
\dot{y}(t)=\gamma
x(t)-d_{2}y(t)-c_{2}y^{2}(t)-\alpha_{2}y(t)C_{o}(t)-\lambda_{2}y(t)g_{2}(C_{e}(t),\beta),\\\\[8.53581pt]
\dot{C_{o}}(t)=kC_{e}(t)-gC_{o}(t)-mC_{o}(t)-bC_{o}(t),\\\\[8.53581pt]
\dot{C_{e}}(t)=u_{e}(t)-hC_{e}(t).\end{array}\right.$ (7)

The study of Bronstein[33], Cantalupo[34], Castro[35], Dzieweczynskia[36]
showed that Betta splendens Regan was very aggressive in courtship, and the
death of adult was mainly caused by intra-specific competition. Hagvar [37]
and Uka [38] showed that the death of adult was mainly caused by intra-
specific competition.

Therefore, according to these biological examples, When $d_{2}y(t)\ll
c_{2}y^{2}(t)$, so $d_{2}y(t)+c_{2}y^{2}(t)\approx c_{2}y^{2}(t)$, we modify
the model (7) and rewrite it as

$\left\\{\begin{array}[]{lll}\vspace{1ex}\dot{x}(t)=by(t)-dx(t)-\gamma
x(t)-c_{1}x^{2}(t)-\alpha_{1}x(t)C_{o}(t)-\lambda_{1}x(t)C_{e}(t),\\\
\vspace{1ex}\displaystyle\dot{y}(t)=\gamma
x(t)-c_{2}y^{2}(t)-\alpha_{2}y(t)C_{o}(t)-\frac{\lambda_{2}y(t)C_{e}^{p}(t)}{1+\beta
C_{e}^{q}(t)},\\\
\vspace{1ex}\dot{C_{o}}(t)=kC_{e}(t)-gC_{o}(t)-mC_{o}(t)-bC_{o}(t),\\\
\vspace{1ex}\dot{C_{e}}(t)=u_{e}(t)-hC_{e}(t).\end{array}\right.$ (8)

###  2.2 Stability analysis of nonlinear time-varying system

From (8), we can get

$\left\\{\begin{array}[]{lll}\vspace{1ex}C_{o}(t)=\displaystyle
k\int_{0}^{t}C_{e}(s)e^{-(g+m+b)(t-s)}\mathrm{d}s+C_{o}(0)e^{-(g+m+b)t},\\\
\vspace{1ex}C_{e}(t)=\displaystyle\int_{0}^{t}u_{e}(s)e^{-h(t-s)}\mathrm{d}s+C_{e}(0)e^{-ht}.\end{array}\right.$
(9)

Denote $X(t)=(x(t),y(t))^{T}$, (9) could be described as follows

$\dot{X}(t)=A(t)X(t)+G(X(t)),$ (10)

where

$\begin{array}[]{l}\vspace{1ex}A(t)=\mbox{$\left(\begin{array}[]{ccc}-d-\gamma-\alpha_{1}C_{o}(t)-\lambda_{1}C_{e}(t)&b\\\
\gamma&-\alpha_{2}C_{o}(t)-\frac{\lambda_{2}C_{e}^{p}(t)}{1+\beta
C_{e}^{q}(t)}\end{array}\right)$},\\\
\vspace{1ex}G(X(t))=\mbox{$\left(\begin{array}[]{lll}-c_{1}x^{2}(t)\\\
-c_{2}y^{2}(t)\end{array}\right)$}\end{array}$ (11)

Omitting the high-order term $G(x(t))$, obtaining the approximate linear
system

$\dot{X}(t)=A(t)X(t),$ (12)

linear system (12) is called the first order approximate equation of (10).
According to Lyapunov’s first approximation theory [39], if $G(X(t))$ is a
higher order infinitesimal of $X(t)$ in the neighborhood of $X=(0,0)^{T}$, the
stability of (10) can often be studied by using the stability of linear system
(12). In fact, we can easily get

$\lim_{t\to\infty}\frac{\|G(X(t))\|}{\|X(t)\|}=0,$ (13)

So next we discuss the stability of linear system (12).

######  Theorem 1.

If $|A(t)|>0$, there exists an $\varepsilon$ such that
$|\dot{a}_{ij}(t)|\leqslant\varepsilon$, then the trivial solution of model
(12) is globally asymptotically stable.

###### Proof.

Let

$A(t)=(a_{ij}(t))_{2\times
2}=\mbox{$\left(\begin{array}[]{ccc}-\Gamma_{1}(t)&b\\\
\gamma&-\Gamma_{2}(t)\end{array}\right)$}$ (14)

where $a_{ij}(t)$ are differentiable with $|a_{ij}(t)|\leqslant
b\vee(d+\gamma+\alpha_{1}+\lambda_{1})\vee(\alpha_{2}+\lambda_{2})\vee\gamma\geqslant
a$. The characteristic equation of $A(t)$ gives the following roots

$\lambda(A(t))=\frac{-(\Gamma_{1}(t)+\Gamma_{2}(t))\pm\sqrt{(\Gamma_{1}(t)-\Gamma_{2}(t))^{2}+4b\gamma}}{2}.$
(15)

When $|A(t)|>0$, each eigenvalue of $A(t)$ has a negative real part, i.e.
$\mbox{Re}\lambda(A(t))<0$. Rewrite the characteristic equation as follows

$f_{A}(\lambda)=\lambda^{2}+p_{1}(t)\lambda+p_{2}(t)=0,$ (16)

where

$p_{1}(t)=-\sum_{i=1}^{2}a_{ii}(t)=\Gamma_{1}(t)+\Gamma_{2}(t),\quad
p_{2}(t)=|A(t)|=\Gamma_{1}(t)\Gamma_{2}(t)-b\gamma.$ (17)

By Routh-Hurwitz theorem,

$\Delta_{1}(t)=p_{1}(t)>0,\quad\Delta_{2}(t)=p_{1}(t)p_{2}(t)>0.$ (18)

Let

$x=x_{1},y=x_{2},$ (19)

to find $v(x_{1},x_{2})$ satisfying

$\sum_{i=1}^{2}\frac{\partial v}{\partial
x_{i}}\sum_{j=1}^{2}a_{ij}x_{j}=-2\prod_{i=1}^{2}\Delta_{i}(t)\sum_{j=1}^{2}x_{j}^{2}.$
(20)

According to Barabashin formula,

$v(x_{1},x_{2})=\frac{\Delta_{1}(t)\Delta_{2}(t)}{\Delta(t)}\left|\begin{array}[]{cccc}0&x_{1}^{2}&2x_{1}x_{2}&x_{2}^{2}\\\
1&a_{11}&a_{21}&0\\\ 0&a_{12}&a_{11}+a_{22}&a_{21}\\\
1&0&a_{12}&a_{22}\end{array}\right|,$ (21)
$\Delta(t)=\left|\begin{array}[]{cccc}a_{11}&a_{21}&0\\\
a_{12}&a_{11}+a_{22}&a_{21}\\\ 0&a_{12}&a_{22}\end{array}\right|,$ (22)

which could be written as

$v(x_{1},x_{2})=C(t)\sum_{i,j=1}^{2}v_{ij}(t)x_{i}x_{j},$ (23)

with $v_{ij}(t)=v_{ji}(t),C(t)=\frac{\Delta_{1}(t)\Delta_{2}(t)}{\Delta(t)}$,
where $v_{ij}$ is the determinant obtained by exchanging the first column and
the $2x_{i}x_{j}$-column ($x_{i}^{2}$-column) and then removing the first row
and the first column in the exchanged determinant. Therefore

$v(x_{1},x_{2})=\sum_{i,j=1}^{2}V_{ij}(t)x_{i}x_{j},$ (24)

with $V_{ij}(t)=C(t)v_{ij}(t)$ and $V_{ij}(t)=V_{ji}(t)$. Next, we will prove
that function $v(x_{1},x_{2})$ is positive definite, given quadratic form

$W=-\prod_{i=1}^{2}\Delta_{i}(t)\sum_{j=1}^{2}x_{j}^{2}$ (25)

From (20), we can find the uniquely determined $v(x_{1},x_{2})$, so the
$v(x_{1},x_{2})$ expressed by (24) should be consistent with the Lyapunov
function in reference [40].

$\begin{array}[]{lll}v(x_{1},x_{2})&=&\displaystyle\Delta_{2}(t)\sum_{j=1}^{2}x_{j}^{2}+\sum_{\rho=1}^{2-1}\sum_{j=1}^{2}\prod_{s=1,s\not=\rho+1}^{2}\Delta_{s}(t)\Delta_{\rho,j}^{2}(x_{1},x_{2})\\\
&=&\displaystyle\Delta_{2}(t)\sum_{j=1}^{2}x_{j}^{2}+\sum_{j=1}^{2}\Delta_{1}(t)\Delta_{1,j}^{2}(x_{1},x_{2})\\\
&=&\displaystyle\Delta_{2}(t)\sum_{j=1}^{2}x_{j}^{2}+\Delta_{1}(t)[(a_{22}x_{1}-a_{12}x_{2})^{2}+(a_{21}x_{1}-a_{11}x_{2})^{2}]\\\
&=&\displaystyle\Delta_{2}(t)\sum_{j=1}^{2}x_{j}^{2}+\Delta_{1}(t)[(\Gamma_{1}(t)x_{1}+bx_{2})^{2}+(\gamma
x_{1}+|\gamma_{2}(t)x_{2})^{2}]\end{array}$ (26)

Mark $\Delta_{\rho,s}(x_{1},x_{2})$ can be found in reference [40], obviously

$v(x_{1},x_{2})\geqslant\Delta_{2}(t)(x_{1}^{2}+x_{2}^{2})=-\lambda_{1}\lambda_{2}(\lambda_{1}+\lambda_{2})(x_{1}^{2}+x_{2}^{2})\geqslant\delta^{3}(x_{1}^{2}+x_{2}^{2})$
(27)

Therefore, the above function $v(x_{1},x_{2})$ is positive definite, so we can
take the positive definite function $v(x_{1},x_{2})$ as the Lyapunov function
of the system (12).

$\left(\frac{\mathrm{d}v}{\mathrm{d}t}\right)_{\eqref{eq:linear
system}}=-2\Delta_{1}(t)\Delta_{2}(t)\sum_{j=1}^{2}x_{j}^{2}+\sum_{i,j=1}^{2}\dot{V}_{ij}(t)x_{i}x_{j}$
(28)

where

$\sum_{i,j=1}^{2}\dot{V}_{ij}(t)x_{i}x_{j}\leqslant\sum_{i,j=1}^{2}|\dot{V}_{ij}(t)|x_{i}x_{j}\leqslant\frac{1}{2}\sum_{i,j=1}^{2}|\dot{V}_{ij}(t)|(x_{i}^{2}+x_{j}^{2})=\sum_{i=1}^{2}\sum_{i=1}^{2}|\dot{V}_{ij}(t)|x_{i}^{2}$
(29)

Because $V_{ij}(t)=C(t)v_{ij}(t)$, so

$\dot{V}_{ij}(t)=\dot{C}(t)v_{ij}(t)+C(t)\dot{v}_{ij}(t)=\Big{(}\sum_{i,j=1}^{2}\frac{\partial
C}{\partial a_{ij}}\dot{a}_{ij}\Big{)}v_{ij}(t)+C(t)\dot{v}_{ij}(t)$ (30)

Let $|\dot{a}_{ij}(t)|\leqslant\varepsilon$, then

$|\dot{V}_{ij}(t)|\leqslant\varepsilon\Big{(}\sum_{i,j=1}^{2}\frac{\partial
C}{\partial a_{ij}}\Big{)}|v_{ij}(t)|+\varepsilon|C(t)|P_{ij}$ (31)

where $P_{ij}$ is the sum of the absolute values of the algebraic cofactors of
the elements whose derivative is not zero in the determinant of $V_{ij}$ (if
an element appears in the form of $a_{ii}+a_{jj}$, multiply by 2 before the
absolute value of its algebraic cofactors), which is obtained from (31)

$|\dot{V}_{ij}(t)|\leqslant\varepsilon
D(t)|v_{ij}(t)|+\varepsilon|C(t)|P_{ij}(t).$ (32)

Among them $\sum_{i,j=1}^{2}\frac{\partial C}{\partial a_{ij}}=D(t)$, so

$\left(\frac{\mathrm{d}v}{\mathrm{d}t}\right)_{\eqref{eq:linear
system}}\leqslant-2\Delta_{1}(t)\Delta_{2}(t)\sum_{j=1}^{2}x_{j}^{2}+\varepsilon\Big{\\{}\Big{[}D(t)Q_{1}(t)+|C(t)|P_{1}(t)\Big{]}x_{1}^{2}+D(t)Q_{2}(t)+|C(t)|P_{2}(t)\Big{]}x_{2}^{2}\Big{\\}},$
(33)

where

$Q_{1}(t)=\sum_{j=1}^{2}|V_{1j}(t)|,Q_{2}(t)=\sum_{j=1}^{2}|V_{2j}(t)|,P_{1}(t)=\sum_{j=1}^{2}P_{1j}(t),P_{2}(t)=\sum_{j=1}^{2}P_{2j}(t),$
(34)

Take

$\varepsilon=\min\Big{\\{}\frac{\Delta_{1}(t)\Delta_{2}(t)}{D(t)Q_{1}+|C(t)|P_{1}},\frac{\Delta_{1}(t)\Delta_{2}(t)}{D(t)Q_{2}+|C(t)|P_{2}}\Big{\\}}>0$
(35)

so when $|\dot{a}_{ij}(t)|\leqslant\varepsilon$, there are

$\left(\frac{\mathrm{d}v}{\mathrm{d}t}\right)_{\eqref{eq:linear
system}}\leqslant-2\Delta_{1}(t)\Delta_{2}(t)\sum_{j=1}^{2}x_{j}^{2}$ (36)

By using the condition $\mbox{Re}\lambda(A(t))\leqslant-\delta<0$

$\left(\frac{\mathrm{d}v}{\mathrm{d}t}\right)_{\eqref{eq:linear
system}}\leqslant-K^{*}\sum_{j=1}^{2}x_{j}^{2}$ (37)

where $K^{*}$ is a constant independent of $t$, so
$\left(\frac{\mathrm{d}v}{\mathrm{d}t}\right)_{\eqref{eq:linear system}}$ is
negative definite.

Since $V(x_{1},x_{2})$) is positive definite and has an infinitesimal upper
limit (obtained immediately from $|a_{ij}(t)|\leqslant a$), combined with the
conclusion proved previously, the zero solution of the system (12) is
asymptotically stable. Because the system (12) is linear, the stability has
global properties. ∎

##  3 Survival analysis of stochastic single population model

###  3.1 Stochastic Model

Because the contact rate between a single population and environmental toxins
is inevitably affected by other conditions, such as weather conditions,
temperature and other nearby noises, it is natural to think that constant
contact rate $\lambda_{i}$ is replaced by random variable
$\lambda_{i}+\sigma_{i}\xi_{i}(t)$, where $\xi_{i}(t)$ is white noise, and
$\xi_{i}(t)=\frac{dB_{i}(t)}{dt}$, $B_{i}(t)$ are two independent standard
Brown motions defined on $(\Omega,\mathcal{F},P)$. Therefore, we get random
single species with psychological effects group model

$\left\\{\begin{array}[]{lll}\vspace{1ex}\displaystyle\mathrm{d}x(t)=[by(t)-dx(t)-\gamma
x(t)-c_{1}x^{2}(t)-\alpha_{1}x(t)C_{o}(t)-\lambda_{1}x(t)C_{e}(t)]\mathrm{d}t-\sigma_{1}x(t)C_{e}(t)\mathrm{d}B_{1}(t),\\\
\vspace{1ex}\displaystyle\mathrm{d}y(t)=\Big{[}\gamma
x(t)-c_{2}y^{2}(t)-\alpha_{2}y(t)C_{o}(t)-\frac{\lambda_{2}y(t)C_{e}^{p}(t)}{1+\beta
C_{e}^{q}(t)}\Big{]}-\frac{\sigma_{2}y(t)C_{e}^{p}(t)}{1+\beta
C_{e}^{q}(t)}\mathrm{d}B_{2}(t),\\\
\vspace{1ex}\displaystyle\mathrm{d}C_{o}(t)=[kC_{e}(t)-gC_{o}(t)-mC_{o}(t)-bC_{o}(t)]\mathrm{d}t,\\\\[8.53581pt]
\displaystyle\mathrm{d}C_{e}(t)=[u_{e}(t)-hC_{e}(t)]\mathrm{d}t.\end{array}\right.$
(38)

For the convenience of subsequent derivation, we propose the following lemmas:

######  Lemma 2.

Under the condition of non-pollution, the system (38) has a unique equilibrium
point $(A,B)$ in the first quadrant.

###### Proof.

Under the condition of non-pollution, the system (38) is simplified as

$\left\\{\begin{array}[]{l}\dot{x}(t)=by(t)-dx(t)-\gamma
x(t)-c_{1}x^{2}(t),\\\ \dot{y}(t)=\gamma
x(t)-c_{2}y^{2}(t),\end{array}\right.$ (39)

The equilibrium point of the system is required, that is, to find the solution
of the following equations

$\left\\{\begin{array}[]{l}0=by(t)-dx(t)-\gamma x(t)-c_{1}x^{2}(t),\\\
0=\gamma x(t)-c_{2}y^{2}(t),\end{array}\right.$ (40)

Since the parameters of the system are all positive, there is only one
solution in the first quadrant of the equation group, which is set as $(A,B)$.
For the convenience of writing, we write the equation group (40) in the
following form:

$\left\\{\begin{array}[]{l}x=\tilde{a}y^{2},\\\
y=\tilde{b}x+\tilde{c}x^{2},\end{array}\right.$ (41)

of which

$\tilde{a}=\frac{c_{2}}{\gamma},\tilde{b}=\frac{d+\gamma}{b},\tilde{c}=\frac{c_{1}}{b}.$
(42)

The solution of equations (41) in the first quadrant is $(A,B)$, where

$A=\tilde{a}B^{2},B=\Lambda-\frac{\tilde{b}}{3\tilde{a}\tilde{c}},\Lambda=\Big{(}\frac{\sqrt{27\tilde{c}+4\tilde{a}\tilde{b}^{3}}}{2\sqrt{27}\tilde{a}^{2}\tilde{c}^{\frac{3}{2}}}+\frac{1}{2\tilde{a}^{2}\tilde{c}}\Big{)}^{\frac{1}{3}}$
(43)

∎

######  Lemma 3.

For models (8) and (38), if $k\leqslant g+m+b$, $u_{e}^{*}\leqslant h$, then,
$0\leqslant C_{o}(t)<1,0\leqslant C_{e}(t)<1$ for $t\in\mathbb{R}_{+}$.

###### Proof.

The explicit solutions $C_{o}(t),C_{e}(t)$ are as follows:

$\begin{array}[]{ll}\vspace{1ex}\displaystyle
C_{o}(t)=k\int_{0}^{t}C_{e}(s)e^{-(g+m+b)(t-s)}\mathrm{d}s+C_{o}(0)e^{-(g+m+b)t},\\\
\vspace{1ex}\displaystyle
C_{e}(t)=\int_{0}^{t}u_{e}(s)e^{-h(t-s)}\mathrm{d}s+C_{e}(0)e^{-ht}.\end{array}$

for the given initial conditions $0\leqslant C_{o}(0)<1,0\leqslant
C_{e}(0)<1$, the explicit solutions can be estimated

$\begin{array}[]{ll}\vspace{1ex}\displaystyle C_{o}(t)\leqslant
1-e^{-(g+m+b)t}+C_{o}(0)e^{-(g+m+b)t}<1,\\\ \vspace{1ex}\displaystyle
C_{e}(t)\leqslant 1-e^{-ht}+C_{e}(0)e^{-ht}<1.\end{array}$

So $0\leqslant C_{o}(t)<1$ and $0\leqslant C_{e}(t)<1$ hold. ∎

###  3.2 Existence and uniqueness of a global positive solution

######  Theorem 4.

For any initial value
$\left(x(0),y(0),C_{o}(0),C_{e}(0)\right)\in\mathbb{R}_{+}^{4}$, there is a
unique solution of model (38), and the solution does not leave
$\mathbb{R}_{+}^{4}$ according to probability 1.

###### Proof.

Our proof refers to the work of Mao [41]. By Lyapunov method mentioned in
Lemma 2 of [42], and the references [43, 44], the coefficients of model (38)
are locally Lipschitz continuous. On $\left[0,\tau_{e}\right)$, for any
initial value $\left(x(0),y(0),C_{o}(0),C_{e}(0)\right)\in\mathbb{R}_{+}^{4}$,
there is a local solution $\left(x(t),y(t),C_{o}(t),C_{e}(t)\right)$, where
$\tau_{e}$ is the blow up time. In order to prove that the solution is global,
we need to prove that $\tau_{e}=\infty$ holds almost surely. Let $m_{0}>1$ be
large enough, for $m>m_{0}$, each component of the initial value
$(x(0),y(0),C_{0}(0),C_{e}(0))\in\mathbb{R}_{+}^{4}$ is in the interval
$[\frac{1}{m_{0}},m_{0}]$. Define the stop time

$\begin{array}[]{lll}\vspace{1ex}\tau_{m}&=&\inf\left\\{t\in[0,\tau_{e}):\min\left\\{x(t),y(t),C_{o}(t),C_{e}(t)\right\\}\leqslant\displaystyle\frac{1}{m}\right.\\\
\vspace{1ex}&&\left.\mbox{or
}\max\left\\{x(t),y(t),C_{o}(t),C_{e}(t)\right\\}\geqslant
m\right\\}.\end{array}$ (44)

In this paper, we set $\inf\phi=\infty$. Obviously, when $m\rightarrow\infty$,
$\tau_{m}$ increases and $\tau_{m}<\tau_{e}$. We denote
$\tau_{\infty}=\lim_{m\rightarrow\infty}\tau_{m}$. We claim that for all
$t\geqslant 0$,
$\left(x(t),y(t),C_{o}(t),C_{e}(t)\right)\in\mathbb{R}_{+}^{4}$, and
$\tau_{\infty}=\infty$. Use the counter evidence method below. Otherwise, for
any $m>m_{0}$, there is a pair of constants $T>0$, $\varepsilon\in(0,1)$ such
that $\mathbb{P}\left\\{\tau_{m}\leqslant T\right\\}\geqslant\varepsilon$. We
define a $C^{2}$-function as follows:
$V_{1}:\mathbb{R}_{+}^{4}\rightarrow\mathbb{R}_{+}$, and let $x=x(t)$,
$y=y(t)$, $C_{o}=C_{o}(t)$, $C_{e}=C_{e}(t)$.

$V_{1}(t)=(x-1-\ln x)+(y-1-\ln y)+(C_{o}-1-\ln C_{o})+(C_{e}-1-\ln C_{e}).$

For all $t\in\left[0,\tau_{e}\right)$, according to Itô’s formula, together
with $V_{1}(t)\coloneqq V_{1}\left(x(t),y(t),C_{o}(t),C_{e}(t)\right)$

$\begin{array}[]{lll}\mathrm{d}V_{1}(t)&=&\displaystyle\left(1-\frac{1}{x}\right)\mathrm{d}x+\left(1-\frac{1}{y}\right)\mathrm{d}y+\left(1-\frac{1}{C_{o}}\right)\mathrm{d}C_{o}+\left(1-\frac{1}{C_{e}}\right)\mathrm{d}C_{e}\\\\[12.0pt]
&&+\displaystyle\frac{1}{2x^{2}}(\mathrm{d}x)^{2}+\frac{1}{2y^{2}}(\mathrm{d}y)^{2}\\\\[12.0pt]
&=&\displaystyle\mathcal{L}V_{1}\mathrm{d}t-(x-1)\sigma_{1}C_{e}\mathrm{d}B_{1}(t)-(y-1)\frac{\sigma_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\mathrm{d}B_{2}(t)\end{array}$ (45)

where

$\begin{array}[]{lll}\mathcal{L}V_{1}(t)&=&(x-1)\left(b\frac{y}{x}-d-\gamma-
c_{1}x-\alpha_{1}C_{0}(t)-\lambda_{1}C_{e}\right)+\frac{1}{2}\sigma_{1}^{2}C_{e}^{2}\\\\[12.0pt]
&&\displaystyle+(y-1)\left(\gamma\frac{x}{y}-c_{2}y-\alpha_{2}C_{o}-\frac{\lambda_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\right)+\frac{\sigma_{2}C_{e}^{2p}}{2\left(1+\beta
C_{e}^{q}\right)^{2}}\\\\[12.0pt]
&&\displaystyle+kC_{e}-(g+m+b)C_{o}-k\frac{C_{e}}{C_{o}}+g+m+b+u_{e}-hC_{e}-\frac{u_{e}}{C_{e}}+h\\\\[12.0pt]
&<&\displaystyle\left(c_{1}-d\right)x-c_{1}x^{2}+\left(b+c_{2}\right)y-c_{2}y^{2}+\left(\alpha_{1}+\alpha_{2}\right)+\left(\lambda_{1}+\lambda_{2}\right)+\frac{\sigma_{1}^{2}+\sigma_{2}^{2}}{2}\\\\[12.0pt]
&&+\gamma+k+g+m+b+d+\left(u_{e}\right)^{*}+h\\\\[12.0pt]
&\leqslant&\displaystyle\frac{\left(c_{1}-d\right)^{2}}{4c_{1}}+\frac{\left(b+c_{2}\right)^{2}}{4c_{2}}+\left(\alpha_{1}+\alpha_{2}\right)+\left(\lambda_{1}+\lambda_{2}\right)+\frac{\sigma_{1}^{2}+\sigma_{2}^{2}}{2}\\\\[12.0pt]
&&+\gamma+k+g+m+b+d+\left(u_{e}\right)^{*}+h\\\\[12.0pt]
&\coloneqq&M>0.\end{array}$ (46)

Hence it may be concluded that

$\mathrm{d}V_{1}(t)<M\mathrm{d}t-(x-1)\sigma_{1}C_{e}\mathrm{d}B_{1}(t)-(y-1)\frac{\sigma_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\mathrm{d}B_{2}(t).$ (47)

Integrate the two sides of (47) from 0 to $\tau_{m}\wedge T$, and then take
the expectation to obtain the following inequality

$\begin{array}[]{lll}\displaystyle\mathbb{E}\left(V_{1}(\tau_{m}\wedge
T\right)&<{}&\displaystyle V_{1}(0)+\mathbb{E}\int_{0}^{\tau_{m}\wedge
T}M\mathrm{d}s\\\\[12.0pt] &&\displaystyle-\mathbb{E}\int_{0}^{\tau_{m}\wedge
T}(x(s)-1)\sigma_{1}C_{e}(s)\mathrm{d}B_{1}(s)\\\\[12.0pt]
&&\displaystyle-\mathbb{E}\int_{0}^{\tau_{m}\wedge
T}(y(s)-1)\frac{\sigma_{2}C_{e}^{p}(s)}{1+\beta
C_{e}^{q}(s)}\mathrm{d}B_{2}(s)\\\\[12.0pt]
&\leqslant{}&V_{1}(0)+MT.\end{array}$ (48)

Note that $\Omega_{m}=\left\\{\tau_{m}\wedge T\right\\}$, then
$\mathbb{P}\left\\{\Omega_{m}\right\\}\geqslant\varepsilon$. For any
$w\in\Omega_{m}$, by the definition of stop time, there is

$x\left(\tau_{m},\omega\right)\wedge y\left(\tau_{m},\omega\right)\wedge
C_{o}\left(\tau_{m},\omega\right)\wedge
C_{e}\left(\tau_{m},\omega\right)=\frac{1}{m},$ (49)

or

$x\left(\tau_{m},\omega\right)\vee y\left(\tau_{m},\omega\right)\vee
C_{o}\left(\tau_{m},\omega\right)\vee C_{e}\left(\tau_{m},\omega\right)=m.$
(50)

Therefore, there are

$V_{1}\left(\tau_{m},\omega\right)\geqslant\min\left\\{m-1-\ln
m,\frac{1}{m}-1+\ln m\right\\}$ (51)

Then inequality (48) is equivalent to

$\begin{array}[]{lll}V_{1}(0)+MT&\geqslant{}&\mathbb{E}\left[1_{\Omega_{m}}(\omega)\cdot
V_{1}\left(\tau_{m}\wedge T\right)\right]\\\\[12.0pt]
&\geqslant{}&\mathbb{P}\left\\{\Omega_{m}\right\\}\cdot\min\left\\{m-1-\ln
m,\frac{1}{m}-1+\ln m\right\\}\\\\[12.0pt]
&\geqslant{}&\varepsilon\min\left\\{m-1-\ln m,\frac{1}{m}-1+\ln
m\right\\}.\end{array}$ (52)

where $1_{\Omega_{m}}(\omega)$ is the demonstrative function of $\Omega_{m}$.
Let $m\rightarrow\infty$, we can get

$+\infty>V_{1}(0)+MT\geqslant+\infty.$ (53)

The contradiction arises, therefore,
$\mathbb{P}\left\\{\tau_{m}=\infty\right\\}=1$. ∎

###  3.3 Weakly persistent in the mean

######  Theorem 5.

Let $X(0)=\left(x(0),y(0),C_{o}(0),C_{e}(0)\right)$ is the initial value,
$X(t)$ be a solution of model (38). If the coefficients of (38) satisfy

$\begin{array}[]{lll}\vspace{1ex}\displaystyle\sigma_{1}^{2}<\frac{bB}{A}+c_{1}A-\frac{b}{2}-\frac{\gamma}{2},&\displaystyle\sigma_{2}^{2}<\frac{\gamma
A}{B}+c_{2}B-\frac{b}{2}-\frac{\gamma}{2},\\\ \vspace{1ex}\displaystyle
k<2(g+m+b),&\displaystyle\sigma_{1}^{2}A^{2}+\sigma_{2}^{2}B^{2}<2h-k-1.\end{array}$
(54)

Then the solution $X(t)$ of (38) has the properties:

$\limsup_{t\rightarrow\infty}\frac{1}{t}\mathbb{E}\int_{0}^{t}\left(\eta_{1}(x(s)-A)^{2}+\eta_{2}(y(s)-B)^{2}+\eta_{3}C_{o}^{2}(s)+\eta_{4}C_{e}^{2}(s)\right)\mathrm{d}s\leqslant
K.$ (55)

In other words, in the long run, the single population near the pollution-free
equilibrium point $(A,B,0,0)$ is weak persistent in the mean.

###### Proof.

We set new variables

$u_{1}(t)=x(t)-A,\quad u_{2}(t)=y(t)-B,\quad v(t)=C_{o}(t),\quad
w(t)=C_{e}(t),$ (56)

where the expressions of $A$ and $B$ are given in Lemma 2. Therefore, the
model (38) is rewritten as follows:

$\left\\{\begin{array}[]{l}\vspace{1ex}\displaystyle\mathrm{d}u_{1}=\left(u_{1}+A\right)\left(\frac{b\left(u_{2}+B\right)}{u_{1}+A}-\frac{bB}{A}-c_{1}u_{1}-\alpha_{1}v-\lambda_{1}w\right)\mathrm{d}t-\left(u_{1}+A\right)\sigma_{1}w\mathrm{d}B_{1}(t),\\\
\vspace{1ex}\displaystyle\mathrm{d}u_{2}=\left(u_{2}+B\right)\left(\frac{\gamma\left(u_{1}+A\right)}{u_{2}+B}-\frac{\gamma
A}{B}-c_{2}u_{2}-\alpha_{2}v-\frac{\lambda_{2}w^{p}}{1+\beta
w^{q}}\right)\mathrm{d}t-\left(u_{2}+B\right)\frac{\sigma_{2}w^{p}}{1+\beta
w^{q}}\mathrm{d}B_{2}(t),\\\
\vspace{1ex}\displaystyle\mathrm{d}v=[ku-(g+m+b)v]\mathrm{d}t,\\\
\vspace{1ex}\displaystyle\mathrm{d}w=\left[u_{e}-hw\right]\mathrm{d}t.\end{array}\right.$
(57)

We define a $C^{2}$-function $V_{2}=u_{1}^{2}+u_{2}^{2}+v^{2}+w^{2}.$ From
Itô’s formula

$\mathrm{d}V_{2}=\mathcal{L}V_{2}\mathrm{d}t-2\left(u_{1}+A\right)\sigma_{1}u_{1}w\mathrm{d}B_{1}(t)-2\left(u_{2}+B\right)\frac{\sigma_{2}u_{2}w^{p}}{1+\beta
w^{q}}\mathrm{d}B_{2}(t),$ (58)

by the elementary equality $2xy\leqslant x^{2}+y^{2}$ for $x>0$ and $y>0$, we
get

$\begin{array}[]{lll}\mathcal{L}V_{2}&=&2u_{1}\left(u_{1}+A\right)\left(\frac{b\left(u_{2}+B\right)}{u_{1}+A}-\frac{bB}{A}-c_{1}u_{1}-\alpha_{1}v-\lambda_{1}w\right)+\left(u_{1}+A\right)^{2}\sigma_{1}^{2}w^{2}\\\\[12.0pt]
&&+2u_{2}\left(u_{2}+B\right)\left(\frac{\gamma\left(u_{1}+A\right)}{u_{2}+B}-\frac{\gamma
A}{B}-c_{2}u_{2}-\alpha_{2}v-\frac{\lambda_{2}w^{p}}{1+\beta
w^{q}}\right)\\\\[12.0pt]
&&+\left(u_{2}+B\right)^{2}\frac{\sigma_{2}^{2}w^{2p}}{\left(1+\beta
w^{q}\right)^{2}}+2v[kw-(g+m+b)v]+2w\left(u_{e}-hw\right)\\\\[12.0pt]
&<&2bu_{1}u_{2}-\frac{2bB}{A}u_{1}^{2}+2\gamma u_{1}u_{2}-\frac{2\gamma
A}{B}u_{2}^{2}-2c_{1}Au_{1}^{2}+\sigma_{1}^{2}u_{1}^{2}\\\\[12.0pt]
&&+2\sigma_{1}^{2}Au_{1}w^{2}+\sigma_{1}^{2}A^{2}w^{2}-2c_{2}Bu_{2}^{2}+\sigma_{2}^{2}u_{2}^{2}+2\sigma_{2}^{2}Bu_{2}w^{2}+\sigma_{2}^{2}B^{2}w^{2}\\\\[12.0pt]
&&+kw^{2}+kv^{2}-2(g+m+b)v^{2}+w^{2}+u_{e}^{2}-2hw^{2}\\\\[12.0pt]
&\leqslant&-\eta_{1}u_{1}^{2}-\eta_{2}u_{2}^{2}-\eta_{3}v^{2}-\eta_{4}w^{2}+K,\end{array}$
(59)

where

$\begin{array}[]{lll}\eta_{1}=\frac{2bB}{A}+2c_{1}A-2\sigma_{1}^{2}-b-\gamma,&\eta_{2}=\frac{2\gamma
A}{B}+2c_{2}B-2\sigma_{2}^{2}-b-\gamma,\\\\[8.53581pt]
\eta_{3}=2(g+m+b)-k,&\eta_{4}=2h-\sigma_{1}^{2}A^{2}-\sigma_{1}^{2}B^{2}-k-1,\\\\[8.53581pt]
K=\sigma_{1}^{2}A^{2}+\sigma_{2}^{2}B^{2}+\left(u_{e}^{*}\right)^{2}.\end{array}$
(60)

Integrating from $0$ to $t$ on both sides of (58) and taking the expectation

$\mathbb{E}\int_{0}^{t}\left(\eta_{1}u_{1}^{2}(s)+\eta_{2}u_{2}^{2}(s)+\eta_{3}v^{2}(s)+\eta_{4}w^{2}(s)\right)\mathrm{d}s<\mathbb{E}(V_{2}(0)-V_{2}(t))+Kt\leqslant
Kt.$ (61)

Further, we obtain

$\limsup_{t\rightarrow\infty}\frac{1}{t}\mathbb{E}\int_{0}^{t}\left(\eta_{1}u_{1}^{2}(s)+\eta_{2}u_{2}^{2}(s)+\eta_{3}v^{2}(s)+\eta_{4}w^{2}(s)\right)\mathrm{d}s\leqslant
K.$ (62)

∎

###  3.4 Stochastic permanence

######  Theorem 6.

If the coefficients of the model satisfy

$\limsup_{t\rightarrow\infty}\left(\left(\alpha_{1}+\alpha_{2}\right)C_{o}(t)+\lambda_{1}C_{e}(t)+\frac{\lambda_{2}C_{e}^{p}(t)}{1+\beta
C_{e}^{q}(t)}+3\sigma^{2}C_{e}^{2}(t)\right)\leqslant
4\sqrt{b\gamma}-\frac{\eta}{\theta}-\frac{c_{1}^{2}+c_{2}^{2}}{2}-d-\gamma,$
(63)

where $\sigma=\max\\{\sigma_{1},\sigma_{2}\\}$, then the total number of the
juveniles and the adults is stochastic permanence.

###### Proof.

We choose positive constants $\theta$ and $\eta$ such that $\theta\geqslant
1$, $\eta\geqslant 1$, and define $C^{2}$-function
$V_{3}:\mathbb{R}_{+}^{4}\rightarrow\mathbb{R}_{+}$ as follows:

$V_{3}=e^{\eta t}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta}.$
(64)

By Itô’s formula, we obtain that

$\begin{array}[]{lll}\vspace{1ex}\mathrm{d}V_{3}&=&\displaystyle\eta e^{\eta
t}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta}\mathrm{d}t-2\theta
e^{\eta
t}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-1}\left(\frac{1}{x^{3}}\mathrm{d}x+\frac{1}{y^{3}}\mathrm{d}y\right)\\\
\vspace{1ex}&&\displaystyle+\theta e^{\eta
t}\left[\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-1}\frac{3}{x^{4}}+(\theta-1)\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-2}\frac{2}{x^{6}}\right](\mathrm{d}x)^{2}\\\
\vspace{1ex}&&\displaystyle+\theta e^{\eta
t}\left[\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-1}\frac{3}{y^{4}}+(\theta-1)\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-2}\frac{2}{y^{6}}\right](\mathrm{d}y)^{2}\\\
\vspace{1ex}&=&\displaystyle\mathcal{L}V_{3}\mathrm{d}t+M(t),\end{array}$ (65)

where

$\begin{array}[]{lll}M(t)&=&\displaystyle 2\theta
e^{\eta\,t}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-1}\left(\frac{\sigma_{1}C_{e}}{x^{2}}\mathrm{d}B_{1}(t)+\frac{\sigma_{1}C_{e}^{p}}{y^{2}\left(1+\beta
C_{e}^{q}\right)}\mathrm{d}B_{2}(t)\right),\\\\[12.0pt]
\mathcal{L}V_{3}&=&\displaystyle\eta e^{\eta
t}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta}\\\\[12.0pt]
&&\displaystyle-2\theta e^{\eta
t}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-1}\left[\frac{1}{x^{2}}\left(b\frac{y}{x}-d-\gamma-
c_{1}x-\alpha_{1}C_{o}-\lambda_{1}C_{e}\right)\right.\\\\[12.0pt]
&&\displaystyle\left.+\frac{1}{y^{2}}\left(\gamma\frac{x}{y}-c_{2}y-\alpha_{2}C_{o}-\frac{\lambda_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\right)\right]\\\\[12.0pt] &&\displaystyle+\theta e^{\eta
t}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-1}\left[\sigma_{1}^{2}C_{e}^{2}\left(\frac{3}{x^{2}}+\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{-1}\frac{2(\theta-1)}{x^{4}}\right)\right.\\\\[12.0pt]
&&\displaystyle\left.+\frac{\sigma_{2}^{2}C_{e}^{2p}}{\left(1+\beta
C_{e}^{q}\right)^{2}}\left(\frac{3}{y^{2}}+\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{-1}\frac{2(\theta-1)}{y^{4}}\right)\right]\\\\[12.0pt]
&\coloneqq&\displaystyle\theta e^{\eta
t}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-2}(I_{31}+I_{32}+I_{33}).\end{array}$
(66)

After simplification, we derive

$\begin{array}[]{lll}I_{31}&=&\displaystyle\frac{\eta}{\theta}\left[1+2\left(\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)+\left(\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{2}\right],\end{array}$
(67)
$\begin{array}[]{lll}I_{32}&<&\displaystyle-2\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)\left[\frac{2\sqrt{b\gamma}}{xy}-\frac{1}{x^{2}}\left(d+\gamma+\alpha_{1}C_{o}+\lambda_{1}C_{e}\right)\right.\\\\[12.0pt]
&&\displaystyle\left.-\frac{1}{y^{2}}\left(\alpha_{2}C_{o}+\frac{\lambda_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\right)\right]+2\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)\left(\frac{c_{1}}{x}+\frac{c_{2}}{y}\right),\end{array}$
(68) $\begin{array}[]{lll}I_{33}&=&\displaystyle
3\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)\left(\frac{\sigma_{1}^{2}C_{e}^{2}}{x^{2}}+\frac{\sigma_{2}^{2}C_{e}^{2p}}{y^{2}(1+\beta
C_{e}^{q})^{2}}\right)+2(\theta-1)\left(\frac{\sigma_{1}^{2}C_{e}^{2}}{x^{4}}+\frac{\sigma_{2}^{2}C_{e}^{2p}}{y^{2}(1+\beta
C_{e}^{q})^{2}}\right)\\\\[12.0pt] &<&\displaystyle
3\sigma^{2}C_{e}^{2}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)\left(\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)+2(\theta-1)\sigma^{2}C_{e}^{2}\left(\frac{1}{x^{4}}+\frac{1}{y^{4}}\right),\end{array}$
(69)

therefore we have

$\begin{array}[]{lll}I_{31}+I_{32}+I_{33}&<&\displaystyle\frac{\eta}{\theta}+\frac{2c_{1}}{x}+\frac{2c_{1}}{x^{3}}+\frac{2c_{2}}{y}+\frac{2c_{2}}{y^{3}}\\\\[12.0pt]
&&\displaystyle+\left(\frac{2\eta}{\theta}+2\left(d+\gamma+\alpha_{1}C_{o}+\lambda_{1}C_{e}\right)+3\sigma^{2}C_{e}^{2}+1\right)\frac{1}{x^{2}}\\\\[12.0pt]
&&\displaystyle+\left(\frac{2\eta}{\theta}+2\left(\alpha_{2}C_{o}+\frac{\lambda_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\right)+3\sigma^{2}C_{e}^{2}+1\right)\frac{1}{y^{2}}\\\\[12.0pt]
&&\displaystyle+\left(\frac{\eta}{\theta}+2\left(d+\gamma+\alpha_{1}C_{o}+\lambda_{1}C_{e}\right)+(3+2(\theta-1))\sigma^{2}C_{e}^{2}\right)\frac{1}{x^{4}}\\\\[12.0pt]
&&\displaystyle+\left(\frac{\eta}{\theta}+2\left(\alpha_{2}C_{o}+\frac{\lambda_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\right)+(3+2(\theta-1))\sigma^{2}C_{e}^{2}\right)\frac{1}{y^{4}}\\\\[12.0pt]
&&\displaystyle+2\left[\frac{\eta}{\theta}+\left(d+\gamma+\alpha_{1}C_{o}+\lambda_{1}C_{e}\right)+\left(\alpha_{2}C_{o}+\frac{\lambda_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\right)+3\sigma^{2}C_{e}^{2}\right.\\\\[12.0pt]
&&\displaystyle\left.+\frac{c_{1}^{2}+c_{2}^{2}}{2}-4\sqrt{b\gamma}\right]\frac{1}{x^{2}y^{2}}-\frac{4\sqrt{b\gamma}}{xy}.\end{array}$
(70)

By condition (63) of Theorem 6, so we get

$\begin{array}[]{lll}I_{31}+I_{32}+I_{33}&<&\displaystyle\frac{\eta}{\theta}+\frac{2c_{1}}{x}+\frac{2c_{1}}{x^{3}}+\frac{2c_{2}}{y}+\frac{2c_{2}}{y^{3}}\\\\[12.0pt]
&&\displaystyle+\left(\frac{2\eta}{\theta}+2\left(d+\gamma+\alpha_{1}+\lambda_{1}\right)+3\sigma^{2}+1\right)\frac{1}{x^{2}}\\\\[12.0pt]
&&\displaystyle+\left(\frac{2\eta}{\theta}+2\left(\alpha_{2}+\lambda_{2}\right)+3\sigma^{2}+1\right)\frac{1}{y^{2}}\\\\[12.0pt]
&&\displaystyle+\left(\frac{\eta}{\theta}+2\left(d+\gamma+\alpha_{1}+\lambda_{1}\right)+(3+2(\theta-1)\sigma^{2}\right)\frac{1}{x^{4}}\\\\[12.0pt]
&&\displaystyle+\left(\frac{\eta}{\theta}+2\left(\alpha_{2}+\lambda_{2}\right)+(3+2(\theta-1)\sigma^{2}\right)\frac{1}{y^{4}}\\\\[12.0pt]
&\coloneqq&F(x,y).\end{array}$ (71)

It is obvious that

$H(x,y)\coloneqq
F(x,y)G(x,y),\quad\,G(x,y)=\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta-2}$
(72)

admits an upper boundary for all $x$ and $y$ in the domain $\mathbb{R}_{+}$.
Let us denote

$H_{1}=\sup_{(x,y)\in\mathbb{R}_{+}^{2}}H(x,y)<\infty.$ (73)

Integrating (65), together with $x_{0}=x(0),~{}y_{0}=y(0)$, which gives that

$e^{\eta
t}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta}<\left(1+\frac{1}{x_{0}^{2}}+\frac{1}{y_{0}^{2}}\right)^{\theta}+\frac{H_{1}\theta}{\eta}(e^{\eta
t}-1)+\int_{0}^{t}M(s)\mathrm{d}s,$ (74)

taking the expectation on both sides, then

$\mathbb{E}\left[\frac{1}{2^{\theta}}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta}\right]<\mathbb{E}\left[\frac{e^{-\eta\,t}}{2^{\theta}}\left(1+\frac{1}{x_{0}^{2}}+\frac{1}{y_{0}^{2}}\right)^{\theta}+\frac{H_{1}\theta}{2^{\theta}\eta}(1-e^{-\eta\,t})+\frac{e^{-\eta\,t}}{2^{\theta}}\int_{0}^{t}M(s)\mathrm{d}s\right].$
(75)

Therefore, together with (75), taking supremum limit, we derive that

$\limsup_{t\rightarrow\infty}\mathbb{E}\left(\frac{1}{x+y}\right)^{2\theta}\leqslant\limsup_{t\rightarrow\infty}\mathbb{E}\left(\frac{1}{x^{2}+y^{2}}\right)^{\theta}<\limsup_{t\rightarrow\infty}\left[\frac{1}{2^{\theta}}\left(1+\frac{1}{x^{2}}+\frac{1}{y^{2}}\right)^{\theta}\right]<\frac{H_{1}\theta}{2^{\theta}\eta}\coloneqq
H_{2}.$ (76)

For any given $\varepsilon>0$, we let

$\xi=\left(\frac{\varepsilon}{H_{2}}\right)^{\frac{1}{2\theta}},$

the Chebyshev inequality gives that

$\mathbb{P}\\{x+y<\xi\\}=\mathbb{P}\left\\{(x+y)^{-\theta}>\xi^{-\theta}\right\\}\leqslant\frac{\mathbb{E}\left((x+y)^{-2\theta}\right)}{\xi^{-2\theta}}=\xi^{2\theta}\mathbb{E}\left((x+y)^{-2\theta}\right).$

Combined with (76), we have

$\lim_{t\rightarrow\infty}\mathbb{P}\\{x+y<\xi\\}<\varepsilon,$ (77)

which implies that

$\liminf_{t\rightarrow\infty}\mathbb{P}\\{x+y\geqslant\xi\\}\geqslant
1-\varepsilon.$ (78)

We define a $C^{2}$-function
$V_{4}=e^{t}\left(x^{m}+y^{n}\right)^{l},0<l<1,m,n\in\mathbb{N}_{+}.$ By Itô’s
formula again

$\begin{array}[]{lll}\vspace{1ex}\mathrm{d}V_{4}&=&\displaystyle\left.e^{t}\left(x^{m}+y^{n}\right)\mathrm{d}t+le^{t}\left(x^{m}+y^{n}\right)\right)^{l-1}\left(mx^{m-1}\mathrm{d}x+ny^{n-1}\mathrm{d}y\right)\\\
\vspace{1ex}&&\displaystyle+\frac{1}{2}l(l-1)e^{t}\left(x^{m}+y^{n}\right)^{l-2}\left(m^{2}x^{2(m-1)}(\mathrm{d}x)^{2}+n^{2}y^{2(n-1)}(\mathrm{d}y)^{2}\right)\\\
\vspace{1ex}&&\displaystyle+\frac{1}{2}le^{t}\left(x^{m}+y^{n}\right)^{l-1}\left(m(m-1)x^{m-2}(\mathrm{d}x)^{2}+n(n-1)y^{n-2}(\mathrm{d}y)^{2}\right)\\\
\vspace{1ex}&=&\displaystyle\mathcal{L}V_{4}\mathrm{d}t-le^{t}\left(x^{m}+y^{n}\right)^{l-1}\left(m\sigma_{1}x^{m}C_{e}\mathrm{d}B_{1}(t)+\frac{n\sigma_{2}y^{n}C_{e}^{p}}{1+\beta
C_{e}^{q}}\mathrm{d}B_{2}(t)\right),\end{array}$

where

$\begin{array}[]{lll}\vspace{1ex}\mathcal{L}V_{4}&=&\displaystyle
e^{t}(x^{m}+y^{n})^{l}+le^{t}(x^{m}+y^{n})^{l-1}\Big{[}\Big{(}b\frac{y}{x}-d-\gamma-
c_{1}x-\alpha_{1}C_{o}-\lambda_{1}C_{e}\Big{)}mx^{m}\\\
\vspace{1ex}&&\displaystyle+\Big{(}\gamma\frac{x}{y}-c_{2}y-\alpha_{2}C_{o}-\frac{\lambda_{1}C_{e}^{p}}{1+\beta
C_{e}^{q}}\Big{)}ny^{n}\Big{]}\\\
\vspace{1ex}&&\displaystyle+\frac{1}{2}l(l-1)e^{t}(x^{m}+y^{n})^{l-2}\Big{(}m^{2}\sigma_{1}^{2}x^{2m}C_{e}^{2}+\frac{n^{2}\sigma_{2}^{2}y^{2n}C_{e}^{2p}}{(1+\beta
C_{e}^{q})^{2}}\Big{)}\\\
\vspace{1ex}&&\displaystyle+\frac{1}{2}le^{t}(x^{m}+y^{n})^{l-1}\Big{(}m(m-1)\sigma_{1}^{2}x^{m}C_{e}^{2}+\frac{n(n-1)\sigma_{2}^{2}y^{n}C_{e}^{2p}}{(1+\beta
C_{e}^{q})^{2}}\Big{)}\\\ \vspace{1ex}&<&\displaystyle
e^{t}(x^{m}+y^{n})^{l-1}\Big{\\{}x^{m}+y^{n}+l\Big{[}\Big{(}b\frac{y}{x}-d-\gamma-
c_{1}x-\alpha_{1}C_{o}-\lambda_{1}C_{e}\Big{)}mx^{m}\\\
\vspace{1ex}&&\displaystyle+\Big{(}\gamma\frac{x}{y}-c_{2}y-\alpha_{2}C_{o}-\frac{\lambda_{1}C_{e}^{p}}{1+\beta
C_{e}^{q}}\Big{)}ny^{n}\Big{]}\\\
\vspace{1ex}&&\displaystyle+\frac{1}{2}l\Big{(}m(m-1)\sigma_{1}^{2}x^{m}C_{e}^{2}+\frac{n(n-1)\sigma_{2}^{2}y^{n}C_{e}^{2p}}{(1+\beta
C_{e}^{q})^{2}}\Big{)}\Big{\\}}\\\ \vspace{1ex}&\leqslant&\displaystyle
e^{t}(x^{m}+y^{n})^{l-1}\Big{\\{}x^{m}+y^{n}\\\
\vspace{1ex}&&\displaystyle+\frac{bml}{2}(x^{2(m-1)}+y^{2})+\frac{ml}{2}(x^{2}+y^{2(n-1)})-c_{1}mlx^{m+1}-c_{2}nly^{n+1}\\\
\vspace{1ex}&&\displaystyle+l\Big{[}\Big{(}-d-\alpha_{1}C_{o}-\lambda_{1}C_{e}\Big{)}mx^{m}+\Big{(}-\alpha_{2}C_{o}-\frac{\lambda_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\Big{)}ny^{n}\Big{]}\\\
\vspace{1ex}&&\displaystyle+\frac{1}{2}l\Big{(}m(m-1)\sigma_{1}^{2}x^{m}C_{e}^{2}+\frac{n(n-1)\sigma_{2}^{2}y^{n}C_{e}^{2p}}{(1+\beta
C_{e}^{q})^{2}}\Big{)}\Big{\\}}\\\ \vspace{1ex}&=&\displaystyle
e^{t}(x^{m}+y^{n})^{l-1}\Big{\\{}-c_{1}mlx^{m+1}-c_{2}nly^{n+1}\\\
\vspace{1ex}&&\displaystyle+\Big{[}1+\frac{1}{2}m(m-1)l\sigma_{1}^{2}C_{e}^{2}-ml(d+\alpha_{1}C_{o}+\lambda_{1}C_{e})\Big{]}x^{m}\\\
\vspace{1ex}&&\displaystyle+\Big{[}1+\frac{n(n-1)l\sigma_{2}^{2}C_{e}^{2p}}{2(1+\beta
C_{e}^{q})^{2}}-nl\Big{(}\alpha_{2}C_{o}+\frac{\lambda_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\Big{)}\Big{]}y^{n}\\\
\vspace{1ex}&&\displaystyle+\frac{bml}{2}(x^{2(m-1)}+y^{2})+\frac{ml}{2}(x^{2}+y^{2(n-1)})\Big{\\}}.\end{array}$
(79)

We take $m=n=2$, then

$\begin{array}[]{lll}\vspace{1ex}\mathcal{L}V_{4}&<&e^{t}\left(x^{2}+y^{2}\right)^{l-1}\left\\{-2c_{1}lx^{3}-2c_{2}ly^{3}\right.\\\
\vspace{1ex}&&+\left[1+l\sigma_{1}^{2}C_{e}^{2}-2l\left(d+\alpha_{1}C_{o}+\lambda_{1}C_{e}\right)+bl+\gamma\right]x^{2}\\\
\vspace{1ex}&&\displaystyle\left.+\left[1+\frac{l\sigma_{2}^{2}C_{e}^{2p}}{\left(1+\beta
C_{e}^{q}\right)^{2}}-2l\left(\alpha_{2}C_{o}+\frac{\lambda_{2}C_{e}^{p}}{1+\beta
C_{e}^{q}}\right)+bl+\gamma\right]y^{2}\right\\}\\\
\vspace{1ex}&\coloneqq&e^{t}J(x,y)\leqslant Ce^{t}\end{array}$ (80)

taking the expectation on both sides of (80), then

$\mathbb{E}(e^{t}(x^{2}+y^{2})^{l})-(x^{2}(0)+y^{2}(0))^{l}<C(e^{t}-1),$ (81)

where $C$ is a constant. Due to $x+y\leqslant x^{2}+y^{2},$ we have

$\liminf_{t\rightarrow\infty}\mathbb{E}\left(x^{2}+y^{2}\right)<C.$ (82)

For any given $\varepsilon>0$, we let

$N=\left(\frac{C}{\varepsilon}\right)^{\frac{\,1}{\,l}},$ (83)

by Chebyshev inequality

$\mathbb{P}\\{x+y>N\\}=\mathbb{P}\\{(x+y)^{\frac{l}{2}}>N^{\frac{l}{2}}\\}\leqslant\frac{\mathbb{E}(x+y)^{l}}{N^{l}}.$
(84)

Combined with (82), we have

$\liminf_{t\rightarrow\infty}\mathbb{P}\\{x+y>N\\}<\varepsilon,$ (85)

which implies that

$\limsup_{t\rightarrow\infty}\mathbb{P}\\{x+y\leqslant N\\}\geqslant
1-\varepsilon.$ (86)

In other words, the single population is random and persistent. ∎

##  4 Numerical simulations

To support the main conclusions of this paper, we selected the appropriate
parameters and initial values for each conclusion. The main method of
numerical simulation is to discretize the original system and use matlab to
observe the trend of the graph. Then verify the conclusions of this paper. At
the same time of verifying random and persistent survival, under the premise
that the remaining parameters are fixed, we take 4 sets of different
psychological effect intensities for numerical simulation.

The parameters and their reference value are shown on this table:

Table 1: Parameters and their reference value Parameter | Value | Reference value | Source
---|---|---|---
$b$ | 0.3 | [0.1,0.3] | [29, 32, 45]
$d$ | 0.01 | 0.01 | [29]
$\gamma$ | 0.5 | [0.2,0.4] | [32, 45]
$c_{1}$ | 0.02 | - | Assumed
$c_{2}$ | 0.05 | [0.1,0.4] | [29, 32, 45]
$\alpha_{1}$ | 0.2 | - | Assumed
$\alpha_{2}$ | 0.1 | [0.1,0.5] | [29, 32, 45]
$\lambda_{1}$ | 0.07 | - | Assumed
$\lambda_{2}$ | 0.05 | 0.05 | [29, 32]
$k$ | 0.1 | [0.1,0.6] | [29, 46, 32, 45]
$g$ | 0.08 | [0.08,0.3] | [29, 46, 32, 45]
$m$ | 0.04 | [0.04,0.3] | [29, 46, 32, 45]
$h$ | 0.8 | [0.1,0.5] | [29, 46, 32, 47, 45]
$u_{e}$ | [0,0.6] | [0,0.8] | [29, 46, 45]
$\beta$ | [0.1,10] | [0.1,10] | [29, 32]
$\sigma_{1}$, $\sigma_{2}$ | [0.1,0.3] | [0.1,0.3] | [29]

###  4.1 Simulation of globally unique positive solution

The solution of equations (41) in the first quadrant is the globally unique
positive solution, noted as $(A,B)$. Take the value in Table 1, then we can
get $A=2.8078$ and $B=5.2989$. We draw and calibrate $(A,B)$ as shown in the
Figure 1.

Figure 1: The globally unique positive solution of equations (41)

###  4.2 Simulations of weakly persistent in the mean

Here we use numerical simulation results to support the conclusion of Theorem
3. We determine the parameters and initial values that meet the conditions of
the theorem. Respectively, we select two toxin input levels for numerical
simulation, through the image branch of numerical simulation to support the
conclusion of subsection 3.3.

Considering that theorem 3 requires non-pollution, according to the actual
situation, we take the smaller toxin input level $u_{e}\left(t\right)=0.1$ and
$u_{e}\left(t\right)=0.1+0.1\sin(t)$, and the psychological effect intensity
$\beta=0.1$, and then take the remaining parameters the same as Table 1 while
the variable parameters are shown on this Table 2:

Table 2: Variable parameters and their value Group | $u_{e}$ | $\beta$ | $\sigma_{1}$ | $\sigma_{2}$
---|---|---|---|---
Group 1 | 0.1 | 0.1 | 0.15 | 0.1
Group 2 | $0.1+0.1\sin(t)$ | 0.1 | 0.15 | 0.1

The diagrams 2 below track the density curves of $x(t)$ and $y(t)$ in Figure
2(a) and Figure 2(b) and toxin concentration $C_{o}(t)$ and $C_{e}(t)$ over
time in Figure 2(c) and Figure 2(d), when this population is weakly persistent
in the mean.

(a) The density curves of $x(t)$

(b) The density curves of $y(t)$

(c) Toxin concentration $C_{o}(t)$

(d) Toxin concentration $C_{e}(t)$

Figure 2: The density curves of $x(t)$ and $y(t)$, Toxin concentration
$C_{o}(t)$ and $C_{e}(t)$ with different $u_{e}$.

###  4.3 Simulations of stochastic permanence

We take the parameters and initial values that meet the conditions of the
Theorem 6, select the input levels of two toxins respectively, and group them
according to the different noise intensity. At the same time, in order to
verify the influence of the intensity of psychological effects on the dynamics
of a single population, we divide them under the premise that the other
parameters are fixed. We determine the variable parameters as in Table 3,
while others are consistent with Table 1.

The variable parameters are shown on this table:

Table 3: Variable parameters and their value Figure | $u_{e}$ | $\beta$ | $\sigma_{1},\sigma_{2}$
---|---|---|---
a | 0.3 | 0.1 | 0.1 or 0.3
b | 0.3 | 1 | 0.1 or 0.3
c | 0.3 | 10 | 0.1 or 0.3
d | $0.3+0.3\sin(t)$ | 0.1 | 0.1 or 0.3
e | $0.3+0.3\sin(t)$ | 1 | 0.1 or 0.3
f | $0.3+0.3\sin(t)$ | 10 | 0.1 or 0.3

The three diagrams in Figure 3, Figure 4 and Figure 5 respectively track the
density curves of $x(t)$ and $y(t)$, toxin concentration $C_{o}(t)$ and
$C_{e}(t)$ over time, when this population is stochastic permanence.

We first take $u_{e}(t)=0.3$, and on the premise of fixing the noise intensity
of the two groups, respectively take four groups of psychological effect
intensity $\beta=0.1,1,10$ for comparison. The results are shown in Figure
3(a)-3(c) and Figure 4(a)-4(c).

Then we take $u_{e}\left(t\right)=0.3+0.3\sin(t)$, on the premise of fixing
the noise intensity of the two groups, respectively take four groups of
psychological effect intensity $\beta=0.1,1,10$ for comparison, the results
are shown in Figure 3(d)-3(f) and Figure 4(d)-4(f):

The results show that when the psychological effect intensity increases, the
adult population size will tend to be stable, which means that adults will
actively evade the living area with high toxin concentration.

(a) $u_{e}=0.3,\beta=0.1$

(b) $u_{e}=0.3,\beta=1$

(c) $u_{e}=0.3,\beta=10$

(d) $u_{e}=0.3+0.3\sin(t),\beta=0.1$

(e) $u_{e}=0.3+0.3\sin(t),\beta=1$

(f) $u_{e}=0.3+0.3\sin(t),\beta=10$

Figure 3: The density curves of $x(t)$ with different $u_{e}$ and $\beta$.

(a) $u_{e}=0.3,\beta=0.1$

(b) $u_{e}=0.3,\beta=1$

(c) $u_{e}=0.3,\beta=10$

(d) $u_{e}=0.3+0.3\sin(t),\beta=0.1$

(e) $u_{e}=0.3+0.3\sin(t),\beta=1$

(f) $u_{e}=0.3+0.3\sin(t),\beta=10$

Figure 4: The density curves of $y(t)$ with different $u_{e}$ and $\beta$

(a) $u_{e}=0.3$

(b) $u_{e}=0.3+0.3\sin(t)$

Figure 5: Toxin concentration $C_{o}(t)$ and $C_{e}(t)$ with different
$u_{e}$.

##  5 Discussion

The results of this paper can be further extended, especially for the single
population model of such stage structure. There are many long-term behaviors
that can be studied, and different methods can be used for theoretical
practice. In reality, it can guide production and life in some fields, such as
fishery production. When the production environment is polluted, it can be
guided by relevant theories to take measures in time. At the same time, for
the biological aquaculture which needs to separate the larva and adult, the
subject can be developed into the theoretical basis for guiding the capture
strategy [10], and also can be introduced into the control theory.

There is still space for further development for the two-stage structure model
established in this paper. We assume that young individuals are transformed
into adult individuals at a certain rate. In reality, the maturity process of
any organism needs time accumulation, that is, the generation of time delay.
In this regard, the model established in this paper also has certain
limitations. Many scholars consider the existence of this factor, which leads
to the time-delay term is used to describe the mature process [16, 15, 17].

##  6 Conclusion

In general, based on the previous studies, we introduce two-stage age
structure, develop and expand some results. We develop the psychological
effect function, select different psychological effect functions for different
age structure, and extract some existing function types.

In subsection 2.1, a single population model with age structure and
psychological effects in polluted environment is established, which is a
nonlinear time-varying system. Then in subsection 2.2, we discuss the
asymptotic stability of the system by Lyapunov first approximation theory, and
give a sufficient condition for the stability.

In subsection 3.1, based on subsection 2.1, the exposure rate of organisms to
environmental toxins is affected by white noise, and then the exposure rate is
modified into a random process, and the corresponding random single population
model is established. The subsequent contents are mostly proved by Lyapunov
function method. In subsection 3.2, we verify the existence of the globally
unique positive solution of the stochastic model. In subsection 3.3, near the
non-pollution equilibrium point, we give the sufficient conditions for the
weak mean persistence of a single population in the expected sense. In
subsection 3.4, we give the sufficient conditions for the random persistence
of a single population.

In section 4, we make some numerical simulation to support our conclusions in
Theorem 5 and Theorem 6.

## Conflict of Interest

All authors stated that they have no conflict of interests in the paper.

## Acknowledgment

The research is supported by National Natural Science Foundation of China
(61911530398), Special Projects of the Central Government Guiding Local
Science and Technology Development (2021L3018), the Natural Science Foundation
of Fujian Province of China (2021J01621), and Scientific Research Training
Program in Fuzhou University (No.26040).

## References

  * [1] T. G. Hallam, C. E. Clark, R. R. Lassiter, Effects of toxicants on populations: a qualitative approach I. Equilibrium environmental exposure, Ecological Modelling 18 (3-4) (1983) 291–304.
  * [2] T. G. Hallam, C. E. Clark, G. S. Jordan, Effects of toxicants on populations: a qualitative approach II. First order kinetics, J Math Biol 18 (1) (1983) 25–37.
  * [3] T. G. Hallam, J. T. de Luna, Effects of toxicants on populations: A qualitative, Journal of Theoretical Biology 109 (3) (1984) 411–429.
  * [4] T. G. Hallam, M. Zhien, Persistence in population models with demographic fluctuations, J Math Biol 24 (3) (1986) 327–39.
  * [5] T. G. Hallam, M. Zhien, On density and extinction in continuous population models, J Math Biol 25 (2) (1987) 191–201.
  * [6] M. Zhien, S. Baojun, T. Hallam, The threshold of survival for systems in a fluctuating environment, Bulletin of Mathematical Biology 51 (3) (1989) 311–323.
  * [7] T. C. Gard, A stochastic model for the effects of toxicants on populations, Ecological Modelling 51 (3-4) (1990) 273–280.
  * [8] T. C. Gard, Stochastic models for toxicant-stressed populations, Bull Math Biol 54 (5) (1992) 827–37.
  * [9] D. M. Thomas, T. W. Snell, S. M. Jaffar, A control problem in a polluted environment, Math Biosci 133 (2) (1996) 139–63.
  * [10] H. Zerong, M. Zhien, On the effect of pollution and catch to a logistic population, Journal of Biomathematics 12 (3) (1997) 230–237.
  * [11] H. Ze-rong, M. Zhi-en, The survival analysis for gallopin’s system in a polluted environment, Applied Mathematics and Mechanics 21 (8) (2000) 921–928.
  * [12] P. Srinivasu, Control of environmental pollution to conserve a population, Nonlinear analysis: real world applications 3 (3) (2002) 397–411.
  * [13] H. D. Landahl, B. D. Hansen, A three stage population model with cannibalism, Bull Math Biol 37 (1) (1975) 11–7.
  * [14] K. Tognetti, The two stage stochastic population model, Mathematical Biosciences 25 (3-4) (1975) 195–204.
  * [15] S. N. Wood, S. P. Blythe, W. S. C. Gurney, R. M. Nisbet, Instability in mortality estimation schemes related to stage-structure population models, Mathematical Medicine and Biology 6 (1) (1989) 47–68.
  * [16] W. G. Aiello, H. I. Freedman, A time-delay model of single-species growth with stage structure, Math Biosci 101 (2) (1990) 139–53.
  * [17] W. G. Aiello, H. I. Freedman, J. Wu, Analysis of a model representing stage-structured population growth with state-dependent time delay, SIAM Journal on Applied Mathematics 52 (3) (1992) 855–869.
  * [18] Y.-g. Yu, S.-c. Zhang, Z.-z. Yang, The effect of dispersal on population growth with stage-structure, Acta Mathematicae Applicatae Sinica, English Series 19 (3) (2017) 499–504.
  * [19] R. Rudnicki, Long-time behaviour of a stochastic prey–predator model, Stochastic Processes and their Applications 108 (1) (2003) 93–107.
  * [20] M. Zhien, C. Guirong, W. Wendi, Persistence and extinction of a population in a polluted environment, Mathematical Biosciences 101 (1) (1990) 75–97.
  * [21] A. Bahar, X. Mao, Stochastic delay Lotka–Volterra model, Journal of Mathematical Analysis and Applications 292 (2) (2004) 364–380.
  * [22] S. Liu, E. Beretta, A stage-structured predator-prey model of beddington-deangelis type, SIAM Journal on Applied Mathematics 66 (4) (2006) 1101–1129.
  * [23] F. Wu, S. Hu, Stochastic functional Kolmogorov-type population dynamics, Journal of Mathematical Analysis and Applications 347 (2) (2008) 534–549.
  * [24] D. T. Lee, D. Sahota, T. N. Leung, A. S. Yip, F. F. Lee, T. K. Chung, Psychological responses of pregnant women to an infectious outbreak: a case-control study of the 2003 SARS outbreak in Hong Kong, J Psychosom Res 61 (5) (2006) 707–13.
  * [25] R. Liu, J. Wu, H. Zhu, Media/psychological impact on multiple outbreaks of emerging infectious diseases, Computational and Mathematical Methods in Medicine 8 (3) (2007) 153–164.
  * [26] D. Xiao, S. Ruan, Global analysis of an epidemic model with nonmonotone incidence rate, Mathematical Biosciences 208 (2) (2007) 419–429.
  * [27] A. K. Misra, A. Sharma, J. B. Shukla, Modeling and analysis of effects of awareness programs by media on the spread of infectious diseases, Mathematical and Computer Modelling 53 (5-6) (2011) 1221–1228.
  * [28] S. Samanta, S. Rana, A. Sharma, A. K. Misra, J. Chattopadhyay, Effect of awareness programs by media on the epidemic outbreaks: A mathematical model, Applied Mathematics and Computation 219 (12) (2013) 6965–6977.
  * [29] F. Wei, L. Chen, Psychological effect on single-species population models in a polluted environment, Math Biosci 290 (2017) 22–30.
  * [30] I. Ghosh, P. K. Tiwari, S. Samanta, I. M. Elmojtaba, N. Al-Salti, J. Chattopadhyay, A simple SI-type model for HIV/AIDS with media and self-imposed psychological fear, Math Biosci 306 (2018) 160–169.
  * [31] S. K. Sasmal, Population dynamics with multiple Allee effects induced by fear factors – A mathematical study on prey-predator interactions, Applied Mathematical Modelling 64 (2018) 1–14.
  * [32] G. Lan, C. Wei, S. Zhang, Long time behaviors of single-species population models with psychological effect and impulsive toxicant in polluted environments, Physica A: Statistical Mechanics and its Applications 521 (2019) 828–842.
  * [33] P. M. Bronstein, Agonistic and reproductive interactions in Betta splendens, Journal of Comparative Psychology 98 (4) (1984) 421–431.
  * [34] C. Cantalupo, A. Bisazza, G. Vallortigara, Lateralization of displays during aggressive and courtship behaviour in the Siamese fighting fish (Betta splendens), Physiology & Behavior 60 (1) (1996) 249–252.
  * [35] N. Castro, A. F. H. Ros, K. Becker, R. F. Oliveira, Metabolic costs of aggressive behaviour in the Siamese fighting fish,Betta splendens, Aggressive Behavior 32 (5) (2006) 474–480.
  * [36] T. L. Dzieweczynski, O. L. Hebert, Fluoxetine alters behavioral consistency of aggression and courtship in male Siamese fighting fish, Betta splendens, Physiol Behav 107 (1) (2012) 92–7.
  * [37] E. B. Hågvar, Multiparasitism of the green peach aphid, Myzus persicae: competition in the egg stage between Aphidius matricariae and Ephedrus cerasicola, Entomologia Experimentalis et Applicata 47 (3) (1988) 275–282.
  * [38] D. Uka, T. Hiraoka, K. Iwabuchi, Physiological suppression of the larval parasitoid Glyptapanteles pallipes by the polyembryonic parasitoid Copidosoma floridanum, J Insect Physiol 52 (11-12) (2006) 1137–42.
  * [39] A. M. Lyapunov, The general problem of the stability of motion, International Journal of Control 55 (3) (1992) 531–534.
  * [40] S. Cai, The formula Of Liapounoff function of system of linear differential equations with constant coefficients, Acta Mathematica Sinica, Chinese Series 9 (4) (1959) 455–467.
  * [41] X. Mao, G. Marion, E. Renshaw, Environmental Brownian noise suppresses explosions in population dynamics, Stochastic Processes and their Applications 97 (1) (2002) 95–110.
  * [42] M. Liu, K. Wang, Survival analysis of stochastic single-species population models in polluted environments, Ecological Modelling 220 (9-10) (2009) 1347–1357.
  * [43] F. Wei, S. A. H. Geritz, J. Cai, A stochastic single-species population model with partial pollution tolerance in a polluted environment, Applied Mathematics Letters 63 (2017) 130–136.
  * [44] Y. Cai, Y. Kang, M. Banerjee, W. Wang, A stochastic SIRS epidemic model with infectious force under intervention strategies, Journal of Differential Equations 259 (12) (2015) 7463–7502.
  * [45] M. Liu, K. Wang, Persistence and extinction of a single-species population system in a polluted environment with random perturbations and impulsive toxicant input, Chaos, Solitons & Fractals 45 (12) (2012) 1541–1550.
  * [46] X. Dai, S. Wang, B. Yan, Z. Mao, W. Xiong, Survival analysis of single-species population diffusion models with chemotaxis in polluted environment, AIMS Mathematics 5 (6) (2020) 6749–6765.
  * [47] S. Sharma, G. P. Samanta, Mathematical analysis of a single-species population model in a polluted environment with discrete time delays, Journal of Mathematics 2013 (2013) 1–18.