Dynamics analysis of Mackey-Glass model with two variable delays

1.   Introduction

Synchronization ^[1,2,3,4,5] and propagation ^[6,7,8,9,10], as the typical interesting and important dynamical behaviors in complex networks ^{[11,12,13,14,15,16,17,18,19,20]}, arouse a hot study in various fields in recent decades. When a disease is discovered and diffused from the affected area, the information about and appraisal of an epidemic is also transmitted along the connective networks. In other words, if there exists relationship between two individuals, one person can get the message from another, so that the message about disease can diffuse along with the relationship networks. While, the transmission of epidemic is through the contact network among individuals. Furthermore, many outstanding works are taken which have looked at epidemic spread as an important aspect of transmission behavior in complex networks from the mathematical biology field.

A thought is taken it for granted that synchronization and propagation almost disappear simultaneously in a complex network, but many examples and phenomena show that situations are opposite in fact. Especially when a kind of epidemic disease transmission, the masses' behavior always reach to consensus and tackle diseases, which can be claimed to achieve the synchronous behaviors for nodes in the complex network. As a matter of fact, we always see or do the same things such as washing hands with soap frequently, doing exercise to enhance our health and avoiding to go to the crowded place and so on, to refrain from infected and to combat the disease when we get the information about epidemic diseases. This kind of consistent behaviors in a complex network is so called behavior synchronization aroused by the epidemic spreading. In recent years, there are many excellent works concerned about such kind of synchronization induced by contagious diseases ^{[21,22,23,24,25]}, while in which synchronization and spreading of epidemic disease are discussed in combination from the aspect of numeric. Therefore, it is so interesting about the interaction among synchronization and epidemic spread in complex networks that we will also go on discussing the behavior synchronization induced by propagation of disease in this work.

From the view of mathematical biology, several interesting models have been analyzed qualitatively on dynamical synchronization aroused by epidemic spread. The adaptive mechanism between dynamical synchronization and epidemic behavior as well as the interplay between collective behavior and spreading dynamics were discussed respectively, where the effect about individuals' awareness was introduced in epidemic models in the later one ^[26,27]. And then mutually restrictive conditions among dynamical synchronization and epidemic transmission were derived from sternly theoretical methods. Follow on, with the assistance of the quenched mean-field theory, behavioral synchronization induced by epidemic spread was addressed under the constructed mathematical models without delay, with a coupling delay and with double delay ^[28]. And Sun et al. ^[29] taken an exploration and simulation of epidemic spread and its control in multiplex networks, and found that the epidemic control strategy from the perspective of behavioral control is extremely relevant for epidemic control.

While those works are all about complete synchronization patterns induced by epidemic spread, the cluster collective behavior in complex networks are more general. And there are many successful efforts on the cluster synchronization in complex networks, such as ^[30,31] and references therein. However, few works paid attention to grouped collective behaviors induced by disease propagation in a complex networks with community structure. The collective behaviors induced by epidemic spread are firstly considered under community structure of complex networks, and obtained corresponding results ^[32]. So universal and significant as the cluster collective behaviors are that we will consider the grouped collective behaviors aroused by epidemic transmission in our work. On the other hand, because of distinct survival conditions and food resources distributed in space, there are agammaegation of individuals in different areas and then forming patches ^[33,34]. And it is common phenomenon that individuals can take dispersal or migration in the patch environment, thus in our discussion there will also consider the effect of dispersal among patches in the epidemic model, and further consider the grouped collective behaviors induced by patchy epidemic spread. Different from other cluster characteristic collective behaviors, the cluster numbers are corresponding to the number of patches in our consideration, and the information about diseases in each patch will be utilized to design the behavioral state feedback pinning controllers in the behavior dynamical systems which familiar with the control strategy ^[35].

Inspired by the above discussions, the cluster collective behaviors in behavioral synchronization systems aroused by epidemic spread under patchy population will be studied in our work. In the problem we considered, the epidemic models with input rates and remove rates under patchy dispersal environment different from the above mentioned, are constructed in view of the quantity of population, while the total population is guaranteed as a constant. In which we assume that the dispersal only appear among the susceptible in patches. Then we consider the relation network of individuals under this constant population, and assume that the relationship is inheritable, even though there are input and remove rates in epidemic models, which means the fixed relation network. In other words, the new input one will fill the lost relation position in the relation network. It should be noticed that the spatial structures of individuals are not concerned, just concerning the relation. And then, dividing the relation network into several subetaoup same as the number of patches, named as behavioral pattern relation group induced by patches. Called the individual belongs to the group $k$ if it lies in the $k$th behavioral pattern relation group. On the basis of above, the main purpose of our work is to achieve behavioral consistence in each group under the transmission of disease by using so called behavior feedback controlling strategy related with the information about diseases and the whole behaviors for all individuals.

The main contribution of this work are listed as follows. First, we construct the coupling system about patch epidemic disease-behavior synchronization models based on population, in which the effect of dispersal only appears in the susceptible while the infected are limited to migration. And a function of contact willing related with information about disease in each patch and whole behavior state is added before the incidence as a kind of cofactor. Then, we design the disease relative feedback pinning controllers by using the diseases' information for each patch and the states of behavior in each group. In other words, a part of nodes in each group are needed to control. Finally, we obtain conditions about global asymptotical stability for disease-free equilibrium by using Lyapunov stability analysis method, and prove the existence of positive equilibrium, and derive the corresponding criteria. On the other hand, we verify the global asymptotical stability of the positive equilibrium for epidemic systems by Lyapunov function, and meanwhile, obtain the criteria that achieve the cluster behavioral synchronization, related with the cluster relation network structures and patch characteristic.

The rest of this paper is organized as follows. The mathematical models and necessary preliminaries are shown in Section 2. In Section 3, the global stability of disease-free equilibrium and the existence of positive equilibrium are analyzed in the patch epidemic-behavior synchronization systems, and the global stability of positive equilibrium as well as the cluster behavioral synchronization under the designed feedback pinning control strategy are discussed. The numerical simulation is performed in Section 4 to demonstrate effectiveness of the theoretical results. Finally, the conclusion is presented in Section 5.

Notations: Trough out this paper, $\mathbb{R}^{+}$ is the set of nonnegative real numbers, $\mathbb{R}^{n}$ represents $n$ dimensional Euclidean space, and $\mathbb{R}^{n\times m}$ denotes the set of $n\times m$ dimensional real matrix. A matrix $P$ is claimed as positive (negative) if $P > 0\; (<0).$ And the unit matrix is denoted as $E_{m\times m}$ with $m\times m$ dimension.

2.   Models and preliminaries

2.1.   Mathematical models

As shown in Figure 1, the mutual influence between patch epidemic and individual behaviors includes two processes, such as, the one is the transmission of information about epidemic diseases, and the behavioral contact among individuals. In which the message about patch epidemic disease will be transferred along with the relation network of the individuals, meanwhile the information about patch disease will be used to design the state feedback controllers for individuals' behavior pattern. The relation network mentioned here is a kind of network about information exchange among individuals, and there exist also reinforcement of transmission, such as more frequent exchange. On the other hand, the spreading for epidemic diseases occurs by the full probability contact network, where considering the bilinear incidence, while the behavior pattern will be as one part for a kind of cofactor function to affect the actual infected population. In this framework, the mutual interaction coupling systems can be formed. It should be noticed that we model the patchy epidemic spreading based on the number about total population of the susceptible and the infected in each patch, and it is taken into consideration that the uniform contact or identical probability contact among the susceptible and the infected. Thus, the concrete topological construction about the contact network leading to disease transmission is not reflected in the population patch epidemic model.

In this paper, we consider the fixed topology structure of the relation network $A = (a_{ij}).$ Because the input and remove are considered in the epidemic model (shown below in the system (2.1)), we assume that the new individuals fill up the leaves one's relation site and remain the origin relation unchange. Even though dispersal is also considered in the epidemic model, we assume that the relation network among individuals is unchanged with the time going, and the subetaoups of relationship are divided from the patch disease information not the individuals' numbers in each patch, which means that the numbers for the subetaoup of relationship are same as the patch numbers. In the same subetaoup of relationship, there will be a enhancer coefficient that reinforces the relationship in the same subetaoup. We name these groups consisting of the individuals in such relationship subetaoup as behavior pattern relation group induced by the corresponding information of patch diseases. Which can be understandable that some individuals in different places that are feedback by one patch so called that they belong the same behavior pattern relation group induced by the information of patch diseases. In each behavior pattern group induced by patch diseases, we will take the patch disease information to design the feedback control tactic to achieve a desired behavior consistence. Maybe these feedbacks are added to a part of individuals in the behavior pattern group.

We now construct the following $m$ patch epidemic models with dispersal

Where $S_{k}(t)$ and $I_{k}(t)$ represent the susceptible and infective population in $k$th patch at $t,$ respectively. And the total population $N(t) = \sum_{k = 1}^{m}(S_{k}(t)+I_{k}(t)).$ In the epidemic models, we assume that the dispersals of the susceptible population are balance, which means that the dispersal matrix $B = (b_{kl})$ satisfies $\sum_{k = 1}^{m}\sum_{l = 1}^{m}b_{kl}S_{l}(t) = \sum_{k = 1}^{m}\sum_{l = 1}^{m}b_{lk}S_{k}(t),$ where $b_{kl}$ is the rate of dispersal with $b_{kl}\geq0,$ and $B = (b_{kl})$ is assumed to be irreducible. $\mu S_{k}^{0}$ is the constant input in the $k$th patch, $\mu$ is the natural death rate, and the diseased death is not considered in this model. And suppose that $\sum_{k = 1}^{m}S_{k}^{0} = N(t).$ In addition, $\phi_{k}(x(t), I_{k}(t))$ represents the contact willing function that is relative with the behavior states and information about patch diseases, which reflects the case that the more similar between individuals' behavior pattern, the more willing to contact with each other, or the more severe the disease, the less likely they are to be exposed. And the actual expression will be defined later.

Discuss the models (2.1), and with the assumption $\sum_{k = 1}^{m}S_{k}^{0} = N(t)$ assistance, we find that

Therefore, the population $N(t)$ is a constant, denoting as $N,$ i.e., $N(t)\equiv N.$

Based on the above finding, we introduce the individuals'(or nodes') behaviors, and then construct the coupling systems about cluster individual behavior network systems and patch epidemic models with dispersal under constant population, as following

In which, corresponding to the number of patches there also be divided into $m$ behavior pattern groups induced by the information of patch diseases under the relation network. Let $\mathbb{K}$ represents the $k$th behavior pattern group induced by $k$th patch epidemic, and $\mathbb{K} = \{r_{k-1}+1, \cdots, r_{k}\}$ is the nodes set of individuals in $k$th behavior pattern group. In the system (2.3), $x_{i}(t)\in\mathbb{R}^{n}$ represents the $i$th individual's behavior state vector, $f_{k}(x_{i}(t)):\mathbb{R}^{+}\times \mathbb{R}^{n}\rightarrow \mathbb{R}^{n}$ denotes the local behavior dynamics in the $k$th group. If the node $i$ and $j$ exist relationship then $a_{ij} > 0,$ otherwise $a_{ij} = 0.$ $c_{ij}$ denotes the intra-group enhancement factor, which means that $c_{ij} = c_{k} > 0$ if $\bar{i} = \bar{j},$ where $\bar{i}$ is the group that $i$ belongs, otherwise $c_{ij} = 1$ when $\bar{i}\neq \bar{j}.$ $u_{i}(t)$ is the state feedback controllers related with information about patch diseases which will be designed later.

And from the models (2.3), one can know that local behavior of nodes in the whole network are non-identical and each subetaoup has distinct behavioral target pattern. What's more, the division for nodes relies on their relationship network, do not on their actual spatial position, which means that if individuals in distinct places are impacted from the same patch's disease information, we will classify them as the same relationship subetaoup. Therefore, in models (2.3), the first one equation represents the node $i$'s dynamical behavior process in $k$'s relationship subetaoup $\mathbb{K}$ that aroused by the $k$'s patch disease. Then when the patch disease spreading, the behavioral control $u_{i}(t)$ for each node $i$ in corresponding subetaoup will be activated to make node's behavior pattern trend the desired state.

So far, we define the expression of contact willing function $\phi_{k}(x(t), I_{k}(t))$ as

where the $\alpha_{k}\in (0, 1)$ is the consciousness strength in the $k$th patch. Let $x(t) = (x_{1}^{T}(t), \cdots, x_{N}^{T}(t))^{T}$ and $y(t) = (y_{1}^{T}(t), \cdots, y_{m}^{T}(t))^{T},$ and denote the behavioral error for each node as $e_{i}(t) = x_{i}(t)-y_{k}(t),$ then there are $e(t) = (e_{1}^{T}(t), \cdots, e_{N}^{T}(t))^{T},$ in which $y_{k}(t)\in\mathbb{R}^{n}$ is the target in each subetaoup. For writing convenient, rewriting the $\phi_{k}(x(t), I_{k}(t))$ as $\phi_{k}(t),$ i.e., there are

The system (2.3) is a kind of compound models mixes with the behavioral equation for individuals' states and patch epidemic spread equations. $x_{i}(t), i\in\mathbb{K}$ represents the individual's state vector at $t$ for the node belongs to the $k$th behavior pattern relation group, which reflects the individual's actions such as washing hands, wearing a mask and having a rest, etc.. $f_{k}(\cdot)$ is the local dynamics for the group $k$ that shows the summation of distinct behaviors for the individuals in the subetaoup, in other words, the nodes are not identical. And $f_{k}(\cdot)$ for each $k$ is active, dynamic, as well as complex. $S_{k}(\cdot)$ and $I_{k}(\cdot)$ are the population of the susceptible and the infective, and bilinear incidence is adopted here. The contact willing function $\phi_{k}$ for each patch, determined by patch information and total behavior state error, is introduced as a factor to affect the incidence. $\phi_{k}$ is affected by the behavior equations, in which the more less the behavior error is, the more willing to contact with others. The controllers in the behavior equation of system (2.3), $u_{i},$ are added to the single node and will use the patch epidemic information and single node's behavior errors, which means that only when epidemic transmission the controllers can work to synchronize behaviors. What's more, the controllers are even pinning that only a part of nodes need to be controlled. It is meaningful to design the controllers $u_{i}$ in this way. Notice a fact that two processes occur, the one is the transmission of information about patch disease along with the relation network $A,$ which will influence individuals' behavior pattern by external controllers. Another is the epidemic spread through the contact network adopting the bilinear incidence affected by behavior related cofactor, which will impact the epidemic spread from contact willing function. Therefore, in models we realize the coupling between behavior equations and patch epidemic spread equations by the external controllers $u_{i},$ as well as the $\phi_{k}.$

2.2.   Preliminaries

Definition 1. It is said a network contains $N$ nodes achieve cluster synchronization if the network consists of $N$ is divided into $m$ subetaoups, such as $i\in\mathbb{K} = \{r_{k-1}+1, \cdots, r_{k}\}, k = 1, \cdots, m,$ for any $i$ and $j$ belong to the same cluster, there are $\lim_{t\rightarrow\infty}||x_{i}(t)-x_{j}(t)|| = 0.$

Remark 1. Usually, we do not compute the errors between any two nodes $x_{i}(t)$ and $x_{j}(t)$ directly. In order to achieve cluster synchronization, it is only needed to obtain $\lim_{t\rightarrow \infty}||x_{i}(t)-y_{k}(t)|| = 0,$ for $i\in\mathbb{K}$ under the given synchronization target $y_{k}(t)$ for each cluster. And in this work, we want to use the patchy disease information $I_{k}(t)$ to design nodes' controllers in each subetaoup to achieve cluster behavior synchronization. Which means only when disease spreading controllers can be activated to work.

Definition 2. A matrix $L$ is claimed to belong the $A1$, denoted as $L\in A1,$ if $A$ satisfies

1. $l_{ij}\geq0, j\neq i, l_{ii} = -\sum_{j = 1, j\neq i}^{N}l_{ij},$

2. $L$ is irreducible.

If the $L\in A1$ and $L$ is symmetric, it is said $L$ belongs to $A2,$ denoted as $L\in A2.$ Furthermore, for any matrix $L\in \mathbb{R}^{n\times m},$ if satisfies the row-sum zero, it can be called that $L_{n\times m}$ belongs to $A3,$ denoted as $L_{n\times m}\in A3.$

Definition 3. A continuous function ^[36] $f(x(t)):[0, \infty]\times \mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is said to be $QUAD$ class of function, denoted as $f\in QUAD(P, \Delta, \eta),$ if there exist positive definite diagonal matrix $P = diag\{p_{1}, \cdots, p_{n}\}$ and the diagonal matrix $\Delta = \{\delta_{1}, \cdots, \delta_{n}\},$ and a positive constant $\eta$ such that for $f$ there holds

for any $x(t), y(t)\in\mathbb{R}^{n}.$

Which also implies that there exists positive constant $\xi$ and a positive diagnose matrix $P = diag\{p_{1}, p_{2}, \cdots, p_{n}\}$ such that

Remark 2. It can be easily verified that almost chaotic systems satisfied the $QUAD$ class of function, for instant, the famous Lorenz system, the Chua's circuit, Rössler system, Chen system, and so on.

Lemma 1. For any vectors $x, y\in\mathbb{R}^{n},$ there exists a positive constant $\theta,$ such that

The relation between the above inequalities is actually derived from the Cauchy inequality for the vector. Under the positive constant $\theta,$ using the Cauchy inequality can easy to obtain above result ^[37].

3.   Main results

In this section, we will consider the boundedness and nonnegativity of the $m$ patches epidemic models in the system (2.1). Obviously, the boundedness of the system (2.1) has been verified from the discussion in Eq (2.2).

For the nonnegativity of system (2.1), which means that for any given initial values $S_{k}(0)\geq0, I_{k}(0)\geq0,$ there will hold $S_{k}(t)\geq0, I_{k}(t)\geq0,$ for all $k = 1, 2, \cdots, m,$ and $t\geq0.$

Assume that for any initial values $S_{k}(0) > 0,$ there exists $k$ and $t_{1},$ such that

In other aspect,

Therefore, when $t\in[0, t_{1}),$ we get the $S_{k}(t)$ satisfies that

then, when $t\rightarrow t_{1},$ we obtain

But from the assumption $S_{k}(t_{1}) = 0,$ there is contradiction. Therefore, for any initial values $S_{k}(0) > 0,$ there are $S_{k}(t) > 0$ for $k = 1, 2, \cdots, m, t\geq0.$ Combined with the continuity of solution to initial value, when $S_{k}(0)\geq0,$ there always hold $S_{k}(t)\geq0$ for any $t\geq0$ and $k = 1, 2, \cdots, m.$

For the $I_{k}(t),$ we get $I_{k}(t) = I_{k}(0)e^{\int_{0}^{t}(\phi_{k}(s)b_{k}S_{k}(s)-\mu)ds}.$ When the initial values $I_{k}(0)\geq0,$ there exists that $I_{k}(t)\geq0$ for all $t\geq0, k = 1, 2, \cdots, m.$

Remark 3. It should be noticed that it is just a kind of qualitative analysis way to discuss the nonnegativity of solutions for the epidemic system, and the concrete solutions expression for the system can not be derived in this process. It is merely a method to analyze character of sign for solutions that we want to know.

3.1.   The disease-free equilibrium of patch epidemic systems

We now analyze the disease-free equilibrium of the system (2.1). We first introduce some notations. Let $S(t) = (S_{1}(t), S_{2}(t), \cdots, S_{m}(t))^{T}\in\mathbb{R}^{m},$ $b = (\mu S_{1}^{0}, \mu S_{2}^{0}, \cdots, \mu S_{m}^{0})\in\mathbb{R}^{m},$ and

Under the situation of the disease-free equilibrium, we discuss the equation

Therefore, we obtain the compact form,

Then the problem about the solution for Eq (3.1) can be transformed into the solution for matrix Eq (3.2). Which means that we want to obtain the matrix equation's solution such that $S(t) > 0$ for Eq (3.2). In other words, the current aim is to get the positive solution for the system of inhomogeneous linear Eq (3.2).

For the Eq (3.2), assume

Then there exists positive solution for Eq (3.2).

The following is simple proof about the existence of positive solution for Eq (3.2).

Proof. At first, let $B_{D} = diag\{\mu+\sum_{l\neq 1}b_{l1}, \mu+\sum_{l\neq2}b_{l2}, \cdots, \mu+\sum_{l\neq m}b_{lm}\},$ since $b_{kl}\geq0,$ then $\mu+\sum_{l\neq k}b_{lk} > 0, k = 1, 2, \cdots, m,$ therefore, $B_{D} > 0$ and $B_{D}$ is invertible as well as $B_{D}^{-1}\geq0,$ where $B_{D}^{-1} = diag\big\{\frac{1}{\mu+\sum_{l\neq1}b_{l1}}, \frac{1}{\mu+\sum_{l\neq2}b_{l2}}, \cdots, \frac{1}{\mu+\sum_{l\neq m}b_{lm}}\big\}.$

Let $D = E-BB_{D}^{-1},$ then $B = (E-D)B_{D}.$ From the $\textbf{(H1)},$ we obtain that $C\equiv(E-D)b > 0.$ And according to $b > 0, D\geq0,$ we get $\rho(D) < 1,$ where $\rho(D)$ represents the spectral radius of $D$ ^[38]. By the Perron-Frobenius theorem there exists a vector $v\geq0,$ such that $D^{T}v = \rho(D)v.$ Thus, there is $\bf{0} = v^{T}\rho(D)-v^{T}D = v^{T}(E\rho(D)-D) < v^{T}(E-D).$ Combining with $b > 0,$ we have $v^{T}(E-D)b = v^{T}C > 0.$ On the other hand, $v^{T}(E-D)b = (1-\rho(D))v^{T}b,$ which implies $1-\rho(D) > 0.$

To sum up, we can claim that $E-D$ is invertible, and then with the representation of $B,$ $B$ is also invertible. Therefore, for the Eq (3.2), there are

Which means there exist positive solution for the Eq (3.2) when $\textbf{(H1)}$ hold.

Then, denote by $E_{0} = (\bar{S}_{1}^{0}, 0, \bar{S}_{2}^{0}, 0, \cdots, \bar{S}_{m}^{0}, 0)$ the disease-free equilibrium, we have

Therefore, we obtain the following results for the disease-free equilibrium $E_{0} = (\bar{S}_{1}^{0}, 0, \bar{S}_{2}^{0}, 0, \cdots, \bar{S}_{m}^{0}, 0).$

Theorem 1. Assumption that $B = (b_{kl})_{m\times m}$ is irreducible and (H1) hold.

1. When $\bar{R}_{k}^{0} = \frac{\bar{S}_{k}^{0}\alpha_{k}b_{k}}{\mu}\leq1, k = 1, 2, \cdots, m,$ then the disease-free equilibrium $E_{0}$ is globally asymptotically stable.

2. When $\underline{R}_{k}^{0} = \frac{\bar{S}_{k}^{0}\alpha_{k}(1-\frac{1}{N})b_{k}}{\mu} > 1, k = 1, 2, \cdots, m,$ then the disease-free equilibrium $E_{0}$ is unstable, and if there are $\widetilde{R}_{k}^{0} = \frac{S_{k}^{0}\alpha_{k}b_{k}}{\mu+\sum_{l\neq k}b_{lk}} > 1, k = 1, 2, \cdots, m,$ then there exist positive equilibrium $E^{*} = (S_{1}^{*}, I_{2}^{*}, S_{2}^{*}, I_{2}^{*}, \cdots, S_{m}^{*}, I_{m}^{*}).$

Proof. Firstly, denote the Lyapunov candidate function

where $v_{k}$ is the $k$th cofactor of the $k$th diagonal entry for the Laplacian matrix of $\widetilde{B} = (\widetilde{b}_{kl})_{m\times m},$ and $\widetilde{b}_{kl} = b_{kl}\bar{S}_{l}^{0}.$ Due to the $B = (b_{kl})_{m\times m}$ is irreducible, so as the $\widetilde{B},$ then $v_{k} > 0, k = 1, 2, \cdots, m.$

Then, take the derivative of $L_{1}(t)$ along the system (2.1), and according to the equation that disease-free equilibrium $E_{0}$ satisfies as shown in Eq (3.3), we obtain

From the fact that

and the equality holds if and only if $S_{k}(t) = \bar{S}_{k}^{0},$ then there is

Let $\widetilde{b}_{kl} = b_{kl}\bar{S}_{l}^{0},$ then $\widetilde{B} = (\widetilde{b}_{kl})_{m\times m}.$ And then denote that

Therefore, from Eq (3.5) we obtain that

On the other hand, there is $1-x+\ln x\leq0$ and the equal hold if and only if $x = 1.$ Which means that $1-x\leq-\ln x$ holds. Then

Then

where $F_{k}(S_{k}(t), \bar{S}_{k}^{0}) = -\frac{S_{k}(t)}{\bar{S}_{k}^{0}}+\ln \frac{S_{k}(t)}{\bar{S}_{k}^{0}}.$

Then we have

Therefore, there exists the following result according to ^[39]

To sum up, we have $\dot{L}_{1}(t)\leq-\sum_{k = 1}^{m}v_{k}\Big(\mu-\bar{S}_{k}^{0}\phi_{k}(t)b_{k}\Big)I_{k}(t).$

In addition, $I_{k}(t)\geq0$ for each $k$ that has been verified, as well as the fact that

Let

then when $\bar{R}_{k}^{0}\leq1,$ we get $\mu-\bar{S}_{k}^{0}\alpha_{k}b_{k}\geq0.$ On the other hands, $\mu-\bar{S}_{k}^{0}\phi_{k}(t)b_{k}\geq\mu-\bar{S}_{k}^{0}\alpha_{k}b_{k}.$ Thus when $\bar{R}_{k}^{0}\leq1,$ there is $\dot{L}_{1}(t)\leq0.$ And $\dot{L}_{1}(t) = 0$ if and only if $I_{k}(t) = 0, S_{k}(t) = \bar{S}_{k}^{0}.$ Which means that the disease-free equilibrium $E_{0}$ is the single point compact invariant set for the $\dot{L}_{1}(t) = 0.$ Therefore by the LaSalle's Invariant Principle, we can claim that the disease-free equilibrium $E_{0}$ is globally asymptotically stable.

From the other aspect, let

Then when $\underline{R}_{k}^{0} > 1,$ we have $-\mu+\bar{S}_{k}^{0}\alpha_{k}\Big(1-\frac{1}{N}\Big)b_{k} > 0,$ which implies that $-\mu+\bar{S}_{k}^{0}\phi_{k}(t)b_{k} > 0.$ It means that $\dot{L}_{1}(t) > 0,$ and the disease-free equilibrium $E_{0}$ is unstable in the neighborhood of $E_{0}.$ Combined with the boundedness of the system, there exists positive equilibrium.

From the system (2.1), we can easily obtain the nonzero solution, denote the nonzero solution as $E^{*} = \Big(S_{1}^{*}, I_{1}^{*}, S_{2}^{*}, I_{2}^{*}, \cdots, S_{m}^{*}, I_{m}^{*}\Big)$ that holds

where $\phi_{k}(t)|_{I_{k}^{*}} = \alpha_{k}$ due to the equilibrium state. In which the $S_{k}^{*}$ and $I_{k}^{*}$ equal to

Denote

Then when $\widetilde{R}_{k}^{0} > 1,$ e.g., $S_{k}^{0}-\frac{\mu}{\alpha_{k}b_{k}}-\sum_{l = 1, l\neq k}^{m}b_{lk}\frac{1}{\alpha_{k}b_{k}} > 0.$ On the other hand,

Therefore, $\widetilde{R}_{k}^{0} > 1$ implies that there exists the positive solution for the system (2.1). The proof of the theorem is then completed.

3.2.   The stability of positive equilibrium and behavior system's behavioral pattern synchronization

In this subsection, we will analyze the behavior pattern cluster synchronization induced by the patch disease under the case of positive equilibrium existence for epidemic systems, implying the epidemic disease spreading, and address the stability of positive equilibrium at the same time. For this purpose, we assume the relationship network is fixed, even though there are input and removing as well as the dispersal among patches, as for the input and death, the new individuals will fill up the removed one's relation site and remain the origin relationship unchanged. Assume there are $m$ patches, and we divide the relation network into $m$ behavior pattern relation groups induced by the patch disease information and let $\mathbb{K} = \{r_{k-1}+1, \cdots, r_{k}\}, k = 1, 2, \cdots, m$ represents the $k$th group with $r_{k}-r_{k-1}$ numbers that individuals belong to, where $r_{0} = 0, r_{m} = N.$ Specially, the individuals' numbers in the $k$th behavioral pattern relation groups do not need to be same with $k$th patch's population. Under this situation, we just focus on the behavior pattern of individuals, not the spatial position, that is to say the individuals from different places can have the same behavior pattern, so that we divide them into the same behavior pattern group. Corresponding to this division, when we divide the relationship network the individuals are also divided in fact and then utilize the patches disease information to design the feedback controllers for individuals in corresponding behavior pattern group, while this feedback control would be added to a part of the individuals in each group, named pinning.

Next, we divide the relationship network. Let

then we have the matrix $L$ as the relationship network $A$'s Laplacian matrix $L = (l_{ij})_{N\times N}$

Blocking $L$ according to the behavioral pattern relation groups, obtain

Additionally, combined with the intra-group enhancement factor, we have

where $\overline{L}_{kk} = c_{k}L_{kk}, k = 1, 2, \cdots, m,$ and $L_{kq_{k}} = L_{kq_{k}}, k = 1, 2, \cdots, m, q_{k} = N-(r_{k}-r_{k-1}).$ And the matrix block $\overline{L}_{kk}$ belongs to $A1$ and $L_{kq_{k}}$ belongs to $A3$ for each $k.$

Now, we use the patch disease information to design the feedback controllers $u_{i}(t)$ in the system (2.3) for individuals in each behavior pattern group to achieve behavior pattern grouped synchronization. And the controllers may add to apart of individuals in each behavior pattern group. We design the $u_{i}(t), i\in\mathbb{K}, k = 1, 2, \cdots, m$ as following

where $c_{\bar{i}}$ is the intra-group enhancement factor in the group that $i$ belongs, and $d_{i}$ represents the control strength with $d_{i} > 0$ for $i = r_{k-1}+1, \cdots, r_{k-1}+l,$ otherwise $d_{i} = 0.$ $e_{i}(t)$ has been denoted in the expression about $\phi_{k}(t)$ in subsection 2.1 called the behavior pattern grouped synchronization error for each individual $i$ in the group $k,$ and $y_{k}(t)$ is the behavior pattern target in the $k$th behavior pattern group that satisfies the uncoupled system

Then the whole error denotes as $e(t) = (e_{1}^{T}(t), e_{2}^{T}(t), \cdots, e_{N}^{T}(t))^{T},$ and $e^{T}(t)e(t) = \sum_{k = 1}^{m}\sum_{i\in\mathbb{K}}e_{i}^{T}(t)e_{i}(t).$ According to the behavior pattern error's definition, we get the patch epidemic-behavior error systems corresponding to the systems (2.3)

in which, from the $\overline{L}_{kk}\in A1$ and $L_{kq_{k}}\in A3$, there are

Then the problem about the stability of positive equilibrium for epidemic systems and the behavior pattern synchronization of the system (2.3) is transferred to the issue about the positive equilibrium stability and the zeros stability for the error system (3.9).

Therefore, we denote the patch epidemic-behavior cluster synchronization manifold $\mathcal{D}$ for patch epidemic-behavior cluster synchronization error network system (3.9) as shown

Then, we obtain the following result.

Theorem 2. Assume that the $B = (b_{kl})_{m\times m}$ is irreducible and the positive equilibrium exists for the epidemic system in (3.9). And assume that the $f_{k}\in QUAD$ for each $k.$ If there exists constant $\theta > 0,$ such that

where $\hat{c}_{k}$ represents the $k$th cofactor of the $k$th diagonal entry for the Laplacian matrix of $\hat{B} = (\hat{b}_{kl})_{m\times m},$ $\hat{b}_{kl} = b_{kl}S_{l}^{*},$ $\delta^{(k)} = \max_{j}\{\delta_{j}^{(k)}\}$ about $\Delta_{k} = \mathit{\text{diag}}\{\delta_{1}^{(k)}, \cdots, \delta_{n}^{(k)}\},$ and $\overline{D}_{k} = \mathit{\text{diag}}\{c_{k}d_{r_{k-1}+1}, \cdots, c_{k}d_{r_{k-1}+l}, 0, \cdots, 0\}.$ Then the patch epidemic-behavior cluster synchronization manifold $\mathcal{D}$ is globally asymptotically stable. In other words, patch epidemic-behavior synchronization network system (2.3) achieve behavior pattern cluster synchronization induced by the patch disease spreading, as well as the positive equilibrium $E^{*}$ for epidemic systems is unique and globally asymptotically stable.

Proof. Denote the Lyapunov function as

where

in which

Firstly, take the derivative of $V_{1}(t)$ along the behavior error equation, and obtain

Next, handling three terms in the above equation.

From the $QUAD$ class of function, we get

where $\Delta_{k} = \mathit{\text{diag}}\{\delta_{1}^{(k)}, \cdots, \delta_{n}^{(k)}\},$ $\delta^{(k)} = \max_{j}\{\delta_{j}^{(k)}\},$

$\widetilde{e}_{j}^{(k)}(t) = \Big(e_{r_{k-1}+1, j}(t), e_{r_{k-1}+2, j}(t), \cdots, e_{r_{k}, j}(t)\Big)^{T}\in\mathbb{R}^{r_{k}-r_{k-1}},$

$\widetilde{e}_{j}^{(q)}(t) = \Big(e_{1, j}(t), \cdots, e_{r_{k-1}, j}(t), e_{r_{k}+1, j}(t), \cdots, e_{N, j}(t)\Big)^{T}\in\mathbb{R}^{N-(r_{k}-r_{k-1})},$

that also be used in later analysis.

For the second term, we have

According to Lemma 1, we get

On the other hand,

Therefore, substitute the Eqs (3.17) and (3.18) into Eq (3.16), we obtain

where $E_{r_{k}-r_{k-1}}$ is identity matrix with dimension $r_{k}-r_{k-1}.$

From the assumption that there exists positive equilibrium for epidemic systems, then there holds limitation

which means that for any $\varepsilon > 0$ there is $t_{0} > 0,$ when $t > t_{0}$ there holds

e.g.,

So there are

where

In conclusion, substitute the Eqs (3.15), (3.19) and (3.20) into the Eq (3.14), we obtain

Then taking the derivative of the $V_{1}^{(k)}(S_{k}(t))$ and $V_{2}^{(k)}(I_{k}(t))$ along with the epidemic model in system (2.3), respectively, and according to Eq (3.6) that positive equilibrium satisfies, we obtain

In addition, there exist the fact that

and the equality holds if and only if $S_{k}(t) = S_{k}^{*}.$

Then, from the expression of $V_{2}(t),$ there are

Be similar to Theorem 1, let $\hat{b}_{kl} = b_{kl}S_{l}^{*},$ and $\hat{c}_{k}$ represents the $k$th cofactor of the $k$th diagonal entry for the Laplacian matrix of $\hat{B} = (\hat{b}_{kl})_{m\times m}.$

Then denote

We get

where

On the other hand,

Thus, according to the result in ^[39], one obtain

From above, we get

In conclusion, combined Eqs (3.14) and (3.25) as well as the $V(t)$'s denote in Eq (3.11), there are

Then, we obtain when Eq (3.10) holds, the $\dot{V}(t) < 0.$ Therefore, from the Lyapunov stability theorem, there are results that when positive equilibrium exists for the epidemic-behavior synchronization model (2.3), it is unique and globally asymptotically stable, meanwhile behavior systems achieve the behavior pattern cluster synchronization induced by patch diseases. In other works, under the conditions of positive equilibrium existence for epidemic systems, combined with conditions (3.10), then the patch epidemic-behavior cluster synchronization manifold $\mathcal{D}$ for the patch epidemic-behavior synchronization network is globally asymptotically stable. Then the theorem is completed.

To achieve the conditions (3.10) in Theorem 2, we can carry out from the following steps.

Step 1. Firstly, we obtain the largest eigenvalue of the $\Big(\frac{-\eta_{k}}{\lambda_{\max}(P)}+\delta^{(k)}+\frac{m-1}{2\theta}+\frac{m}{N}\frac{\hat{c}_{k}\alpha_{k}b_{k}}{\lambda_{\min}(P)}\Big)E_{r_{k}-r_{k-1}}+\frac{\theta}{2}L_{kq_{k}}L^{T}_{kq_{k}},$ denote as $\rho_{k}.$

Step 2. Then, we also get the least eigenvalue of $L_{kk}+D_{k}(I_{k}^{*}-\varepsilon)$ when given the strength of controlling, denote as $w_{k},$ where

Step 3. Naturally, in order to making condition (3.10) be satisfied, it is needed to select the appropriate $c_{k},$ which implies that the $c_{k}$ is just selected larger than $\frac{\rho_{k}}{w_{k}},$ i.e., $c_{k} > \frac{\rho_{k}}{w_{k}}.$

From the steps we conclude that, on the one hand, the internal group enhance factors $c_{k}$ are working as an important role in synchronizing the same group individuals, in other words, behavior cluster synchronization and stability of positive equilibrium positive can achieve needed to select the $c_{k}$ as large as possible. Which means that we can strengthen subetaoup connection to realize the desired target. On the other hand, it is unrealistic to choose any large $c_{k}$ due to the limited resource, therefore this procedure also provides a feasible method to select the advisable $c_{k}$ working among the internal group individuals.

Remark 4. Notice that from the $QUAD$ class of function in the Definition 3 the conditions (3.10) can also be represented as

which is convenient to compute in the numerical simulation. Furthermore, it is obvious that conditions (3.10) are relative to the relation sub-group network's structure, susceptible population's dispersal rates, and the scale of patchy epidemic diseases. In other words, in view of cluster conditions (3.10), the behavior pattern cluster synchronization is codetermined by topological structure of individuals' relation network and the migration for the susceptible population as well as the appropriate external controlling based on the epidemic disease spreading.

Remark 5. In conditions (3.10), $\overline{L}_{kk}$ is inner-group relation network structure containing relation enhance factor, that is to say the relation sub-networks among nodes in the same group are weighted networks. Thus, combined with appropriate the strength of controlling, it is easily to achieve conditions (3.10). What's more, from the designed control strategy shown in Eq (3.7), the external control can work well only when epidemic disease spreading, thus, conditions (3.10) will be satisfied when giving appropriate strong of control.

4.   Numerical simulations

In this section, we give the numerical simulation to demonstrate the effectiveness of our results obtained in the previous sections.

As for the local dynamics in the behavior system, we consider the unified system as shown in the early work ^[40] and its generation ^[41]

From the reference ^[40], we know that for all $\tau\in[0, 1],$ the unified systems are all chaotic. And for $0\leq \tau < \frac{1}{29},$ there is the positive invariant set for the above system shown in ^[40],

Then,

We denote that $f_{1}(x, \tau) = (25\tau+10)(x_{2}-x_{1}), f_{2}(x, \tau) = (28-35\tau)x_{1}-x_{1}x_{3}-(29\tau-1)x_{2}, f_{3}(x, \tau) = x_{1}x_{2}-\frac{8+\tau}{3},$ then $f(x, \tau) = (f_{1}(x, \tau), f_{2}(x, \tau), f_{3}(x, \tau))^{T}.$ Let $e = x-y,$ when $p_{2} = p_{3},$ then there are

where $\xi^{(\tau)}_{1} = -(25\tau+10)p_{1}+\epsilon \frac{(25\tau+10)p_{1}+(M+56-70\tau)p_{2}}{2}+v\frac{p_{2}M}{2},$ $\xi^{(\tau)}_{2} = (29\tau-1)p_{2}+\frac{(25\tau+10)p_{1}+(M+56-70\tau)p_{2}}{2\epsilon},$ $\xi^{(\tau)}_{3} = -\frac{8+\tau}{3}p_{2}+\frac{p_{2}M}{2v}.$

In this numeric simulation, we set $\tau = 0, 0.01, 0.03$ to get three chaotic systems as the local dynamical behavior in the behavior dynamic systems. And from the $QUAD$ class of function, we just need to obtain the positive value $\xi_{k}$ for $k = 1, 2, 3.$ Thus, we set $P = diag\{0.099, 0.1, 0.1\}$ and $\epsilon = 1,$ $v = 1.$ Then we get these values

Then we let $\xi_{1} = \max_{k = 1, 2, 3}\{\xi^{(0)}_{k}\} = 5.1923,$ $\xi_{2} = \max_{k = 1, 2, 3}\{\xi^{(0.01)}_{k}\} = 5.3821,$ $\xi_{3} = \max_{k = 1, 2, 3}\{\xi^{(0.03)}_{k}\} = 8.4266.$

Therefore, we obtain the following relationship for every local behavior dynamics when $k = 1, 2, 3$ corresponding to the Definition 3

for any $x(t), y(t)\in\mathbb{R}^{n}.$

We consider 18 nodes network and each subetaoup contains six nodes, and the local behavior dynamics for each group are shown in above, then we denote the relationship network structure's Laplacian matrix corresponding to the relationship network is defined as

where {{

}} {{

}} Then we select the intra-group enhancement factors $c_{1} = 32, c_{2} = 35, c_{3} = 50.$

Then, we select the parameters in the epidemic system, and let the natural death rate as $\mu = 0.5,$ the awareness coefficients for each group denote as $\alpha_{1} = 0.2, \alpha_{2} = 0.15, \alpha_{3} = 0.2,$ and the infective rates $b_{1} = 0.06, b_{2} = 0.05, b_{3} = 0.06.$ And the dispersal matrix for the susceptible individuals denote as

Select $S^{0}_{1} = 5, S^{0}_{2} = 7, S^{0}_{3} = 6.$ Then we can verify the condition $(H1)$ holds as

Then, the disease-free equilibrium can be derived as

so that one can get

thus they are all less than $1,$ so $E_{0}$ is globally asymptotically stable from the Theorem 1. When select any 4 different initial values for the epidemic systems, then the state graphs for the epidemic system under the different initial values are shown in Figure 2, in which we use different colors to reveal the different state variables for susceptible and infective individuals, that is to say, the same color curves represent the same state variable curves with different initial values. In Figure 2, the infective states curves all trend to the zeros states and the susceptible states curve trend to positive stable states, which means that the states of susceptible and the infective all trend to the stable one, i.e., the disease-free equilibrium $E_{0}.$ Meanwhile, the behavior states are not synchronization under this time due to the controllers unworked ultimately, and the errors between states and target are shown in Figure 3, in which neither the intra-group behaviors nor between-group behaviors trend to synchronization, i.e., the corresponding errors curves toward nonzero, where the intra-group errors and the between-group errors are defined as $E_{k} = \sqrt{\sum_{i, j\in\mathbb{K}}\|x_{i}(t)-x_{j}(t)\|^{2}}, k = 1, 2, 3,$ $E_{12} = \sqrt{\sum_{i = 1}^{6}\sum_{j = 7}^{12}\|x_{i}(t)-x_{j}(t)\|^{2}},$ $E_{13} = \sqrt{\sum_{i = 1}^{6}\sum_{j = 13}^{18}\|x_{i}(t)-x_{j}(t)\|^{2}},$ $E_{23} = \sqrt{\sum_{i = 7}^{12}\sum_{j = 13}^{18}\|x_{i}(t)-x_{j}(t)\|^{2}},$ respectively.

When we reset the $\alpha_{k},$ $\mu$ and $b_{k}, k = 1, 2, 3,$ such as $\alpha_{1} = 0.2, \alpha_{2} = 0.1, \alpha_{3} = 0.15,$ $\mu = 0.4,$ $b_{1} = 0.6,$ $b_{2} = 0.8,$ $b_{3} = 0.7.$ In this situation, the $(H1)$ also holds. And we can compute that

which means that $\underline{R}^{0}_{k} > 1$ for $k = 1, 2, 3.$ Then from the Theorem 1, one can claim that the disease-free equilibrium is unstable. Further, one obtain that $\widetilde{R}^{0}_{k} = \frac{S^{0}_{k}\alpha_{k}b_{k}}{\mu+\sum_{j\neq k}b_{jk}} > 1, k = 1, 2, 3,$ where

Therefore, according to the Theorem 1, there are positive equilibrium, and $E^{*} = (S^{*}_{1}, I^{*}_{1}, S^{*}_{2}, I^{*}_{2}, S^{*}_{3}, I^{*}_{3}) = (3.3333, 1.756, 5, 2.9286, 3.8095, 1.1726).$

Select the $\theta = 2,$ and only added controllers to first three nodes in each subetaoup, and the control strength for each group's nodes are selected as

In addition, the cofactor can also be computed as $\hat{c}_{1} = 0.0083, \hat{c}_{2} = 0.0101, \hat{c}_{3} = 0.0871.$

Then, according to the Remark 4, the conditions in the Theorem 2 can be verified easily, and when $k = 1, 2, 3,$ there are

With these parameters, the epidemic states' curves can be shown as in Figure 4, meanwhile, the behavior errors curves are also shown in Figure 5. In this situation, we also use 4 different initial values for the epidemic systems as given in the case of disease-free, and as shown in Figure 4, one can see that the susceptible states' curves and infective states' curves under different initial values all trend to the same stable nonzero constant states with the time goes on, ultimately. In other words, the positive equilibrium is globally asymptotically stable. As well as the behavior states errors in Figure 5 can achieve the desired states, which means that the behavior errors in the same subetaoup achieve the behavior synchronization, such as the first subfigure in Figure 5 with each subetaoup's behavior errors states curve trending to zeros, ultimately, and the other three subfigures representing the errors among different subetaoup do not trend to zeros. It is to say that the nodes intra-group achieve behavior synchronization, while the between's ones are not.

Besides, the simple comparison about the effectiveness of behavior awareness to epidemic states is also shown. Under the situation of positive equilibrium existing, and the parameters are same as in corresponding cases, we compare three cases about the infected term in each patch, such as $b_{k}S_{k}(t)I_{k}(t),$ $\alpha_{k}b_{k}S_{k}(t)I_{k}(t)$ and $\phi_{k}(t)b_{k}S_{k}(t)I_{k}(t),$ for $k = 1, 2, 3,$ respectively. And plotting their states' curves as shown in Figure 6. There are different marked states' curves corresponding to represent distinct cases, and one can see the stable infected states with cofactor $\alpha$ and $\phi$ before incidence are lower than the without cases, i.e., the cofactor is $1.$ And comparing the effect among $\alpha$ and $\phi,$ we see that their curves are almost the same, but the susceptible curves with $\phi$ is upper than the $\alpha$ curves and the infected curves with $\phi$ are lower than the one with $\alpha.$ Which means that there exists tiny difference in effect between $\alpha$ and $\phi$ that can not ignore. Those two cases trend to the same stable states, ultimately. This also shows that the behavior of individuals can impact the process of epidemic.

Remark 6. Notice that from the range of $\alpha_{k}$ and $\phi_{k}$, we get $b_{k}S_{k}(t)I_{k}(t) > \alpha_{k}b_{k}S_{k}(t)I_{k}(t) > \phi_{k}(t)b_{k}S_{k}(t)I_{k}(t),$ which means the numbers of infected individuals in the unit time under the later two cases are less than the first case. From Figure 6, and under the influence about parameters $\alpha_{k}$ and $\phi_{k},$ the infected number level are becoming lower. Furthermore, when errors of behavior trend to zero, then $\phi_{k}$ trend to $\alpha_{k},$ so, their ultimate states will be same, while the process may not the same due to the $\phi_{k}$ contains behavior pattern.

Remark 7. To be notice that in our simulation part, we just list some parameters values in systems to demonstrate our main results and it is not based on the practical data. In the future work, there will be a trying direction to introduce the real epidemic disease data to simulation and even forecast the disease trend, such as the COVID-19 spread. And structure about the real relationship network is a vast work that will be our aim in the next work.

5.   Conclusions

In this paper, we investigate the problem about behavior pattern cluster synchronization aroused by patchy epidemic transmission based on external feedback pinning control strategy, in which the coupled models combined with behavior synchronization systems and patchy dispersal epidemic systems with constant population are constructed, where so called contact willing functions are introduced in the epidemic systems and the dispersal is only concerned among the susceptible. Firstly, in view of the constant population the boundedness is obvious, and the nonnegativity is also discussed by the qualitative analysis for epidemic systems. And then we concern the stability of disease-free equilibrium for the patchy epidemic systems and obtain the patchy related conditions that guarantee globally asymptotically stable by using the Lyapunov analysis method. As well as, conditions that insure the existence of positive equilibrium for epidemic systems are also derived. Next, behavior pattern cluster synchronization aroused by patchy epidemic disease under the case of disease spreading, as well as the global stability of positive equilibrium for epidemic systems are addressed at the same time by means of Lyapunov function and the Lyapunov stable theory, meanwhile, sufficient criteria, related to the structure of inner-group relation, dispersal matrix of the susceptible and scale of diseases, making individuals to achieve behavior pattern cluster synchronization, are derived. The simulation example is exhibited to demonstrate the effectiveness of the main results obtained here.

Acknowledgments

This work was jointly supported by the NSFC grants under Grant Nos. 11572181 and 61873171.

Conflict of Interests

All authors declare no conflicts of interest in this paper.