Turing Patterns for a Nonlocal Lotka–Volterra Cooperative System

This paper is devoted to investigating the pattern dynamics of Lotka–Volterra cooperative system with nonlocal effect and finding some new phenomena. Firstly, by discussing the Turing bifurcation, we build the conditions of Turing instability, which indicates the emergence of Turing patterns in this system. Then, by using multiple scale analysis, we obtain the amplitude equations about different Turing patterns. Furthermore, all possible pattern structures of the model are obtained through some transformation and stability analysis. Finally, two new patterns of the system are given by numerical simulation.


Introduction
Pattern is a kind of non-uniform macrostructure with certain regularity in space or time, which exists widely in nature. All kinds of pattern structures constitute a colorful and charming world. Therefore, it is of great significance to understand the cause and mechanism of pattern formation for revealing the mystery of nature formation. In order to understand its formation mechanism, many scholars turned their work to study various mathematical models with different factors, such as delay [14,16,36], noise [19], network [25,28,34], season [29] and so on, among which, the reactiondiffusion equation was often used in [1,22] and the references therein.
In particular, nonlocal effect has attracted more and more attention because it can describe a variety of morphogenetic phenomena. Then, how does the nonlocal effect play an important role in the formation of patterns? Actually, the research on the nonlocal model can be traced back to Britton [5,6]. He found that introducing nonlocal effects into a single reaction diffusion equation would lead to the change of stability about the positive equilibrium point and the disturbance around it (it should be noted that when considering spatiotemporal patterns, the instability of the equilibrium point and the appearance of branches are often accompanied). Subsequently, scholars further considered this problem in a limited region, and studied the Neumann boundary condition [10] and Dirichlet boundary condition [27]. Recently, Guo [11,12] and Yang and Xu [32] investigated the spatiotemporal patterns for a single species reaction diffusion equation. Of course, the introduction of nonlocal has great influence not only on the dynamics of patterns, but also on the traveling wave solution [3,4,17] and the well posedness about the solution of Cauchy problem [15,16,33].
The above considerations are the influence of nonlocal effect on the single population reaction-diffusion model, so a very natural question, how does it affect the pattern dynamics of reaction-diffusion system? In this paper, we consider the influence of nonlocal effect on the pattern dynamics of Lotka-Volterra cooperative system. Specifically, we will study the following model where and for (x, t) ∈ ℝ 2 × (0, ∞) . Here, u and v indicate the population density of two species. d an r denote the diffusion coefficient and the natural growth rate of species v respectively. And a 1 , a 2 represent the intensity of interaction between two species. In addition, the coefficients d, r, a 1 , a 2 , and the population density are both positive. For more biological interpretation about model (1.1), we can refer to [20,23,24].
Obviously, system (1.1) has at least three equilibrium points (0, 0), (0, 1), (1,0). And there is a fourth equilibrium point when a 1 a 2 < 1 . The appearance of the fourth equilibrium point means that two groups finally will reach the state of coexistence, so from the biological point of view, we always assume that a 1 a 2 < 1 . In addition, we know that the classic Lotka-Volterra cooperative system corresponding to system (1.1) will not have Turing patterns under normal conditions. Under the influence of Predator-Prey and Lotka-Volterra competition system, we will study whether and what kind of patterns will appear in system (1.1) after introducing nonlocal effect? In fact, the study of nonlocal effect on the pattern dynamics of Predator-Prey and Lotka-Volterra competition systems is a hot topic in recent years. Xu et al. [30] considered the following Predator-Prey system and gave the possible Turing patterns about system (1.3) by using stability analysis and numerical simulation. Following, Chen and Yu [8] researched the model with nonlocal intraspecific competition of prey for resources, and found the steady state can lose the stability when conversion rate passes through some Hopf bifurcation and the bifurcating periodic solutions near such bifurcation value can be spatially nonhomegeneous. It should be pointed out that this is very different from the classic results. Guo and Yan [13] investigated the following Lotka-Volterra competition system with nonlocal delay effect where By means of Lyapunov-Schmidt reduction, the normal form theory and the center manifold reduction, the existence and stability of spatially nonhomogeneous steadystate solutions and bifurcation direction of Hopf bifurcating periodic orbits are given. Then, Yan and Guo [31] further extended the above conclusion to Lotka-Volterra model with cross-diffusion and delay effect. And Han and Wang [14] also studied the following Lotka-Volterra competition diffusion system with nonlocal term where * * u and * * v are similar to model (1.1) with (1.2). Through using the multi-scale analysis and numerical simulation, they obtained the pattern of system (1.4) under what conditions and gave the verification. For more research about the influence of nonlocal (or delay) effects on the system, see [2,7,17,18,21,35].
However, we find that few studies consider the influence of nonlocal effect on Lotka-Volterra cooperative system. It may be that in a general sense, people think that branch is more likely to occur in predator-prey or competing relationship, which leads to the appearance of Turing patterns. Of course, from a mathematical point of view, it is easier to deal with the first two cases than the cooperative 1 3 system. In this paper, we will try to study the pattern dynamics of model (1.1). Specifically, we always assume a 1 a 2 < 1 in this paper. Firstly, by analyzing the characteristic equation of the corresponding system, the conditions of Turing bifurcation are obtained. Further, according to these conditions, we give Turing space which is made up of diffusion coefficient d and time delay by numerical simulation. Secondly, through multi-scale analysis, we get the amplitude equations, and moreover with the help of some transformation and solving steadystate solution, we obtain all possible pattern structures of system (1.1). Finally, through numerical simulation, we give two different patterns of system (1.1).
The remainder of the paper is arranged as follows. In Sect. 2, we discuss the conditions of Turing bifurcation caused by nonlocal effect in Lotka-Volterra cooperative system, and give Turing space by numerical simulation. Then we use multi-scale analysis, stability analysis and so on to give the possible pattern of system (1.1), and through numerical simulation , two kinds of specific pattern forms in Sect. 3 are given. Finally, a brief discussion is drawn in Sect. 4.

Turing Bifurcation
In this section, we research the Turing bifurcation of (1.1) and seek out the corresponding Turing space. It's worth noting that here we mainly consider the influence of nonlocal effect. Firstly, we prove that the equilibrium (u * , v * ) of the classic Lotka-Volterra cooperative system is stable. Then, we show that (u * , v * ) will become unstable after the introduction of nonlocal effect. Lastly, by analyzing the characteristic equations of diffusion and non diffusion, some conditions for Turing bifurcation occurring near the positive equilibrium point (u * , v * ) are given.
By linearizing the classic Lotka-Volterra cooperative system near (u * , v * ) , we can get Select the following test functions where i is the imaginary number unit, is the growth rate of the disturbance in time t, ⋅ = k 2 and k is the wave number, = (X, Y) is the two-dimensional space vector. By substituting Eq. (2.2) into Eq. (2.1), we can get the characteristic equation which is equivalent to It is obviously that the characteristic Eq. (2.3) has two negative characteristic roots, so the system is stable near the equilibrium point (u * , v * ) . Next, we consider the influence of nonlocal effect. In order to simplify our research system, let w(x, t) = ( * * u)(x, t) , the system (1.1) is able to be converted to the following form Similar to the above discussion, the system (2.4) has three equilibrium points (0, 0, 0), (0, 1, 0), (1, 0, 1). And when a 1 a 2 < 1 , the system (2.4) has the fourth equilibrium point Next, linearize the system (2.4) near the equilibrium point (u * , v * , w * ) , and we can obtain where Select the following test functions (2.5) where i is the imaginary number unit, is the growth rate of the disturbance in time t, ⋅ = k and k is the wave number, = (X, Y) is the two-dimensional space vector. By substituting Eq. (2.6) into equation (2.5), we can get the characteristic equation which is equivalent to where Obviously, from the characteristic Eq. (2.7), we can see that the eigenvalue is not always negative, which further shows that under the influence of nonlocal delay, the system (2.4) may lose stability near the positive equilibrium point (u * , v * , w * ). Next, with the help of analysis the characteristic equation for non-diffusion and diffusion, we look for some conditions of Turing bifurcation arising around (u * , v * , w * ).
When the system (2.4) has no diffusion term, its characteristic equation is as follows where So, we can get the proposition below.

Proposition 1
The positive equilibrium point (u * , v * , w * ) of system (2.4) is stable only when these conditions are satisfied: Proof According to the ordinary differential equations, we know that the necessary and sufficient conditions for the stability of equilibrium point (u * , v * , w * ) are: 2,3). Then, according to the Routh-Hurwitz criterion, we can get Re( i ) < 0 only when the above conditions are satisfied. This completes the proof. ◻ When the system (2.4) has diffusion term, similar to Proposition 1, we can obtain that the positive equilibrium point (u * , v * , w * ) of (2.4) is stable only when these conditions are satisfied: When any of these conditions is broken, the coexistence equilibrium state of (2.4) will become unstable. In the next step, we will analyze the characteristic Eq. (2.7) to obtain the necessary condition for Turing bifurcation. Under this condition, the system will appear spatiotemporal pattern. Therefore, we consider the following three cases: Let b 3 (k) = F(k 2 ) and z = k 2 , we can obtain where Obviously, we can observe that f 1 > 0 , f 2 > 0 , and therefore, Turing bifurcation can only appear under the condition of f 3 < 0 . Furthermore, we can know the following properties about F(z).
dz is a quadratic function, and Δ > 0 , so F(z) has two poles: If z 1 , z 2 are the poles, we can see that the maximum value of the equation is keeps to the left of the minimum value, that is, z max < z min , so z min = z 1 .
According to the previous discussion and analysis, we can obtain that when b 3 (k) ≤ 0 and the condition (2.8) is true, Turing bifurcation will appear in system (2.4). Therefore, we can obtain the conclusion below. Conclusion 1 If the condition (2.8) is true, Turing bifurcation will be produced only when the following conditions are established.
Proof From the above analysis, f 3 < 0 is necessary. Then, from the properties of F(z), we can easily get z 1 > 0 , and z 1 is the minimum point (for the convenience of understanding, we drew the sketches of dF(z) dz and F(z), as shown in Fig. 1). Therefore, where Obviously we can get that h 1 > 0 , h 2 > 0 . Similar to the discussion of F(z), we can know the following properties about H(z).
is an increasing function. If Δ > 0 , then H(z) has two poles: Fig. 1 This graph shows the general trend of dF(z) dz and F(z) 1 3 If z 1 , z 2 are the poles, we can see that the maximum value of the equation is keeps to the left of the minimum value, that is, z max < z min , so z min = z 1 .
According to the previous discussion and analysis, we can obtain that when b 1 (k)b 2 (k) ≤ b 3 (k) and the condition (2.8) is true, Turing bifurcation will appear in the system. Therefore, we can obtain the conclusion below. Conclusion 2 If the condition (2.8) is true, Turing bifurcation will be produced only when the following conditions are true Proof It can be seen from the condition (2.8) that h 1 > 0, h 2 > 0, h 4 > 0 , therefore,Turing bifurcation can only appear under the condition of h 3 < 0 . So h 3 < 0 is necessary. Then, from the properties of H(z), we can easily get z 1 > 0 , and z 1 is the minimum point (for the convenience of understanding, we drew the sketches of dH(z) dz and H(z), as shown in Fig. 2). Therefore,the condition of H z 1 ≤ 0 is necessary. ◻ Obviously, we can see that the conditions for conclusions 1 and 2 are a little complex. For the convenience of understanding, we have plotted Fig. 3, which directly shows the overview of these conditions. The green area indicate the value range of the parameters that may appear Turing bifurcation. In other words, if we set the values of r , a 1 , a 2 to 0.1, 0.8 and 0.9, and use the conditions in conclusion 1 and conclusion 2, we can obtain the Turing space composed of parameters d and , as shown in Fig. 3.
It should be noted that w(x, t) satisfies the third formula in system (2.4) when (x) with form (1.2). In this case, the system (1.1) and the system (2.4) are completely equivalent, and the local stability about the equilibrium point is completely equivalent by using the characteristic method, that is to say, the stability of (u * , v * , w * ) in system (2.4) is same with the stability of (u * , v * ) in system (1.2). This is different from the stability of traveling wave in article [20]. Furthermore, even if we consider the global stability (including the influence of initial value), since the initial of w(x, t) in system (2.4) can be calculated by the initial of u(x, t) in system (2.4) (see P1977 in [17] or P16 in [16]), system (1.1) and system (2.4) both only need the initial values of u and v. Therefore, to a certain extent, we can use system (2.4) to study system (1.1), so as to overcome the difficulties caused by nonlocal effects.

Multiple Scale Analysis and Conclusions
In this section, we plan to make use of the multiple scale analysis method to acquire the amplitude equations that form the different Turing patterns. Then, by analyzing the stability of the amplitude equations, we can obtain the sufficient conditions for the appearance of different Turing patterns, so that different Turing patterns can be constructed by numerical simulation.
Since the coexistence equilibrium will not be stable until the wave number disturbance approaches the critical value k = k T , the multiple scale analysis method is mostly used near the bifurcation threshold. Therefore, the multiple scale analysis method is widely applied and becomes an effective method to obtain the amplitude equations. Specifically, we will first expand the nonlinear term N, variable U and control parameters with a small parameter . Then we can substitute them into the equation and compare the coefficients of , 2 , 3 to obtain three new equations. And then using Fredholm solvability conditions, we will acquire the amplitude equations. Finally, we will explore the stability of amplitude equations in order to construct different Turing patterns.
In the next, we first look for the amplitude equations of the system (1.1). Because in the vicinity of = T , the eigenvalues are close to zero near the critical modules, we can know that they are slowly changing modules. However, those modules that deviate from the critical mode will soon relax. Therefore, we are supposed to consider the perturbations of k in the vicinity of k T . To do this, we first look for k T . Therefore, we are supposed to find the extreme point of b 3 (k) = F(k 2 ) for k 2 , that is, we are supposed to obtain the solution of dF(z) dz = 3f 1 z 2 + 2f 2 z + f 3 . From the discussion in section 2, we can get: Substituting k = k T into b 3 (k) = 0 , the value of T can be obtained. It should be noted that k T and T do not need to give a specific formulas. We only need to give the values of parameters a 1 , a 2 , d and r to get the values of k T and T . Next, we rewrite the system (2.4) near the equilibrium point (u * , v * , w * ) to the following form where Since close to onset = T , the solutions of system (2.4) are able to be represented as are a couple of oscillatory wave vectors with different directions. Analogously, the solutions of system (3.1) are able to expand to (3.2) as follows where U s indicates the homogeneous steady state and U 0 = (l 1 , l 2 , 1) T (the value of l 1 , l 2 will be determined below) is the characteristic value of the linearized operator.
U 0 A j e ik j ⋅r +Ā j e −ik j ⋅r , A j and the conjugate Ā j are the amplitudes integrated with the modes k j and −k j , respectively. Let U = (u, v, w) T , N = (N 1 , N 2 , N 3 Using the multiple scale method, we can expand the nonlinear term N, the variable U and the control parameter in the following forms with a small parameter . where h 2 and h 3 express the second and the third order of in the expansion of the nonlinear term N, respectively, and  For the first order of as (3.12), because L T is the linear operator close to the onset in the system, (u 1 , v 1 , w 1 ) T is the linear combination of characteristic vector corresponding to the characteristic value 0. Solving (3.12), we can get where W j represents the amplitude about the mode e ik 1 r (j = 1, 2, 3) of the system on the first order perturbation, and c.c. represents the conjugate of the previous term. Its form is determined by the perturbation terms of the higher order.
For the second order of as (3.13), we have In order to make sure the existence of the nontrivial solution of this Eq. (3.16), let the adjoint operator of L T be L * T . And then according to the Fredholm solubility conditions, we can know that the vector-valued function on the right hand side of (3.16) need be orthogonal with the zero characteristic vector of operator L * T . In this system, the zero eigenvector of L * T is where Therefore, we can obtain where F j u , F j v and F j w correspond to the coefficients of e ik j r in F u , F v and F w , respectively, it means that According to (3.15) and (3.16), we can get (3.17) Using the Fredholm solubility condition and taking (3.18)-(3.20) into (3.17), we can get the following results By substituting (3.15) into (3.16), we can get At the same time, we can obtain the coefficients in (3.22) through solving the linear equations about e 0 , e ik j r , e i2k j r , e i(k j −k m )r , and we can know where For the third order of as (3.14), we have

3
They exist under the following conditions: .
(d) A mixed structure solution with g 2 > g 1 . Then several different Turing patterns of Eq. (1.1) are given by numerical simulation. These simulation results reflect the spatiotemporal distribution of these studied species, and the associated spatiotemporal patterns are shown here. For this, we consider the Neumann boundary conditions. We set the values of parameter as follows: a 1 = 0.8, a 2 = 0.9, r = 0.1 and take the size of the space area as 1000 × 1000 , space steps as Δx = 1, Δy = 1 and time steps as Δt = 0.05 . In addition, we consider that From the above numerical simulation, we find that system (1.1) will appear strip pattern (see Figs. 4,5) and spot pattern (see Figs. 6, 7) under certain conditions. In addition, these patterns also show that in real life, the distribution of species is uneven and affected by many factors. In particular, in this paper we consider the combined effect of the diffusion coefficient of species v and time delay. From Figs. 4 and 5, we find that the system (1.1) will produce strip pattern under the combined effect of d and , that is to say, in some areas, the population density over (or falls below) the maximum (or minimum) value. This phenomenon is similar to the distribution of forest after desertification control (see [26]). From Figs. 6 and 7, we find that the system (1.1) finally appeared spot pattern. In the spots pattern (see [9]), we find that the population density about species u and v may has maximum (or minimum) values in some places. Obviously, because we only consider the effect of d and , the results obtained are often simpler than those in real nature.

Discussion
In this paper, we systematically study the effect of diffusion coefficient d and time delay on the pattern of nonlocal Lotka-Volterra cooperative system (1.1). It should be noted that in order to simplify the calculation and learn more about the influence of nonlocal effects on Lotka-Volterra cooperative system, we only consider nonlocal effects on species u here. In fact, there are nonlocal effects on species u and species v, and the research ideas are basically similar.
In addition, the study of this paper breaks our inherent thinking. Generally, we think that the pattern is generated in the competition or predation process, but through the study of this paper, we found that pattern can also appear in the nonlocal Lotka-Volterra system. The research in this paper enriches the theory of pattern dynamics to some extent. We know that many natural phenomena can be explained by using the theory of pattern dynamics. Here we just make a preliminary attempt. In the future, we will further study this aspect and expect to use our theory to explain more natural phenomena.