STABILITY ANALYSIS OF A DELAYED SOCIAL EPIDEMICS MODEL WITH GENERAL CONTACT RATE AND ITS OPTIMAL CONTROL

. In this paper, we formulate an alcohol quitting model in which we consider the impact of distributed time delay between contact and infection process by characterizing dynamic nature of alcoholism behaviours, and we generalize the infection rate to the general case, simultaneously, we consider two diﬀerent control strategies. Next, we discuss the qualities on the model, the existence and boundedness as well as positivity of the equilibrium are involved. Then, under certain proper conditions, we construct appropriate Lyapunov functionals to prove the global stability of alcohol free equilibrium point E 0 and alcoholism equilibrium E ∗ respectively. Furthermore, the optimal control strategies are derived by proposing an objective functional and using classic Pontryagin’s Maximum Principle. Numerical simulations are conducted to support our theoretical results derived in optimal control.


(Communicated by Kok Lay Teo)
Abstract. In this paper, we formulate an alcohol quitting model in which we consider the impact of distributed time delay between contact and infection process by characterizing dynamic nature of alcoholism behaviours, and we generalize the infection rate to the general case, simultaneously, we consider two different control strategies. Next, we discuss the qualities on the model, the existence and boundedness as well as positivity of the equilibrium are involved. Then, under certain proper conditions, we construct appropriate Lyapunov functionals to prove the global stability of alcohol free equilibrium point E 0 and alcoholism equilibrium E * respectively. Furthermore, the optimal control strategies are derived by proposing an objective functional and using classic Pontryagin's Maximum Principle. Numerical simulations are conducted to support our theoretical results derived in optimal control. 1. Introduction. In this paper, we will propose a mathematical model to depict social epidemics. To be precise, we propose delayed and optimal-control consideration on mathematical models to analyze human behaviors related to addictions. As we all know, alcohol which is addictive has done great harm to human health and also seriously affected people's normal life since its advent in the world [3,5,4], nonetheless, alcohol abuse has intensified for many reasons, especially in recent years [40,25]. So more and more health and nutrition experts as well as sociologists committed to the study of this phenomenon from pathological angle or cultural perspective respectively [6,10,44,36,43,35]. In the past two decades, many researchers on infectious diseases and applied mathematics had a keen interest in this phenomenon, they mainly constructed mathematical models and used mathematical methods to carry out macro research on alcohol abuse, for example, theoretical analysis of the stability in related models [41,37,23,38,45].
We formulate an alcoholics quitting model as the follows (1) In this model, the total population N (t) concerned is partitioned into four compartments: the susceptible compartment S(t) which refers to the persons never drink or drink moderately without affecting the physical health; the alcoholisms compartment A(t) which refers to the persons who binge drinking and affect the physical health seriously; the treatment compartment T (t) which refers to the persons who have been received treatments by taking pills or other medical interventions after alcoholism, and the quitting compartment Q(t) which refers to the persons who recover from alcoholism after treatment and keep off alcohol hereafter. Π is the birth number/recruitment of the population; µ i , i = 1, 2, 3, 4 is the natural death rate of the people in the above respective compartments (µ i are not equal due to the impact of alcoholism); u 1 is the fraction of the susceptible individuals who successfully avoid to keep off the alcoholisms; u 2 is the fraction of the alcoholisms who take part in treatments, here, 0 ≤ u i ≤ 1, i = 1, 2, and they will be considered as two control variables in section 5; β is the transmission coefficient of the infection for the susceptible individuals from the alcoholisms individuals; δ is the rate coefficient of the person who have received effective treatment and recover from alcoholisms forever.
As professor Kuang pointed out in his outstanding monograph [24], "time delays occur so often, in almost every situation, that to ignore them is to ignore reality", time delay is very important in characterizing the dynamic behavior of related systems. The refs [32,15,33,2,31,42,46] focused on the research of epidemic models with discrete or distributed delays. Due to delays, the authors discussed the effect on stability of the corresponding models and a variety of dynamic behaviors such as hopf bifurcation phenomenon, etc. We notice that there also exists a time delay when susceptible persons are exposed to infection from alcoholics and then become new alcoholics, it's necessary to take the delay into account when modelling alcohol abuse. In model (1), h represents the supremum of the time delay which means the time between the alcoholics come into contact with susceptible persons and susceptible persons are then infected to become new alcoholics, we assume the susceptible persons will leave the susceptible compartment immediately when the effective contact takes place with the alcoholism persons. The expression h 0 p(τ )f (S(t − τ ), A(t − τ ))A(t − τ )dτ means the alcoholics have different infection intensities on the susceptibles at different time t when t ∈ [0, h], which is embodied by the delay kernel function p(τ ) [33,42]. The delay kernel function p(τ ) works just as a probability density function, so we further assume that p(τ ) is continuous on [0, h] satisfying h 0 p(τ )dτ = 1. This less demanded property of p(τ ) is critical for the subsequent mathematical analysis.
At the same time, contact rate is also very crucial in characterizing the nature of internal dynamics in the modeling of alcohol abuse. We notice that the authors almost all adopted specific incidence in modelling the dynamics of alcoholism, for example, ref. [41] used bilinear incidence rate, while refs [37,23] used standard contact rate to describe the infection between the susceptibles and the alcoholisms. But as we all know, these specific contact rates have their limitations in application, because of the different ways of contact and demographic differences [46]. In our model, we extend the contact rate into the form of an abstract function of which the aforementioned incidences can be involved, in order that our proposed model can be applied in different environments where drinking behavior takes place.
Optimal control theory and methods play increasingly important role in biological and mathematical applications, many experts on biomathematics are committed to this fascinating field. Specific work in this regard, please refer to [8] and the widely cited monograph [28] as well as the references therein. For practical needs, we still premeditate two treating methods as [45] did, namely, prevention of susceptibles from alcoholisms and treatment on alcoholisms as control variables, hence, we will derive a delayed SAT Q-type model including controls. We will propose an objective functional which considers not only alcohol quitting effects but also the cost of alcoholics. Then, we consider a range of issues related to the delay optimal control with the method of Pontryagin's Maximum Principle similar to the way as what the documents [11,12] did.
In short, compared with above-mentioned alcoholism models, our model has the following main characteristics which can be seen some alterations or improvements: • We consider the distributed time delay; • We generalize the infection rate into the form of abstract function; • We consider two different control strategies, i.e, prevention and treatment, for the sake of practical significance.
Furthermore, it's necessary for us to clarify the difference of our model from the ones in references [32,15,33]. Obviously, the first but not the most important aspect is, our model in this paper focus on alcoholism dynamics while the others are not, additionally, our model involves four compartments with two different control strategies while they considered only three ones without any control strategy. In [32,33], the author considered an SIR epidemic model which involved distributed time delay between S(t) and I(t), but the time lag occurred only in the infected persons, and the author thought that the susceptible persons of being infected was a gradual process, which is different from alcohol abuse in our paper. While in [15], the authors considered only discrete delay between S(t) and I(t) in a threedimensional SIR model of infectious diseases with a general contact rate just as we do, and they considered mortality of the susceptible persons who later became infected during infection period. We don't consider the above-mentioned mortality because this death rate is very low in the dynamic behavior of alcoholics and can be almost negligible. In a certain sense, our model is a promotion in [15].
The organization of this paper is as follows. In the next section, firstly, we give some necessary conditions needed later, then we discuss some fundamental properties of our model and get the expression of the basic reproduction number R 0 , which carries much important biological and mathematical significance [7]. In section 3, we prove the global stability of alcohol free equilibrium E 0 by constructing an appropriate Lyapunov functional. In section 4, under additional reasonable conditions, we construct another appropriate Lyapunov functional to prove the global stability of alcoholism equilibrium E * . In section 5, optimal control strategies by the classic method of PMP (Pontryagin's Maximum Principle) are implemented on our model in which we change the general incidence rate into bilinear one and consider only discrete time delay for the sake of feasibility and simplicity on calculations. In section 6, we conduct numerical simulations to verify the theoretical results in section 5. We come to conclusions and some discussions which are necessary for follow-up work in the last section.

2.
Preliminaries. In this section, we will give some conditions and assumptions which are necessary for the analysis, and then study its important properties such as the boundness and positivity of it's solutions.
According to the facts and the need of analysis later, the incidence function f (S, A) is assumed to be continuously differentiable in the interior of IR 2 + and satisfies the following fundamental hypotheses: ∂A ≤ 0, for all S ≥ 0 and A ≥ 0, which include the case of mass action βSA (or so-called bilinear incidence, see [2,31]), saturating incidence βS A 1+cA [42], and also standard (or proportional) incidence β SA S+A [21]. Actually, these hypotheses are also reasonable and consistent with the reality. The first hypothesis is obvious, the second one means the more susceptible persons, the more infectious events will occur. It is necessary to explain the rationality of the third hypothesis, as we all know, the suspected persons are to be vigilant to avoid effective contact with the alcoholics when the alcoholics around them gradually increased.

Fundamental properties of solutions.
It is important to show positivity and boundedness for the system (1) as they represent populations. Firstly, we show the positivity of the solutions. As we can see that the first two equations have nothing to do with the variables T (t) and Q(t), in other words, when we have investigated the properties (here mainly refers to positivity and boundedness and stability) of variables S(t) and A(t), through the last two equations in (1), we can derive the corresponding results for other two variables T (t) and Q(t). For the sake of convenience, we only need to consider the following subsystem: (3) Let C = C([−h, 0], IR 2 ) be the Banach space of continuous functions mapping the intervals [−h, 0] into IR 2 with the topology of uniform convergence. By the basic theory of functional differential equations [14], it is easy to see that there exists a unique solution (S(t), A(t)) of the system (3) with initial value (ϕ 1 , ϕ 2 ) ∈ C. For the well-known reasons that the initial population is non-negative, we assume the initial conditions for (3) satisfies: Theorem 2.1. The solutions of system (3), subject to condition (4), remain nonnegative and bounded for all t ≥ 0.
Proof. Firstly, we prove the solution S(t) is positive for all t ≥ 0. We will show this fact by the method of contradiction. If not, we assume t 1 ≥ 0 is the first time such that S(t 1 ) = 0, then by the first equation of (3) under the hypothesis i), we can get S (t 1 ) = Π > 0, it means there exists a sufficient small number > 0, such that . This contradicts the fact that S(t) > 0, for t ∈ [0, t 1 ). It follows that S(t) > 0, for all t > 0. Next, we begin to prove the positivity of A(t). From the second equation of (3), we get Then we come to prove boundness of the solutions. It follows from the first So we can get Thus, we complete the proof.
To make sure what we concern about is meaningful, one another important issue that we should prove is the existence of the solutions to system (3).
Obviously, system (3) has always an alcohol free equilibrium E 0 ( Π µ1 , 0). According to the definition [15,7], we can define the basic reproduction number R 0 of our model as represents the number of susceptible individuals at the beginning of the "infectious" process and f ( Π µ1 , 0) represents the value of the function f when all individuals are susceptible. Hence, the above-mentioned R 0 is biologically well defined.
It is easy to see that if R 0 ≤ 1, the alcohol free equilibrium E 0 ( Π µ1 , 0) is the unique steady state, which corresponds to the extinction of alcoholism. The following theorem presents the existence and uniqueness of endemic equilibrium if R 0 > 1. (3) always exists for all parameter values.
(2) If R 0 > 1, then the system (3) has a unique alcoholism equilibrium of the form Proof. The conclusion of issue (1) is obvious, so we only need to prove issue (2). The steady state of the system (3) satisfies the following equations From the second equation above, we know A = 0, or f (S, A) = u2+µ2 1−u1 . When A = 0, it corresponds to alcohol free equilibrium point E 0 ( Π µ1 , 0). When A = 0, from the above two equalities, we can get ∂f ∂A > 0 based on the hypothesis ii) and iii) in subsection 2.1. Therefore, there exists a unique alcoholism equilibrium E * (S * , A * ) with S * ∈ (0, Π µ1 ) and A * > 0. The proof is completed.
3. Global stability of the alcohol free equilibrium. In this section, motivated by refs [22,32,15,33], we will construct an appropriate Lyapunov functional to prove global stability of the alcohol free equilibrium.
Theorem 3.1. The alcohol free equilibrium E 0 of the system (3) is globally asymptotically stable when R 0 ≤ 1.
Proof. Consider a Lyapunov functional where S 0 = Π µ1 , we can easy to verify that V 1 (t) is a non-negative functonal. Next, we give detailed calculations on the total derivative of V 1 (t) with respect to system (3) as the following: The last inequality is based on hypothesis iii). Under the hypothesis of ii), we can easily derive the following result: Therefore, when R 0 ≤ 1, we have dV1 dt | (3) ≤ 0. We note that dV1 dt | (3) = 0 iff S = S 0 and A(R 0 − 1) = 0 are both established. Thus, we will make discussions according to the following two cases: 1) when R 0 < 1, then obviously S = S 0 , A = 0; 2) when R 0 = 1, then by first equation of (3) and hypothesis ii), we can easily derive the same conclusion: S = S 0 , A = 0.

4.
Global stability of the alcoholism equilibrium. Now, we should establish a set of condition which is sufficient for the global stability of the alcoholism equilibrium E * (S * , A * ). Thus, we give the following further assumption It is important to note the following remark. , where β, δ, δ 1 , δ 2 , δ 3 ≥ 0 are constants. Recently, this more generalized incidence function is used in [1,29,19].
Firstly, we give a lemma which is a useful tool in the following derivations.
Proof. The proof is obvious, so we omit it here. Proof. In order to facilitate the calculation we introduce the following curt notations: Based on the Lyapunov functionals constructed in [15,17,18], we consider the following dX has the global minimum at x = x * and ψ(x * ) = 0. Then, ψ(x) ≥ 0 for any x > 0. From Lemma 4.2, we deduce that V 2 (t) ≥ 0.

5.
Delay optimal control in a common case. In this section, we begin to investigate the optimal control strategies to limit alcohol abuse and purify the social environment based on system (1), in which we will assign specific incidence of infection as f (S(t), A(t))A(t) = βS(t)A(t), i.e., the so-called bilinear incidence rate, since it is almost the simplest case to be calculated and satisfies the conditions given in section 2.1. Furthermore, we will let p(τ ) = δ(t − τ ), here δ denotes Dirac function, so the distributed delay is reduced to discrete one. The reason for doing these simplified process is in order to ensure the feasibility and convenience for follow-up calculations. To do that, we set prevention rate u 1 and treatment rate u 2 be variables u 1 (t) and u 2 (t) whose values range are both in [0, 1]. u 1 (t) is the first control variable usually implemented by psychological impact and behavioral constraints, specifically, we can take measures such as propaganda and education, etc. u 2 (t) is the second control variable which corresponds to takeing appropriate treatment measures, such as taking pills or seeking other medical help. Our goals are to minimize the number of alcoholisms and to maximum the number of the recovered and the susceptibles. Based on consideration of the reality of the situation on alcohol abuse, our desired target is slightly different from the general disease control problems in which reducing the number of the the susceptibles is also involved [8,28]. However, just as a coin has two sides, while we control alcoholism, it will also generate much associated costs. Therefore, it's advisable to balance between the costs and the alcohol effects. In view of this, our optimal control problem is to maximize the objective functional which is formulated as which subjects to system with initial conditions Here, u i (t) ∈ U = {(u 1 , u 2 )|u i (t) is measurable and 0 ≤ u i (t) ≤ 1, ∀t ∈ [0, t f ]}, t f is the end time to be controlled, U is an admissible control set, c i , i = 1, 2, are weight factors (positive constants) that adjust the intensity of two different control measures. The meaning of ϕ i (θ), i =, 1, 2, 3, 4, is the same as before.
Next, we will investigate the existence of the optimal control of the abovementioned problem. From aforementioned statement in section 2, we can easily know the existence of solutions to system (7) with the given initial values (8), so we just need to focus on solving control-related problems, for example, the existence of optimal control and the specific characteristic of the optimal control, etc. 5.1. The existence of optimal control. Theorem 5.1. There exists an optimal control pair u * = (u * 1 , u * 2 ) ∈ U such that J(u * 1 , u * 2 ) = max{J(u 1 , u 2 ) : (u 1 , u 2 ) ∈ U } subjects to the control system (7) with initial conditions (8).
Proof. We denote X(t) = (S(t), A(t), T (t), Q(t)) T and L(t; X(t), u 1 (t), u 2 (t)) = S(t) To prove the existence of an optimal control, according to the classic literatures [9,30], we have to show the following: (1) The set of controls and corresponding state variables is nonempty.
(2) The control set U is convex and closed.
(3) The right hand side of the state system is bounded by linear function in the state and control variables.
(4) The integrand of the objective functional is concave on U .
(5) There exist constants d 1 , d 2 > 0 and α > 1 such that the integrand L(t; X(t), u 1 ; u 2 ) of the objective functional satisfies In order to verify these conditions, we use the result given in Section 2.2 to establish the existence of solutions of system (7), which gives condition 1. The control set is convex and closed by definition, which gives condition 2. The state system is bilinear in u 1 , u 2 and the solutions are bounded, hence, the right hand side of system (7) satisfies condition 3. Note that the integrand of our objective functional is concave, so the condition 4 is satisfied. It remains to verify the last condition. By Theorem 2.1, we know the solutions to system (7) are bounded, so there exists a positive number M suth that |X(t)| ≤ M , for any t ≥ 0. Therefore, we have Then the last condition is satisfied by choosing d 1 = min{ c1 2 , c2 2 }, d 2 = 2M and α = 2. The proof is completed.

5.2.
The characterization of the optimal control. With the existence of the optimal control pairs established, we now continue to discuss the theorem that relates to the characterization of the optimal control. For the sake of convenience, we note X(t) as before and X τ (t) = X(t − τ ), U (t) = (u 1 (t), u 2 (t)) T , λ(t) = (λ 1 (t), λ 2 (t), λ 3 (t), λ 4 (t)) T . To do this, we begin by defining a Hamiltonian H with penalty terms for the control constraints as follows: where w ij (t) ≥ 0 are the penalty multipliers satisfying w 11 (t)u 1 (t) = w 12 (t)(1 − u 1 (t)) = 0 at optimal control u * 1 , and w 21 (t)u 2 (t) = w 22 (t)(1 − u 2 (t)) = 0 at optimal control u * 2 . Our main method to study the optimal control problem with delay in the state (6)-(8) is given by the Pontryagin Maximum Principle [39]. A necessary condition for optimal control problems can be found in [13] and it's applications [11,12]. There exists a continuous function λ(t) on [0, t f ] satisfying the following three equations, that is, the state equation the optimality condition and the adjoint equation where H λ , H U , H X and H Xτ denotes the derivative with respect to λ(t), U (t), X(t), and X τ (t) respectively. Further, it should be noted that χ I (t) is usually called indicator function or characteristic function which is defined as χ I (t) = 1, when variable t is in interval I, 0, when variable t isn t in interval I.
Now we apply the necessary conditions to the Hamiltonian H, we can derive the following theorem.
Theorem 5.2. Given an optimal control pairs (u * 1 , u * 2 ) and solutions S(t), A(t), T (t), Q(t) of the corresponding state system (7), there exist adjoint variables λ i (t), i = 1, 2, 3, 4, satisfying Furthermore, (u * 1 , u * 2 ) are represented by Proof. According to the formula for the adjoint equation, we first differentiate the Hamiltonian operator H, with respect to states. Then the adjoint system can be written as with a series of calculations, we can derive the costate system (11). Conditions (12) of adjoint equations is obviously given by λ i (t f ) = 0, i = 1, 2, 3, 4.
To obtain the necessary conditions of optimality (13), by formula (10), we differentiate the Hamiltonian operator H, with respect to U = (u 1 (t), u 2 (t)) and set them equal to zero, then Solving for the optimal control, we obtain To determine an explicit expression for the optimal control without w 11 and w 12 , a standard optimality technique is utilized [9]. We consider the following three cases. (i) On the set {t|0 < u * 1 (t) < 1}, we have w 11 (t) = w 12 (t) = 0. Hence the optimal control is (ii) On the set {t|u * 1 (t) = 1}, we have w 11 (t) = 0. Hence This implies that (iii) On the set {t|u * 1 (t) = 0}, we have w 12 (t) = 0. Hence This implies that Combining these results, the optimal control u * 1 (t) is characterized as Using the similar arguments, we can also obtain the other optimal control function u * 2 (t) = min 1, max 0, The proof is completed. We point out that the optimality system consists of the state system (7) with the initial conditions S(0), A(0), T (0), Q(0), the adjoint (or costate) system (11) with the terminal conditions (12), and the optimality condition (13). Any optimal control pairs must satisfy this optimality system. Finally, for the sake of convenience in numerical simulation which will be implemented in the next section, we get the optimality system as the follows and u * 2 (t) = min 1, max 0, 6. Numerical simulation in optimality system. In this section, we use the efficient numerical method developed by Hattaf and Yousfi in [16] in order to solve the optimality system (15) because this system is represented by delay differential equations (DDEs) and is a two-point boundary-value problem with the initial conditions are given for state variables at time t = 0, while terminal conditions are given for adjoint variables at time t = t f . Hence, the state system is solved forward in time and the adjoint system is solved backward in time. For more details, see the algorithm of the method in [page 6, [16]]. Furthermore, this iterative method is recently used in [26] to simulate a delayed SIRS epidemic model with vaccination and treatment. For this simulation, we take initial values S(0) = 400, A(0) = 50, T (0) = 20, Q(0) = 10 and we use the parameters values given in the following table: The period of the control course considered is 50 weeks(almost one year) and the two weight factors are c 1 = 10 3 and c 2 = 10 2 respectively. Remark 2. For the death rate of susceptible individuals µ 1 we cite [37]; for β, δ and τ , we take average values referring to [20], in which a range of values are pointed out. Other parameter values and initial values are estimated at a reasonable sense.
In the following, firstly, we will give the graphes of evolution with time in different compartment populations with or without controls to reveal the effect of the control strategies considered in our paper. Secondly, we will show the graphes of the optimal control strategies which have been derived in section 5 to give specific guidance in our actual operating.   In figure 1, with respect to S(t) compartment, the population without controls will decrease monotonically for the infection from alcoholics; while with the controls, almost in the beginning five weeks, the population will still be reduced because of the delay between the effective contact and infection, but then the susceptible population will be increased till the end of the experiment since the control strategies play a role. With respect to A(t) compartment, the population without controls will increase monotonically, conversely, the population with controls will decrease continuously. With respect to T (t) compartment, the population without controls will decrease continuously to almost zero, while with controls, the population will increase to a peak value in the beginning eight weeks and then reduce to a fixed point. The law of population variation in Q(t) compartment is similar to the case of T (t) compartment, so we will not repeat here. In short, the two control actions are both contribute to effectively reducing alcoholism, but control effect will be affected by the time delay between the effective contact and infection, which confirms the objective role of time delay in a certain sense.
In figure 2, we can derive the following conclusions: (1). from the perspective of control effectiveness and economic factors, we should intensify effort of the second control u 2 (t)(i.e., treatment control) and weaken the first control u 1 (t)(i.e., prevention control), to which one possible explanation is that treatment is more operational compared with prevention. (2). the optimal control strategy of u 1 (t) is that we should gradually enhance the effort to 0.1 in the first seven weeks and then weaken the effort to 0.034 thereafter; while the optimal control strategy of u 2 (t) is that we should gradually enhance the effort to 0.56 in the first seven weeks and then weaken the effort to 0.38 till end.
7. Conclusions and discussions. In this paper, we formulate a four-dimensional SAT Q-type kinetic model of alcoholism behaviour which includes three important characteristics, i.e., the general contact rate between the suspected persons S(t) and the alcoholism persons A(t), the time distributed delay between effective contact and infection process of them, and also two different control strategies to reduce or eliminate the alcohol abuse which is of great importance nowadays. Due to the characteristics of our model, we reduce it to the two-dimensional case to simplify the analysis and computation involved. Under appropriate conditions, we obtain the existence of two equilibrium, and prove their global asymptotic stability respectively by virtue of Lyapunov method. Comparing with the results from [45], the conclusions derived in this paper tell us that the time delay above-mentioned does not affect the stability of the system under some suitable conditions, in other words, the system is so-called absolutely stable. In the end, we use classic Pontryagin Maximum Principle to obtain the optimal control strategies in the context that we only consider bilinear incidence rate as well as discrete delay to facilitate the subsequent calculations. From the results of optimal control, we can see that the time delay considered contributes to relieving alcohol abuse. We also carry out numerical simulation to verify our optimization results.
Much work to continue research and analysis on the dynamics of alcoholism is still left, for example, we can consider the effect of awareness programs by adding a population compartment of aware population into our model just as the published paper [34] did, which is very effective and objective in studying and controlling alcohol abuse for the purposes of economic and social development in today's world, so it's an important and prospective work to be continued. On the other hand, the same to other infectious diseases dynamic behavior, alcoholism behaviour is also filled with a lot of uncertainties, hence, in our model, we can consider white noise or colored noise which are caused by the social environment, cultural customs or human-induced, etc, to do that, we can derive the corresponding stochastic dynamics model, then, we can further use various stochastic approaches to in-depth study. Our follow-up work will commence in accordance with the above topics.