Modelling and analysis of an alcoholism model with treatment and effect of Twitter

1.   Introduction

The phenomenon of vegetation patterns has been found in a great many semiarid areas in recent decades ^[1,2,3,4]. Klausmeier ^[5] has proposed a semiarid vegetation model to explain the formation of striped vegetable patterns from the Turing instability and pointed out that the regular pattern would be self-organized from an irregular initial state. Klausmeier also showed that the diffusion induced by Brownian movement can act on the formation of regular patterns of plant vegetation. In addition, some more complicated spatiotemporal patterns were founded in the semiarid vegetation systems (see for example ^[6,7,8]). Some partial diferential equations mathematical modles were used to describe the vegetation patterns (see for example ^[9,10,11]). The stable steady solution were studied to demonstrate the pattern formation in ^[12,13]. Some frameworks of analysise were proposed to investigate the vegetation phenomena (see for example ^[14,15,16]). Recently, in ^[17,18], semiarid vegetation models have been characterized by fractional diffusion equations. In ^{[19,20,21,22]}, the modified Klausmeier models were investigated by hyperbolic equations.

However, most of the research involved a reaction-diffusion system with plant seeds spreading through a continuous space. Frequent human activities have divided the space into many fine-scale habitats that form a network connected with spreading seeds in reality. The network can be depicted by a weighted graph. The graph is denoted by ${\mathcal G} = (V, E)$, including the vertices $V = \{1, 2, \cdots, n\}$ and the edges $E$. If the vertex $y$ is adjacent to the vertex $x$, then we set $y\sim x$. If every adjacent $x$ and $y$ is given by a weight function $\omega$, then we get a weighed graph. Here $\omega:\; V\times V\rightarrow [0, \infty)$ is a positive function satisfying $\omega(x, y) = \omega(y, x)$ and $\omega(x, y) > 0$ only when $x\sim y$. Hence, the degree of vertex $x$ is $D_\omega (x) = \sum_{y\sim x, y\in V} \omega (x, y)$. The weighted graph Laplacian operator is defined as follows:

Recent research work of pattern formation for the networked dynamics in physics illustrates that Turing patterns can occur with the large network ^[23]. The further recent research work exhibits that Turing-like waves can take place in a networked one-species spatiotemporal dynamics with delay ^[24], which is different from the convectional theory that Turing instability can only take place in the two-component spatiotemporal dynamics ^[25,26]. We focus on unexplored Turing patterns that occur in the weighed networked semiarid vegetation models representing fragmented habitats.

We consider the vegetation dynamics with the following weighted network:

where $U$ and $V$ represent densities of water and plants, respectively. Water is supplied at the rate of $A$ uniformly and loses at the rate of $LU$ because of evaporation. Plants absorb water at the rate of $RU^2V$. $J$ is the amount of plants biomass per unit of water consumed. Plants lose at the mortal rate of $MV$. The diffusion of water and plants are constructed by the graph Laplacian diffusion coefficient $B$ and $D$, respectively. The original Klausmeier model ^[5] contains the sloped terrain, while our model, a modified Klausmeier model, only studies flat terrain. Various techniques were proposed to show the existence and uniqueness of solutions for graph Laplacian equations (see for example ^{[27,28,29,30]}). In ^[31,32], the Lyapunov functions were constructed to show the global stability of reaction-diffusion systems where the space was confined in a graph. In ^[33,34,35], the method of stablity analysis were developed to deal with the networked reaction-diffusion systems. In ^[36,37], the pattern formation was investigated to describe the complex dynamical behavior of graph Laplacian equations.

For the sake of simplicity, scaling the variables $u = R^{1/2}L^{-1/2}JU, \; v = R^{1/2}L^{-1/2}V, \; \tilde t = Lt$ and then dropping the tildes, the dimensionless, weighted networked spatiotemporal dynamics can be taken as:

where $a = AJR^{1/2}L^{-3/2}$, $b = M/L$, $c = BJR^{1/2}L^{-3/2}$ and $d = DJR^{1/2}L^{-3/2}$.

This paper is organized as follows. In Section 2, we show the global existence of solutions to the system (1.1). In Section 3, we obtain the Turing parameter space from the analytic results of the linear stability of the positive equilibrium and ensure that the Turing bifurcation takes place before Hopf bifurcation. In Section 4, an amplitude equations near the Turing instability critical point is derived by a weakly nonlinear analysis. Further, the stability of Turing patterns are considered by analyzing the amplitude equations. In Section 5, we establish the framework for deriving optimal control strategies of rainfall by deriving the adjoint equations and employing optimal control theory. In Section 6, some numerical simulations are presented to verify the theoretical analysis and explored the effects of the graph Laplacian diffusion on Turing patterns. Finally, some discussions and conclusions are given in Section 7.

2.   Existence and uniqueness of solutions

In this section, we derive the existence and uniqueness result for the solution of system (1.1) by constructing the monotone iterative sequence. We need to give the definition of coupled upper and lower solutions, which was initially proposed by ^[38].

Definition 2.1. If $\tilde u ,\; \hat u, \; \tilde v ,\; \hat v \in C[0, T]$ are differentiable in $(0, T]$, then pair of functions $\mathbf{\tilde u} = (\tilde u, \tilde v)$ and $\mathbf{\hat u} = (\hat u, \hat v)$ are called coupled upper and lower solutions of (1.1) if $\mathbf{\tilde u}\geq\mathbf{\hat u}\geq 0$ and if

Since the reaction terms of system (1.1) are Lipschitz continuous, we denote $K$ as the Lipschitz constant. By setting $\mathbf{\overline u}^{(0)} = \mathbf{\tilde u}$ and $\mathbf{\underline u}^{(0)} = \mathbf{\hat u}$ as the initial iterations, we construct the sequences $\mathbf{\overline u}^{(m)}$ and $\mathbf{\underline u}^{(m)}$ from the following iteration system:

Here $\mathbf{\overline u}^{(m)} = ({\overline u}^{(m)}, {\overline v}^{(m)})$ and $\mathbf{\underline u}^{(m)} = ({\underline u}^{(m)}, {\underline v}^{(m)})$. Since system (2.2) is a scalar, ordinary differential equation, the sequences $\mathbf{\overline u}^{(m)}$ and $\mathbf{\underline u}^{(m)}$ exist and are unique for some $T$. We have obtained the following monotonic properties of these sequences.

Lemma 2.1. If $\mathbf{\overline u}^{(m)}$ and $\mathbf{\underline u}^{(m)}$ are defined by (2.2), then the following holds:

Moreover, $\mathbf{\overline u^{(m)}}$ and $\mathbf{\underline u^{(m)}}$ are coupled upper and lower solutions of (1.1) for each $m = 1, 2, \cdots$.

Proof. Considering system (2.2) is a monotone system, we can extend the method of monotone semiflows in ^[38] to our graph Laplacian equations.

Let $\underline p^{(1)} = \underline u^{(1)}-\underline u^{(0)}$. Then, by (2.2), for $(x, t)\in V\times(0, T]$, $\underline p^{(1)}$ satisfies

where the last inequality comes from (2.1). Meanwhile, $\underline p^{(1)}(x, 0) = 0$ for $x\in V$. Using the maximum principle (Lemma 2.1 in ^[39]), we have $\underline p^{(1)}(x, t)\geq0$ for $(x, t)\in V\times[0, T]$. Therefore, $\underline u^{(0)}(x, t)\leq \underline u^{(1)}(x, t)$ for $(x, t)\in V\times[0, T]$.

Let $\underline q^{(1)} = \underline v^{(1)}-\underline v^{(0)}$. Then, by (2.2), for $(x, t)\in V\times(0, T]$, $\underline q^{(1)}$ satisfies

where the last inequality comes from (2.1). Meanwhile $\underline q^{(1)}(x, 0) = 0$ for $x\in V$. Using the maximum principle (Lemma 2.1 in ^[39]), we have $\underline q^{(1)}(x, t)\geq0$ for $(x, t)\in V\times[0, T]$. Therefore, $\underline v^{(0)}(x, t)\leq \underline v^{(1)}(x, t)$ for $(x, t)\in V\times[0, T]$.

Thus, we have

By a similar argument as above, we also have

Now, we set $\phi^{(1)} = \overline u^{(1)}-\underline u^{(1)}$. By (2.2), it follows that for $(x, t)\in V\times(0, T]$,

where the last inequality comes from $0\leq\underline v^{(0)}(x, t)\leq \overline v^{(0)}(x, t)$. Moreover, $K$ is the Lipschitz constant of (1.1). Hence, we have

The initial condition is $\phi^{(1)}(x, 0) = 0$ for $x\in V$. Applying the maximum principle (Lemma 2.1 in ^[39]) yields that ${\overline u}^{(1)}\geq{\underline u}^{(1)}$. By a similar argument for $v$, we have ${\overline v}^{(1)}\geq{\underline v}^{(1)}$. Therefore,

Combining (2.4)–(2.6) yields that

Then, we show that $\mathbf{\overline u}^{(1)}$ and $\mathbf{\underline u}^{(1)}$ are coupled upper and lower solutions of (1.1). It remains to show that $\mathbf{\overline u}^{(1)}$ and $\mathbf{\underline u}^{(1)}$ satisfy (2.1). By (2.2) and (2.7), for $(x, t)\in V\times(0, T]$, we have

where the last inequality holds because $K$ is defined in (1.1). In a similar argument, we obtain that $\mathbf{\overline u^{(1)}}$ and $\mathbf{\underline u^{(1)}}$ are coupled upper and lower solutions of (1.1).

By an induction method, we extend the above result to some arbitrary $m$. We choose $\mathbf{\tilde u} = \mathbf{\overline u}^{(1)}$ and $\mathbf{\hat u} = \mathbf{\underline u}^{(1)}$. After the above process, we have

Thus, we can obtian (2.3). □

Theorem 2.1. If $\mathbf{\tilde u}$ and $\mathbf{\hat u}$ are a pair of coupled upper and lower solutions of system (1.1), then there exists $\mathbf{u^*}$ in $< \mathbf{\hat u}, \mathbf{\tilde u} >$ the unique solution of system (1.1) for $t\in [0, T]$.

Proof. In terms of Lemma 2.1, the sequences $\mathbf{\overline u}^{(m)}$ and $\mathbf{\underline u}^{(m)}$ are convergent for $t\in [0, T]$. We set

We directly solve the solution of system (2.2) as follows:

Since $\mathbf{\overline u^{(m)}}$ and $\mathbf{\underline u^{(m)}}$ are bounded in $< \mathbf{\hat u}, \mathbf{\tilde u} >$, by the dominated convergence theorem, it follows that $\overline u(x, t)$ and $\underline u(x, t)$ satisfy

Since $K$ is the Lipschitz constant, we obtain

By (2.2), it follows that

Since $\overline v^{(m)}$ are bounded in $(\hat v, \tilde v)$ for $(x, t)\in V\times [0, T]$, applying the dominated convergence theorem yields that

In term of the definition of $K$, we have

We now estimate the boundedness ofthe above inequality. We have

where $n$ is the number of the vertices of graph. Thus, we have

where $C_1$ and $C_2$ only depends on $K$, $c$, $n$, and $\max_{x\in V} D_\omega (x)$.

Thus, there exists a constant $T_0: = 1/2(C_1+C_2)$ such that $||\overline u-\underline u||_\infty = 0$ and $||\overline v-\underline v||_\infty = 0$, that is $\overline u\equiv \underline u$ and $\overline v\equiv \underline v$ hold for $t\in (0, T_0]$. Since $C_i$ does not depend on the initial value, we can extend $T_0$ to the time $T$. By setting $\overline u = u^*$ and $\overline v = v^*$, we obtain

Hence, $(u^*, v^*)$ is the unique solution of system (1.1) for $t\in [0, T]$. □

Theorem 2.2. The solution of system (1.1) exists and is unique for all $t\in [0, \infty)$.

Proof. In view of Theorem 2.1, system (1.1) has an unique solution $\mathbf{u^*}$ for $t\in [0, T]$. In order to show that $T = \infty$, we need to show a priori upped bound of the solution by the approach of upper and lower solutions. We set $(\hat u, \hat v) = (0, 0)$ and $(\tilde u, \tilde v)$. Thus,

Since they satisfy

$(\tilde u, \tilde v)$ and $(\hat u, \hat v)$ are a pair of upper and lower solutions of system (1.1). By Theorem 2.1, the solution of system (1.1) is the uniformly upper bounded. Hence, $T = \infty$. □

3.   Stability analysis

This section is devoted to obtaining the existence of the Turing bifurcation around the uniform equilibrium of system (1.1). System (1.1) has three equilibria: $E_0$, $E_1$ and $E_2$. Here $E_0 = (a, 0)$, $E_1 = (\frac{a+\sqrt{a^2-4b^2}}{2}, \frac{2b}{a+\sqrt{a^2-4b^2}})$, and $E_2 = (\frac{a-\sqrt{a^2-4b^2}}{2}, \frac{2b}{a-\sqrt{a^2-4b^2}})$. Considering the biological significance of the system (1.1), we only focus on the case of the positive equilibrium existing. Thus, the case $a > 2b$ always holds. By performing linear stability analysis, the second positive equilibrium $E_1$ is an unstable saddle. Therefore, we only consider the case of the third positive equilibrium $E_2\equiv(u_s, v_s)$ existing. In what follows, we give the space decomposition of $[L^2(V)]^2$ with respect to graph Laplace.

Lemma 3.1. The eigenvalue problem

for

admits eigenvalues

where the associated eigenfunctions are $\{\phi_i\}_{i = 1}^n$. Moreover, $[L^2(V)]^2$ space decomposition

holds. Here, $\mathbf{E_{i}}: = \{\mathbf{c}\cdot \phi_{i}:\ \mathbf{c}\in \mathbb{R}^2\}$.

Proof. Let $\lambda$ and $\phi(x)$ be a solution of (3.1). By multiplying $\phi(x)$ and integrating it over $V$, we have

We define a bilinear functional

the domain $D(F)$ of $F$ is $D(F): = \{u:\; u\in L^2(V), ||u||_{L^2(V)} = 1\}.$

Step 1: Determine $\lambda_1$. In light of the definition of $F$ in (3.3) and (1.1), $F$ is continuous in $D(F)$ and bounded from below, that is, there exists $\phi_1\in D(F)$ such that

Here, $\lambda_1$ and $\phi_1$ are a solution of the eigenvalue problem (3.1). By substituting $\lambda_1 = 0$ and $\phi_1(x) = 1$ into (3.1), it easy to see that $\lambda_1 = 0$ is the first eigenvalue. Moreover, $\lambda_1$ is simple.

Step 2: Determine $\lambda_2$. In view of the continuous property of $F$, the second minimization problem admits a solution, that is, $\phi_2(x)\in L^2(V)$ with $||\phi_2||_{L^2(V)} = 1$ and $\phi_2\perp\phi_1$ such that

where $u\perp\phi_1$ means that $u$ and $\phi_1$ are orthogonal, that is, $\int_V u\phi_1dx = 0$. By the definition of (3.5), we have $\lambda_2 > \lambda_1$.

By performing the above steps $n$ times, we can construct a sequence $\{\lambda_k\}_{k = 1}^n$ whose associated eigenfunction is $\{\phi_k\}_{k = 1}^n$. Owing to the finite nature of $V$, we find that $n$ is finite.

Step 3: $[L^2(V)]^2$ space decomposition. We denote $\mathbf{E}_i$ as $\mathbf{E_{i}}: = \{\mathbf{c}\cdot \phi_{i}:\ \mathbf{c}\in \mathbb{R}^2\}$. Owing to the orthogonality of $F(\phi_i)$ for each $i = 1, 2, \cdots, n$, we obtain $[L^2(V)]^2 = \bigoplus_{i = 1}^n \mathbf{E}_i$. □

Theorem 3.1. If

holds, and $(u_s, v_s)$ is locally asymptotically stable.

Proof. The linearization of system (1.1) at the equilibrium $(u_s, v_s)$ is

with

In view of Lemma 3.1, the eigenvalue of $({\mathbf J}+{\mathbf D}\Delta_\omega)$ is equilivalent to the matrix $({\mathbf J}-\lambda_i{\mathbf D})$. Then, we can write the solution of (3.7) as follows:

Substituting (3.8) into (3.7), we obtain

Since $\phi_i$ and $e^{\sigma t}$ are nonzero vectors, we can obtain

Hence, $\sigma$ is sastified

This leads to the characteristic polynomials of (1.1)

where $g_i = 1+b+v_s^2-2u_sv_s-c\lambda_i-d\lambda_i$ and $h_i = (1+v_s^2-c\lambda_i)(b-2u_sv_s-d\lambda_i)+2u_sv_s^3$.

In view of (3.6), we obtain that $g_i > 0$ and $h_i > 0$. Hence, $\sigma < 0$. We conclude that $(u_s, v_s)$ is locally asymptotically stable. □

Theorem 3.2. System (1.1) presents a Turing bifurcation at $b = 2u_sv_s+d\lambda_i-\frac{2u_sv_s^3}{1+v_s^2-c\lambda_i}$.

Proof. According to Theorem 3.1, the characteristic polynomial of (1.1) is

where $g_i = 1+b+v_s^2-2u_sv_s-c\lambda_i-d\lambda_i$ and $h_i = (1+v_s^2-c\lambda_i)(b-2u_sv_s-d\lambda_i)+2u_sv_s^3$.

Turing bifurcation point refers to $g_i > 0$ and $h_i = 0$. Thus, the Turing bifurcation critical value $b_c$ is

Then, system (1.1) has a Turing bifurcation at $b = 2u_sv_s+d\lambda_i-\frac{2u_sv_s^3}{1+v_s^2-c\lambda_i}$. □

Remark 3.1. Although there were many standard techniques for Turing bifurcation of reaction-diffusion equations in the past, our method is aimed at network diffusion equations, which is a major advancement in bifurcation dynamics in network equations.

4.   Weakly nonlinear analysis

In this section, we investigate the bifurcation behavior in the area of the Turing bifurcation point $b_c$ by using a weakly nonlinear analysis of system (1.1).

Theorem 4.1. If (3.6) and $\Gamma+4v_s\rho_2+2u_s\rho_2^2 > 0$ hold, the Turing bifurcation of system (1.1) is stable.

Proof. By transforming $(u_s, v_s)$ into $(0, 0)$, we first rewrite the system (1.1) into the following form:

where

We introduce the parameter $\varepsilon$ satisfying $\varepsilon^2 = \frac{b-b_c}{b_c}$ by the method of perturbation technique. Let $T = \varepsilon^2 t$ be the slow time of the Turing pattern bifurcation and expand $\mathbf{U}$ with respect to the parameter $\varepsilon$ as follows:

Therefore, $\mathbf{J}$ is

where

Expanding (4.1) with respect to different orders of $\varepsilon$, we can get three equations at orders $\varepsilon^j(j = 1, 2, 3)$:

where

Since $(u_1, v_1)^T$ is the linear combination of the eigenvector of system (3.9) associated $\sigma = 0$. Thus, the solution at $O(\varepsilon)$ is $(u_1, v_1)^T = {\mathbf \rho}A(T) \phi_i$, and ${\mathbf \rho} = (\rho_1, \rho_2)^T = (1, -\frac{1+v_s^2-c\lambda_i}{2u_s v_s})^T$. Here, $\phi_i$ is the eigenvector corresponding to $\lambda_i$ and $\phi_i$ is the spatial wave function of Turing pattern. $A(T)$ is the amplitude of the solution determined by the higher order perturbation term of $\varepsilon$ and is still unknown.

Next, we have the equation associated $O(\varepsilon^2)$

As a result, we can obtain that

where $\rho_3 = \frac{b_c-\lambda_id}{\Theta}$ and $\rho_4 = \frac{\lambda_ic-1}{\Theta}$, $\Theta = (1+v_s^2-\lambda_ic)(b_c-2u_sv_s-\lambda_id)+2u_sv_s^3$.

Further, we consider $O(\varepsilon^3)$. The equation associated $O(\varepsilon^3)$ is written as follows:

where $\Gamma = \rho_2^2+2(2v_s\rho_2+u_s\rho_2^2)(u_s\rho_2\rho_4+v_s\rho_2\rho_3+v_s\rho_4)$.

According to Fredholm solubility condition, $({\mathbf J_c}+{\mathbf D}\Delta_\omega)$ has the adjoint operator ${\mathbf L_c^+}$, and the operator ${\mathbf L_c^+}$ has the nontrivial kernel $(1, (v_s^2+1-c\lambda_i)/v_s^2)^T\phi_i$. Multiplying (4.5) by $(1, (v_s^2+1-c\lambda_i)/v_s^2)^T\phi_i$ and integrateing it on $V$, we can induce Turing bifurcation's amplitude equation

where $\sigma = -\frac{b_cA(v_s^2+1-c\lambda_i)}{v_s^2}$, and $L = \frac{1-c\lambda_i}{v_s^2}(\Gamma+4v_s\rho_2+2u_s\rho_2^2).$

Therefore, the Turing bifurcation is stable if

holds, which is called supercritical bifurcation. □

5.   The rainfall effect of system (1.1)

This section discusses the effect of rainfall on the networked system (1.1). We will study how rainfall determine the dynamics of the networked structure. Hence, we involve the decreasing rainfall due to climate change in system (1.1). We use the term $W(t)$ to represent the probablity of decrease in rainfall. When including the probablity of decreasing in rainfall $W(t)$, the vegetation system takes the following form:

thus, the objective functional is written as

over the control set:

Here, we define the positive constant $A$ to assign weight to the probablity of rainfall in the objective function. The optimal rainfall control problem can be expressed as follows:

Find $W^*\in V$ such that

subject to system (5.1).

For the above optimal control problem, after a similar argument, it is easy to show the solution of system (5.1) exists. To apply Pontryagin's maximum principle ^[40], we need to establish the existence of an optimal control.

Theorem 5.1. There exists an optimal control $W^*$, and corresponding solution $(u^*, v^*)$ of system (5.1) that minimizes $J(W)$ defined by (5.2) over $V$.

Proof. We refer to the conditions described in Theorem III.4.1 by Fleming and Rishel ^[41]. The requirements on the set of admissible controls $V$ and on the set of end conditions are clearly satisfied here. Moreover, for any function $W(t)$ in $V$ given in (5.3), using a similar argument in Theorem 2.2, there exists a unique solution $(u, v)$ of system (5.1). It remains to show that $(u, v)$ and $W(t)$ are globally existing.

After an argument similar to the proof of 2.1, we have

By using the dominated convergence theorem, it follows that for $t\in [0, T]$ the limits $\overline u(x, t)$ and $\underline u(x, t)$ satisfy

Since $K$ satisfies the Lipschitz condition, we have

In a similar way, we have

In terms of the definition of weighed Laplacian operator, we have

Since $0\leq W(t)\leq 1$, we have

where $C_1$ and $C_2$ only depends on $K$, $c$, $n$, and $\max_{x\in V} D_\omega (x)$.

Thus, there exists a constant $T_0: = 1/2(C_1+C_2)$ such that $||\overline u-\underline u||_\infty = 0$ and $||\overline v-\underline v||_\infty = 0$, that is, $\overline u\equiv \underline u$ and $\overline v\equiv \underline v$ hold for $t\in (0, T_0]$. Since $C_i$ is independent on the initial value, we can extend $T_0$ to $T$. Thus, $(u, v)$ and $W(t)$ are globally existing. □

Theorem 5.2. If the optimal control $W^*$ and corresponding solution $(u^*, v^*)$ of system (5.1) minimize the objective functional (5.2), then there eixsts adjoint variables $\lambda_1$ and $\lambda_2$ satisfying

Furthermore, this optimal control is characterized by

Proof. We shall use Pontryagin's minimum principle ^[40] to complete the proof. We define the following Hamiltonian:

It is easy to check that $\frac{\partial^2 H}{\partial W^2} = 2 > 0$. Thus, the critical point $W^*$ is indeed a minimum. By applying Pontryagin's minimum principle ^[40], the adjoint equations are given by

and must satisfy transversality condition $\lambda_i(T) = 0$ for $i = 1, 2$. Hence, (5.5) holds.

Next, since $W^*$ is a critical point of the Hamiltonian, it follows that $\frac{\partial H}{\partial W} = 0$ at $W^*$. This yields the following condition on the optimal control:

Since $W^*$ must belong to $W$, we have $W^* = \max\left(\min\left(1, \frac{\lambda_1u}{2}\right), 0\right)$. □

Remark 5.1. The biological interpretation of the optimal rainfall $W^*$ in (5.6) is read as: When $W(t)$ satisfies $0 < W(t) < 1$, then the optimal rainfall $W^* = \frac{\lambda_1u}{2}$. The optimal rainfall probability is determined only by the density of water, not by the density of plants.

6.   Numerical simulations

In this section, we shall carry out some numerical simulations to validate and extend analytical results of supercritical Turing patterns. In view of the conditions (3.6) and (4.7) determining the stability of the Turing bifurcation, we take the following set of parameters coming from the semiarid vegetation model ^[5]:

We show the real part of the eigenvalue characteristic polynomials (3.9) corresponding to the wave number $\lambda_i$ in Figure 1. When the solid line is above the dashed line, Turing bifurcation occurs. In addition, we describe that the behavior of system (1.1) is determined by the input rate of water $a$ and the death rate of plants $b$. By using $b_c$ as the critical value of the Turing bifurcation, when the plant Laplacian diffusion coefficient $d = 242.5$, the Turing parameter space is shown in Figure 2. Considering the decrease of water input or the increase of plants loss, a transition from homogeneous vegetation to no vegetation is predicted in Figure 2.

Furthermore, we construct a Watts-Strogatz network consisting of a ring grid with 100 vertices connected to the nearest 4 neighbors from the left and right directions. In view of the fact that the motion of water usually obeys random walk, we add edges to the graph to satisfy the edge with a rewired probability of 0.15 for each pair of vertices. Through the use of numerical approximation, we obtain the solution of $u$ and $v$. Figure 3 displays the plant density with network at different time. We describe the initial plant density with a small perturbation near the equilibrium, which can not effect the formation of spatiotemporal patterns as shown on the first panel of Figure 3. By means of the time integration of Runge-Kutta scheme, we exhibit the evolution of the spatial density of the plant with network at $t = 50, 75$, and $100$ as shown in the $2$nd, $3$rd and $4$th panels of Figure 3, respectively. From Figure 3, we can see that as the time $t$ increases, the density of plant exhibits different value for each fixed vertex, that is, the system shows that a Turing pattern occurs under Turing parameter space $a = 2$ and $b = 0.45$. Though spatial patterns for different time are illustrated in these figures, we can not get the asymptotic behavior of the plant density for each vertex. Therefore, we describe that the asymptotic behavior converges to the temporal solution of system (1.1) for each vertex, where each node is arranged in a vertical column without considering the network structure. In Figure 4, the Turing pattern is stable and the solution converges to the temporally homogeneous behavior for each node in the left panel, the average density of plant converges to a constant with time in the middle panel, and the average density of plant over time is inhomogenous with respect to each node in the right panel. Hencen Figure 4 illustrates that the Turing pattern has a stable dynamical behavior under a large diffusion rate $d = 242.5$.

In Figure 5, we plot the 2-dimensional continuous spatial Turing pattern for the plant density, the solution of which has generated a spatial regular pattern based on a self-organizational continuous diffusion mechanism. In its right panel, the average density of plants over continuous space converges to a constant with time. Comparing to the different spatial type, we find that the speed of the Turing self-organizational process is not different. The 2-dimensional continuous spatial system (middle panel of Figure 4) is faster than the networked system (right panel of Figure 5).

7.   Discussions

The reaction-diffusion system of a semiarid vegetation model is studied in this paper with a network considering the fragmented habitats due to the dispersion of plant seeds. The existence and uniqueness of solutions are shown by the monotone iterative technique. The stability of the Turing bifurcation is studied by the technique of weakly nonlinear analysis, which brought about the occurrence of the stable Turing patterns. We also propose a nonstandard finite difference method for the numerical solution. Furthermore, we demonstrate that the graph Lapacian diffusion plays an important role in the pattern formation.

Our results indicate that both the plant loss and rainfall have an impact on the Turing parameter space of a semiarid vegetation model with weighted networks. By Theorem 2.2, we have shown that the system with weighted networks admits an unique global solution. Theorems 3.1 and 3.2 indicate that if the plant loss $b$ crosses to $b_c$, the positive equilibrium of the system has a different dynamical behavior, which generates Turing instability. Theorem 4.1 indicates that when the Turing instability occurs, the Turing pattern is asymptotic stable.

It is shown that optimal rainfall probability is determined only by the density of water, not by the density of plants. Moreover, we have performed the numerical simulations to find the difference between these two self-organizational Turing patterns. The reaction-diffusion system has a self-organizational speed faster than that of the graph Laplacian diffusion model. In summary, this paper has proposed a new way to depict the movement of surface water and biomass with the graph Laplacian diffusion term defined in the semiarid vegetable model.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Research was partially supported by NSF of China (No. 12471476) and the NSF of Jiangsu Province (No. BK20231357).

Conflict of interest

The authors declare there is no conflict of interest.