Turing Patterns in a Competition Lotka-Volterra Model with Nonlinear-Diffusion

1. Introduction

In the subsequent decades, researchers developed and refined various reaction-diffusion models, exploring their applications in different fields, especially in Ecology ([2,6,18]). The Lotka- voltera model is defined by Lotka (1925) in [27] and Volterra (1926) in [28] describes the interaction of two species in competition models for the same resources and is mathematically explained with a negative sign. The Lotka-Volterra competition model with linear diffusion term has been studied widely ([12,13,20,25,31,32,33,34]). Besides that, the dynamics of predator-prey interactions were also interactive for studying some complex phenomena. The integration of nonlinear diffusion into these models added a new layer of complexity and realism to pattern formation studies. The foundation for reaction-diffusion models can be traced back to the work of Alan Turing in 1952 ([4]). Turing’s seminal paper, “The Chemical Basis of Morphogenesis,” demonstrated how reaction-diffusion systems could generate complex patterns. This work laid the groundwork for the study of pattern formation in nonlinear systems. In [30], with a typical Lotka–Volterra–Gause reaction-diffusion model and constant diffusion coefficient the producing of a pattern was not achieved, however, the existence of non-linear diffusion in Shigesada-Kawasaki-Teramoto (SKT) model was enough to build patterns and establish Turing instability. In the context of ecology, non-linear diffusion manifests as either cross-diffusion or self-diffusion. Cross-diffusion refers to the phenomenon where the movement of one species affects the diffusion of another species, whereas self-diffusion describes the spread of individuals within the same species. Both theoretical models and experimental studies have demonstrated that cross-diffusion significantly contributes to the formation of spatial patterns [1,3,11]. Specifically, cross-diffusion mechanisms are crucial in explaining how species interactions can lead to complex pattern formation, such as stripes, spots, or labyrinthine structures in population distributions. This model extends the focus on pattern formations as in SKT model ecology system as in [25] and to bridges the ecology with solid state system by adapting nonlinear diffusion effect in crystal growth. One of the pivotal concepts in understanding pattern formation through diffusion processes is Turing instability. This instability arises when the diffusion rates of interacting species lead to the spontaneous emergence of spatial patterns. To theoretically verify the presence of Turing instability in a given ecological model, four specific conditions must be satisfied. These conditions, detailed in the works of Barrio et al. (1999), Othmer et al. (2009), and Page et al. (2005) provide a rigorous framework for predicting and analyzing pattern formation due to non-linear diffusion [8,9,10]. The system we will study is

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial u_1}{\partial t}}& =& \nabla .J_1+u_1(1+a_1{u}_1-(1+a_1)u_1^2-b_1{u}_2),\hfill \\ \hfill {\displaystyle \frac{\partial u_2}{\partial t}}& =& \nabla .J_2+u_2(1+a_2{u}_2-(1+a_2)u_2^2-b_2{u}_1).\hfill \end{array}$$

(1)

where

$$\begin{array}{c}\hfill J_1=\nabla ({\alpha }_1{u}_1+{\beta }_1{u}_2)u_1,\\ \hfill J_2=\nabla ({\alpha }_2{u}_1+{\beta }_2{u}_2)u_2\end{array}$$

(2)

$u_1$ and $u_1$ are the interaction species as a function of x and t. The diffusion coefficients are nonlinear, while the reaction term includes logistic growth of second order $u_1(1-(1+a_1)u_1^2)$ and $u_2(1-(1+a_2)u_2^2)$, the cooperation within species are defined by the terms ${\alpha }_1{u}_1$ and ${\alpha }_2{u}_2$, while the competition terms are $b_1{u}_2$ and $b_2{u}_1$. All parameters are positive. Note that this model is studied before in [25] in the case of constant diffusion coefficients. Pattern formation is absent for constant diffusion in [25] and the reaction terms with different parameter values was not sufficient to produce any type of pattern solutions. Therefore, we extend the study to nonlinear diffusion coefficients and we include various forms of non-linear diffusion.

The initial conditions are

$$u_1(x,0)=g_0\left(x\right),u_2(x,0)=h_0\left(x\right),$$

(3)

and boundary conditions

$${\displaystyle \frac{\partial u_1}{\partial x}}(0,t)={\displaystyle \frac{\partial u_1}{\partial x}}(L,t)=0,{\displaystyle \frac{\partial u_2}{\partial x}}(0,t)={\displaystyle \frac{\partial u_2}{\partial x}}(L,t)=0,$$

(4)

where L is the length of domain x.

2. Steady State Solutions

Steady-state solutions are the intersection points of the two reaction parts in (1)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{u}_1}{dt}}& =& u_1(f\left(u_1\right)-b_1{u}_2),\hfill \\ \hfill {\displaystyle \frac{d{u}_2}{dt}}& =& u_2(g\left(u_2\right)-b_2{u}_1),\hfill \end{array}$$

(5)

where $f\left(u_1\right)=1+a_1{u}_1-(1+a_1)u_1^2$ and $g\left(u_2\right)=1+a_2{u}_2-(1+a_2)u_2^2$.

Our work focuses on studying the existence of a pattern for (1) at the positive coexistence steady state solutions $E(u_1^{*},u_2^{*})$ of (5). The existence and stability of $E(u_1^{*},u_2^{*})$ are studied in detail in [25] and have been categorized into seven topological regions. The maximum number of coexistence steady states are three. Furthermore, two types of bifurcation have been also shown for this model in [25]. Namely, Saddle-Node and Transcritical bifurcation.

3. Conditions of Turing Instability

In this section, we will examine the probability of the emergence of patterns for 1. We will follow the hypothesis, that the change of stability of the steady state solution $E(u_1^{*},u_2^{*})$ under the impact of diffusion term may lead to the existence of patterns. The linear stability analysis will be used in the following sections to find the critical wave number and critical value for the bifurcation parameter. The Turing conditions will be derived for the cases, without and with diffusion terms. This can be done using linearisation concepts.

4. Linearization of Reaction System

Firstly, we study the stability of $E(u_1^{*},u_2^{*})$ in (5) and we exclude the diffusion term. We follow the same technique as in (Cross & Greenside, 2009; Turing, 1990) [14,16]. We introduce small perturbations around the steady state so that $u_1={u_1}^{*}+\stackrel{⌢}{u_1}$ and $u_2=u_2^{*}+\stackrel{⌢}{u_2}$, where $\left|\stackrel{⌢}{u_1}\right|,\left|\stackrel{⌢}{u_2}\right|\ll 1$. Now we replace $u_1$ and $u_2$ in (5) with $exp(\lambda t$ and differentiate the left side, then we get

$$\dot{W}=KW,W=\left[\begin{array}{c}{u}_1-u_1^{*}\\ u_2-u_2^{*}\end{array}\right]$$

$$K=\left[\begin{array}{cc}{K}_{11}& K_{12}\\ K_{21}& K_{22}\end{array}\right],$$

(6)

The characteristic equation of 6 arises from

$$det(K-\lambda I)=0\Rightarrow {\lambda }^2-tr\left(K\right)\lambda +det\left(K\right)=0$$

(7)

where $tr\left(K\right)$ and $det\left(K\right)$ are refers to the trace and the determinant of K, respectively and mathematically are defined to be

$$tr\left(K\right)=K_{11}+K_{22},$$

(8)

$$det\left(K\right)=K_{11}{K}_{22}-K_{12}{K}_{21},$$

(9)

where

$$\begin{array}{cc}\hfill & K_{11}=1+2a_1{u_1}^{*}-3(1+a_1){{u_1}^{*}}^2-b_1{u_2}^{*},K_{12}=-b_1{u_1}^{*},\hfill \\ \hfill & K_{21}=-b_2{u_2}^{*},K_{22}=1+2a_2{u_2}^{*}-3(1+a_2){{u_2}^{*}}^2-b_2{u_1}^{*}.\hfill \end{array}$$

The steady state $E(u_1^{*},u_2^{*})>0$ is stable if

$$tr\left(K\right)<0,det\left(K\right)>0.$$

(10)

In the rest of the paper, we focus on steady-state solutions $E(u_1^{*},u_2^{*})$ which satisfy the two conditions in (10), i.e. we study stable $E(u_1^{*},u_2^{*})$ in the case of missing diffusion part. Next, we try to find if it is possible to break at least one of the conditions in (10) by adding the diffusion term. This will help us to investigate the region where may pattern solution exists in (1).

5. Linearization of Reaction-Diffusion System (RDS)

The linearization in this section will cover three cases depending on the form we give to $J_{1,2}$ in the diffusion term. The linearization of (1) and (2) is

$$\dot{W}=KW+D{\nabla }^{2}W,W=\left[\begin{array}{c}{u}_1-u_1^{*}\\ u_2-u_2^{*}\end{array}\right]$$

where

$$D=\left[\begin{array}{cc}{D}_{1,1}& D_{1,2}\\ D_{2,1}& D_{2,2}\end{array}\right]=\left[\begin{array}{cc}2{\alpha }_1{u}_1^{*}+{\beta }_1{u}_2^{*}& {\beta }_1{u}_1^{*}\\ {\alpha }_2{u}_2^{*}& 2{\alpha }_2{u}_2^{*}+{\beta }_2{u}_1^{*}\end{array}\right]$$

and K is the same as in 4. We assume solutions of the form $exp(\lambda t+ik.x)$ where $k\ge 0$, then the dispersion relation can be obtained in the form

$${\lambda }^2+[k^{2}tr\left(D\right)-tr\left(K\right)]\lambda +det(K-k^{4}D)+s\left(k^2\right)=0$$

(11)

where

$$\begin{array}{ccc}\hfill s\left(k^2\right)& =& [(-b_1{u_2}^{*})\left({\alpha }_2{u}_2^{*}\right)+(-b_2{u_1}^{*})\left({\beta }_1{u}_1^{*}\right)\hfill \\ \hfill & -& (1+2a_1{u_1}^{*}-3(1+a_1){{u_1}^{*}}^2-b_1{u_2}^{*})(2{\alpha }_2{u}_2^{*}+{\beta }_2{u}_1^{*})\hfill \\ \hfill & -& (2{\alpha }_1{u}_1^{*}+{\beta }_1{u}_2^{*})(1+2a_2{u_2}^{*}-3(1+a_2){{u_2}^{*}}^2-b_2{u_1}^{*})]k^2.\hfill \end{array}$$

(12)

Turing instability states that the values of $\lambda$ in (11) should be at least one positive. This will break the stability of $E(u_1^{*},u_2^{*})$ and support the hypothesis of the existence of patterns. This can be done whenever either $T{U}_1\left(k\right)<0$ or $T{U}_2\left(k\right)>0$, where:

$$\begin{array}{ccc}\hfill T{U}_1\left(k\right)& =& det(K-k^{4}D)+s\left(k^2\right)\hfill \\ \hfill T{U}_2\left(k\right)& =& k^{2}tr\left(D\right)-tr\left(K\right).\hfill \end{array}$$

(13)

In the following, we will focus on satisfying the first Turing condition $T{U}_1\left(k\right)$. Next, we try to define regions of the Turing instability and Turing bifurcation Threshold depending on the parameters in (2). For this purpose we fixed the parameters in reaction term 5, i.e. we put $a_1=0.5$, $b_1=0.7$, $a_2=0.4$ and $b_2=0.7$ to get the numerical steady state $E(u_1^{*},u_2^{*})=(0.76,0.74)$. A simple calculation using MATLAB program shows that the positive steady state $(0.76,0.74)$ satisfies the conditions in 10 and therefore is stable. In Figure 1a–c, three regions of Pattern formation are defined for the parameter values of ${\beta }_2$ as the x-axis and ${\alpha }_1$ as the y-axis. In the black regions, Turing instability conditions are satisfied, and unstable eigenvalues can be found.

5.1. Turing Bifurcation Threshold

The Turing bifurcation Threshold is the minimum values of the wave number $k^2$ which make $T{U}_1\left(k\right)=0$ for the parameter values in the nonlinear diffusion term (2). Figure 2 explains the region of Threshold, where $x-axis$ represented by $k^2$ and $y-axis$ by ${\beta }_2$ and, when ${\alpha }_1=-2$, ${\alpha }_2={\beta }_1=1$. The unstable region is above the Threshold curve while the stable region is under the threshold curve. Furthermore, we show positive eigenvalues of (11) whenever $T{U}_1\left(k\right)>0$. This leads to change of the stability of steady state solutions with the present of diffusion term as can be seen in Figure 3, the values are $-10\le {\beta }_2<0$ and ${\alpha }_1=-2$, ${\alpha }_2={\beta }_1=1$, $a_1=0.5$, $b_1=0.7$, $a_2=0.4$ and $b_2=0.7$. The steady-state $(0.76,0.74)$ loses stability and becomes unstable for some values of wave number $k^2$. Thus, we conclude that the stability of $E(u_1^{*},u_2^{*})$ can be broken under the effect of the non-linear diffusion term. Next, we show the type of pattern solutions for our model and general nonlinear diffusion, cross-diffusion, and self-diffusion forms.

6. Pattern Formations for RDS (1) and (2) in Two Dimensions

In this section, we have used one of the efficient numerical methods, the Finite element method to find high-quality forms of patterns in the defined regions. The Turing instability conditions have been applied in the previous sections to our model in (1) and (2) and support the idea of existing patterns in our model. The Finite element software package namely COMSOL Multiphysics will applied next to study some type of patterns depending on the initial conditions and the domain size. The domain we choose is square with $-L\le x\le L$ and $-L\le y\le L$, where $L=200$. The mesh size is of type extra fine resolution. The method of time stepping Backward differentiation formula(BDF) with adoptive mesh is used, and the convergence criteria; relative tolerance of order ${10}^{-6}$ is applied to all the variables. The boundary conditions are of type zero-flux boundary conditions, while the initial conditions are chosen to be near the steady state solution and with small perturbation

$$u_1(x,y,0)=\stackrel{⌢}{u_1}{e}^{-(x^2+y^2)}\mathrm{and}{u}_2(x,y,0)=\stackrel{⌢}{u_2}{e}^{-(x^2+y^2)}.$$

The boundary conditions are expressed as ${\displaystyle \frac{\partial u}{\partial x}}=0$, ${\displaystyle \frac{\partial u}{\partial y}}=0$, ${\displaystyle \frac{\partial w}{\partial x}}=0$, ${\displaystyle \frac{\partial w}{\partial y}}=0$. Simulations use COMSOL Multiphysics with time $(T=40000)$ to ensure pattern convergence. For the consistency with theoretical results in the previous section, we choose the same parameters and the numerical steady-state value of $E(u_1^{*},u_2^{*})$. However, we vary the parameters in (2) to show the patterns for the case of general nonlinear diffusion and both cross-diffusion and self-diffusion. In Figure 4a,b, we assume the case of self-diffusion by setting the parameters as follows: ${\alpha }_2={\beta }_1=0$, ${\alpha }_1=-6$ and ${\beta }_2=-8$. The pattern solutions for both $u_1$ and $u_2$ show different forms of nice decoration of patterns. In Figure 4c,d, the nonlinearity of diffusion is defined to be cross-diffusion by setting ${\alpha }_1={\beta }_2=0$, ${\alpha }_1=-2$ and ${\alpha }_2=-6$.

The similarity in the patterns can be seen for both $u_1$ and $u_2$. Another form of patterns have been produced as shown in Figure 5a–d. This time we choose ${\beta }_2$ to be negative in Figure 5a,b and positive in Figure 5c,d.

In the solid state, the diffusion rates of different species generally differ by several orders of magnitude. This is because they are typically governed by nonlinear Arrhenius escape rates, $exp(-E/k_{B}T)$, where the migration barrier E can range from fractions of an electron volt to several electron volts. In other words, the diffusion term in the solid state model can be written as a nonlinear term (similar to the case in our model)during the diffusion of different materials in the solid form till it ends up producing some shape of patterns. For example, the competition between Crystal materials that produce some types of Patterns forms [35,36]. Figure 5 and Figure 6 illustrate the pattern formation that emerges from unstable fixed points, or stationary states, which have been spatially perturbed using Gaussian functions. The numerical simulations display regions characterized by periodic structures, which are interspersed with irregular and complex boundaries. This pattern bears a resemblance to grain boundaries observed in multi-crystal materials. The intriguing complexity of these patterns stems from the unique characteristics of the specific model under study.

7. Conclusions

In this paper, we extend the study of a Reaction-Diffusion model in Ecology, with two interaction species. The reaction terms represent competition between two species and cooperation within species. The spatial uniform solutions for this reaction equation have been previously studied but only with constant diffusion terms. However, we expand the study of this model assuming that the diffusion coefficient is nonlinear. We have shown that this model satisfied the Conditions of Turing instability for specific steady-state solutions and for some values of wave number theoretically. Then we studied the pattern formation solutions for this model with nonlinear Diffusion coefficients using COMSOL MULTIPHYSICS software. We demonstrate diverse Patterns for different values of parameters in the Diffusion coefficient terms specially when we have the cases of cross- and self-diffusion. The regions where patterns can be formulated depending on two parameters in Diffusion coefficient parts have been classified. Finally, the Turing bifurcation Threshold has been found and plotted for a numerical steady-state solution. Unlike the constant diffusion model ([25]), our extended model has exhibits Turing instability and produce some types of Pattern solutions.