Changes in soft coral <em>Sarcophyton</em> sp. abundance and cytotoxicity at volcanic CO<sub>2</sub> seeps in Indonesia

Hedi Indra Januar; Neviaty Putri Zamani; Dedi Soedarma; Ekowati Chasanah; Hedi Indra Januar; Neviaty Putri Zamani; Dedi Soedarma; Ekowati Chasanah

doi:10.3934/environsci.2016.2.239

AIMS Environmental Science

2016, Volume 3, Issue 2: 239-248. doi: 10.3934/environsci.2016.2.239

Previous Article Next Article

Research article

Changes in soft coral Sarcophyton sp. abundance and cytotoxicity at volcanic CO₂ seeps in Indonesia

1.
Indonesian Research and Development Center for Marine and Fisheries Products Processing and Biotechnology, KS Tubun Petamburan VI Street, Slipi, Central Jakarta 10260, Indonesia
2.
Department of Marine Science, Faculty of Fisheries and Marine Science, Bogor Agricultural University, Kampus IPB, Darmaga Raya Street, Bogor 16680, Indonesia

Received: 18 January 2016 Accepted: 19 April 2016 Published: 25 April 2016

This study presents the relationship between benthic cover of Sarcophyton sp. living on coral reefs and their cytotoxicity (an assumption of soft coral allelochemical levels) along acidification gradients caused by shallow water volcanic vent systems. Stations with moderate acidification (pH 7.87 ± 0.04), low acidification (pH 8.01 ± 0.04), and reference conditions (pH 8.2 ± 0.02) were selected near an Indonesian CO₂ seep (Minahasa, Gunung Api Island, and Mahengetang Island). Cover of the dominant soft coral species (Sarcophyton sp.) was assessed and tissue samples were collected at each site. The cytotoxicity tissue extracts were analyzed using the 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolinon bromide (MTT) method. Levels of cytotoxicity were strongly correlated with Sarcophyton sp. cover (p < 0.05; R² = 0.60 at 30 ppm and 0.56 at 100 ppm), being highest at mean pH 8.01 where the soft corals were most abundant. This finding suggests that Sarcophyton sp. can be expected to survive ocean acidification near Indonesia in the coming decades. How the species might be adversely affected by further ocean acidification later in the century unless CO₂ emissions are reduced remains a concern.

Keywords:

Citation: Hedi Indra Januar, Neviaty Putri Zamani, Dedi Soedarma, Ekowati Chasanah. Changes in soft coral Sarcophyton sp. abundance and cytotoxicity at volcanic CO2 seeps in Indonesia[J]. AIMS Environmental Science, 2016, 3(2): 239-248. doi: 10.3934/environsci.2016.2.239

Related Papers:

[1]	Eid Elghaly Hassan, Diping Zhang . The usage of logistic regression and artificial neural networks for evaluation and predicting property-liability insurers' solvency in Egypt. Data Science in Finance and Economics, 2021, 1(3): 215-234. doi: 10.3934/DSFE.2021012
[2]	Jawad Saleemi . Political-obsessed environment and investor sentiments: pricing liquidity through the microblogging behavioral perspective. Data Science in Finance and Economics, 2023, 3(2): 196-207. doi: 10.3934/DSFE.2023012
[3]	Paarth Thadani . Financial forecasting using stochastic models: reference from multi-commodity exchange of India. Data Science in Finance and Economics, 2021, 1(3): 196-214. doi: 10.3934/DSFE.2021011
[4]	Aihua Li, Qinyan Wei, Yong Shi, Zhidong Liu . Research on stock price prediction from a data fusion perspective. Data Science in Finance and Economics, 2023, 3(3): 230-250. doi: 10.3934/DSFE.2023014
[5]	Obinna D. Adubisi, Ahmed Abdulkadir, Chidi. E. Adubisi . A new hybrid form of the skew-t distribution: estimation methods comparison via Monte Carlo simulation and GARCH model application. Data Science in Finance and Economics, 2022, 2(2): 54-79. doi: 10.3934/DSFE.2022003
[6]	Yaya Su, Yi Qu, Yuxuan Kang . Online public opinion and asset prices: a literature review. Data Science in Finance and Economics, 2021, 1(1): 60-76. doi: 10.3934/DSFE.2021004
[7]	Xiaoling Chen, Xingfa Zhang, Yuan Li, Qiang Xiong . Daily LGARCH model estimation using high frequency data. Data Science in Finance and Economics, 2021, 1(2): 165-179. doi: 10.3934/DSFE.2021009
[8]	Bilal Ahmed Memon, Hongxing Yao . Correlation structure networks of stock market during terrorism: evidence from Pakistan. Data Science in Finance and Economics, 2021, 1(2): 117-140. doi: 10.3934/DSFE.2021007
[9]	John Leventides, Evangelos Melas, Costas Poulios . Extended dynamic mode decomposition for cyclic macroeconomic data. Data Science in Finance and Economics, 2022, 2(2): 117-146. doi: 10.3934/DSFE.2022006
[10]	John Leventides, Evangelos Melas, Costas Poulios, Paraskevi Boufounou . Analysis of chaotic economic models through Koopman operators, EDMD, Takens' theorem and Machine Learning. Data Science in Finance and Economics, 2022, 2(4): 416-436. doi: 10.3934/DSFE.2022021

Abstract

1. Introduction

Due to the difficulty in measuring the biomass of plankton, mathematical models of plankton populations are important methods for understanding the physical and biological processes of plankton. Various phytoplankton-zooplankton models (PZMs) have been studied: An NPZD model ^[1], a stochastic PZM ^[2], a plankton-nutrient system ^[3], a marine phytoplankton-zooplankton ecological model (PZEM) ^[4], Trophic phytoplankton-zooplankton models (PZM) ^[5], a bioeconomic PZM with time delays ^[6], a two-zooplankton one-phytoplankton PZM ^[7], a phytoplankton and two zooplankton PZM with toxin-producing delay ^[8], planktonic animal and plant PZM ^[9], and other PZMs ^{[10, 11, 12]}. In ^[12], Wang et al. investigated the following PZEM:

$\begin{eqnarray} \left\{ \begin{array}{l} \frac{dP}{dt} = \overset{plant\ growth}{\overbrace{rP\left( 1-\frac{P}{K} \right) }}-\overset{plant\ loss}{\overbrace{\frac{e_1PZ}{1+e_1h_1P+e_2h_2A}}}, \\ \frac{dZ}{dt} = \underset{animal\ gain}{\underbrace{\frac{n_1e_1PZ+n_2e_2AZ}{1+e_1h_1P+e_2h_2A}}}-\underset{natural\ death}{\underbrace{\tau Z}}-\underset{animal\ death}{\underbrace{\sigma \frac{PZ}{h+P}}}. \end{array} \right. \end{eqnarray}$ $

(1.1)

Fractional-order reaction-diffusion equation can describe the dynamic behavior of real systems more accurately. In recent years, fractional-order models have made remarkable progress in fields such as biomedicine, ecology, and public health. In biomedicine, fractional-order models have been used to study tumor growth ^[13] and hepatitis B virus (HBV) treatment ^[14], revealing their advantages in describing the complex dynamic behaviors of biological systems. In ecology, fractional-order models have demonstrated their potential in the study of ecosystem complexity by analyzing the dynamic behaviors of plant-herbivore interactions ^[15]. In public health, the fractional model has been applied to analyze the transmission of smoking behavior ^[16] and the correlation between human papilloma virus (HPV) and cervical cancer ^[17], highlighting their utility in public health research. These studies not only enrich the theoretical application of fractional models but also have made significant progress in numerical analysis, providing new tools and methods for solving complex biomedical and ecological problems ^{[13, 17]}. For example, in ecology, fractional PZEMs can more accurately describe the complex dynamics of ecosystems, including population fluctuations and interactions. Li et al. ^[18] studied an established fractional-order delayed zooplankton-phytoplankton model. Javidi and Ahmad ^[19] studied dynamic analysis of time fractional-order phytoplankton-toxic PZM. Kumar et al. ^[20] studied a fractional plankton-oxygen modeling. The proposed fractional-order PZEM extends traditional integer-order models by incorporating memory effects through fractional derivatives, which better capture long-term ecological interactions. This model accounts for toxic substances and additional food transmission, factors critical to understanding algal blooms, and population collapse. Models often neglect cross-diffusion and fractional dynamics, limiting their ability to simulate complex spatiotemporal patterns. Our work bridges this gap, offering insights into chaotic attractors and pattern formation under realistic ecological constraints. In this paper, we investigate the following fractional-order PZEM:

$\begin{eqnarray} \left\{ \begin{array}{l} \frac{\partial^{\alpha_1}P}{\partial t^{\alpha_1}} = \overset{plant\ growth}{\overbrace{rP\left( 1-\frac{P}{K} \right) }}-\overset{plant\ loss}{\overbrace{\frac{e_1PZ}{1+e_1h_1P+e_2h_2A}}}+ \overset{plant\ diffusion}{\overbrace{d_{11} \nabla^{2} P}}, \\ \frac{\partial^{\alpha_2}Z}{\partial t^{\alpha_2}} = \underset{animal\ gain}{\underbrace{\frac{n_1e_1PZ+n_2e_2AZ}{1+e_1h_1P+e_2h_2A}}}-\underset{natural\ death}{\underbrace{\tau Z}}-\underset{animal\ death}{\underbrace{\sigma \frac{PZ}{h+P}}}+\underset{plant\ and \ animal \ diffusion}{\underbrace{ d_{21} \nabla^{2} P+d_{31} \nabla^{2} Z}}, \end{array} \right. \end{eqnarray}$ $

(1.2)

where $ P $ and $ Z $ represent the population densities of phytoplankton and zooplankton, respectively, at time t. $ \nabla ^2 $ stands for the Laplacian operator of the two-dimensional plane, $ \varOmega $ is a bounded connected region with smooth boundaries $ \partial\varOmega $, and $ d_{11} $ and $ d_{31} $ represent the diffusion coefficients of phytoplankton and zooplankton, respectively. These parameters characterize the stochastic movement of individuals within a population from areas of high concentration to those of lower density. $ d_{21} $ is a cross-diffusion coefficient that indicates the effect of the presence of phytoplankton on the movement of zooplankton populations. $ \alpha_i (i = 1, 2) $ is the fractional-order orders, $ 0 < \alpha_i < 1 $ and $ \frac{\partial^{\alpha_1}P}{\partial t^{\alpha_1}}, \frac{\partial^{\alpha_2}Z}{\partial t^{\alpha_2}} $ are the Grünwald-Letnikov (GL) fractional derivative. $ t $ is the time variable. The meanings of the parameters in the model are shown in Table 1.

Table 1. Parameter and biological significance in the model.

Parameter	Ecological Significance
$ r $	The intrinsic growth rate of phytoplankton
$ K $	The carrying capacity of the environment for phytoplankton
$ h_{1} $	The time zooplankton need to process one unit of phytoplankton
$ h_{2} $	The time zooplankton need to process one unit of supplementary food
$ e_{1} $	The capacity of zooplankton to prey on phytoplankton
$ e_{2} $	The capacity of zooplankton to forage for supplementary food
$ n_{1} $	The nutritional value of phytoplankton
$ n_{2} $	The nutritional value of supplementary food
$ A $	The biomass of supplementary food
$ \tau $	The mortality rate of zooplankton
$ \sigma $	The proportion of toxins released by phytoplankton
$ h $	The half-saturation constant

| Show Table

DownLoad: CSV

Several numerical schemes have been used to solve the Fractional-order reaction-diffusion equation, the Galerkin method ^{[21, 22]}, the Fourier spectral method ^{[23, 24]}, the finite difference method ^{[25, 26]}, the operational matrix method ^{[27, 28]}, and so on ^{[29, 30, 31, 32, 33, 34]}. In this paper, we investigate dynamic properties and numerical solutions of the fractional-order PZEM (1.2). The major contributions of this paper are as follows:

(a) The concept of the phytoplankton-zooplankton ecological model is first extended to fractional-order PZEM that incorporates the effects of toxic substances and additional food transmission in the environment. This model leverages the memory effect of fractional-order derivatives to more effectively capture biological significance compared to traditional integer-order models.

(b) Stability, Turing instability, Hopf bifurcation, and weakly nonlinear property are analyzed for the PZEM. a new high-precision numerical method is developed for the fractional PZEM without diffusion term, and a discretization method is established for the PZEM with a diffusion term.

The paper is organized as follows: In Section 2, we give the stability and Hopf bifurcation analysis. In Section 3, we detail weakly nonlinear analysis. In Section 4, we give the numerical algorithm and the numerical simulation results of the system. Finally, in Section 5, we the summarize paper.

2. Stability analysis

In this section, we investigate the stability of the equilibrium point and the emergence of Hopf bifurcation. Let

$\begin{eqnarray*} \begin{aligned} &\bar{u} = e_1 h_1 P, \quad \bar{v} = \frac{e_1}{r} Z, \quad t = r T, \quad \kappa = e_1 h_1 K, \quad \alpha = \frac{n_1 h_2}{n_2 h_1}, \\ & \tau = \frac{n_1}{r h_1}, \quad \zeta = \frac{n_2 e_2 h_1}{n_1} A, \quad m = \frac{\tau}{r}, \quad \mu = \frac{\tau}{r}, \quad G = e_1 h h_1. \end{aligned} \end{eqnarray*}$ $

Equation (1.2) is converted to the following form:

$\begin{eqnarray} \left\{ \begin{array}{l} \frac{\partial^{\alpha_1}}{\partial t^{\alpha_1}} \bar{u} = \bar{u}\left(1 - \frac{\bar{u}}{\kappa}\right) - \frac{\bar{u}\bar{v}}{1 + \alpha\zeta + \bar{u}} + d_{1} \nabla^{2} \bar{u}, \\ \frac{\partial^{\alpha_2}}{\partial t^{\alpha_2}} \bar{v} = \frac{\tau(\bar{u} + \zeta) \bar{v}}{1 + \alpha\zeta + \bar{u}} - m \bar{v} - \mu \frac{\bar{u}\bar{v}}{G + \bar{u}} + d_{2} \nabla^{2} \bar{u} + d_{3} \nabla^{2} \bar{v}, \end{array} \right. \end{eqnarray}$ $

(2.1)

where, $ d_1 = e_1h_1d_{11}, d_2 = e_1h_1d_{21}, d_3 = \frac{e_1}{r}d_{31} $.

2.1. Stability analysis

We analyze the stability of Eq (2.1) without the diffusion term.

$\begin{eqnarray} \left\{ \begin{array}{l} \frac{\partial^{\alpha_1}}{\partial t^{\alpha_1}} \bar{u} = \bar{u}\left(1 - \frac{\bar{u}}{\kappa}\right) - \frac{\bar{u}\bar{v}}{1 + \alpha\zeta + \bar{u}}, \\ \frac{\partial^{\alpha_2}}{\partial t^{\alpha_2}} \bar{v} = \frac{\tau(\bar{u} + \zeta) \bar{v}}{1 + \alpha\zeta + \bar{u}} - m \bar{v} - \mu \frac{\bar{u}\bar{v}}{G + \bar{u}}. \end{array} \right. \end{eqnarray}$ $

(2.2)

Solving the following equation:

$\begin{eqnarray} \left\{ \begin{array}{l} \bar{u}\left(1 - \frac{\bar{u}}{\kappa}\right) - \frac{\bar{u} \bar{v}}{1 + \alpha \zeta + \bar{u}} = 0 , \\ \frac{\tau (\bar{u} + \zeta) \bar{v}}{1 + \alpha \zeta + \bar{u}} - m \bar{v} - \mu \frac{\bar{u} \bar{v}}{G + \bar{u}} = 0. \end{array} \right. \end{eqnarray}$ $

(2.3)

We can get the following equilibrium points of Eq (2.2):

$\begin{eqnarray} E_{0} = (0, 0), \quad E_{1} = (\kappa, 0), \quad E^{*} = \left(\bar{u}^{*}, \bar{v}^{*}\right), \end{eqnarray}$ $

(2.4)

where

$\begin{eqnarray*} \bar{v}^* = \left(1 + \alpha \zeta + \bar{u}^*\right) \left(1 - \frac{\bar{u}^*}{\kappa}\right). \end{eqnarray*}$ $

$ u^{*} $ represents the positive solution to the quadratic equation (2.5)

$\begin{eqnarray} \bar{D}_{1} \bar{u}^{2} + \bar{D}_{2} \bar{u} + \bar{D}_{3} = 0, \end{eqnarray}$ $

(2.5)

where

$\begin{eqnarray*} \bar{D}_{1} = \tau - m - \mu, \bar{D}_{3} = \tau \zeta G - m G - m \alpha \zeta G, \bar{D}_{2} = \tau \zeta + \tau G - m \alpha \zeta - m G - \mu \alpha \zeta - \mu - m. \end{eqnarray*}$ $

Solving Eq (2.5), we can get:

$\begin{eqnarray} \left\{ \begin{array}{l} \bar{u}_{1, 2}^{*} = \frac{-(\tau \zeta + \tau H - m \alpha \zeta - mH - \mu \alpha \zeta - \mu - m)}{2(\tau - m - \mu)} \pm \frac{\sqrt{(\tau \zeta + \tau H - m \alpha \zeta - mH - \mu \alpha \zeta - \mu - m)^{2} - 4(\tau - m - \mu)(\tau \zeta H - mH - m \alpha \zeta H)}}{2(\tau - m - \mu)}, \\ \\ \bar{v}^{*} = \left(1 + \alpha \zeta + \bar{u}^{*}\right) \left(1 - \frac{\bar{u}^{*}}{\kappa}\right). \end{array} \right. \end{eqnarray}$ $

(2.6)

The Jacobian matrix associated with $ E^{*} = \left(\bar{u}^{*}, \bar{v}^{*}\right) $ is formulated as follows

$\begin{eqnarray} J = \left[ \begin{matrix} a_{11}& a_{12}\\ a_{21}& a_{22}\\ \end{matrix} \right] , \end{eqnarray}$ $

(2.7)

where

$\begin{gather*} a_{11} = 1 - \frac{2 \bar{u}^*}{\kappa} - \frac{\bar{v}^* (1 + \alpha \zeta)}{(1 + \alpha \zeta + \bar{u}^*)^2}, \ a_{12} = -\frac{\bar{u}^*}{1 + \alpha \zeta + \bar{u}^*}, \\ a_{21} = \frac{\tau \bar{v}^* (1 + \alpha \zeta - \zeta)}{(1 + \alpha \zeta + \bar{u}^*)^2} - \frac{\mu \bar{v}^* G}{(G + \bar{u}^*)^2}, \ a_{22} = \frac{\tau (\bar{u}^* + \zeta)}{1 + \alpha \zeta + \bar{u}^*} - m - \frac{\mu \bar{u}^*}{G + \bar{u}^*}. \end{gather*}$ $

When $ (\alpha_{1}, \alpha_{2}) = (1, 1) $, the characteristic equation at the equilibrium point is as follows:

$\begin{eqnarray} \lambda ^2+Tr_{0}\lambda +\bigtriangleup_{0} = 0, \end{eqnarray}$ $

(2.8)

where,

$\begin{eqnarray*} Tr_{0} = a_{11}+a_{22}, \bigtriangleup_{0} = a_{11}a_{22}-a_{21}a_{12}. \end{eqnarray*}$ $

Let $ \alpha = 0.55, \tau = 4.1, \zeta = 0.72, G = 8.43, \mu = 0.79, \kappa = 2.38 $, and $ m = 2.45 $. The system's equilibrium point and its associated eigenvalues are presented in Table 2.

Table 2. The calculation results of the equilibrium point, eigenvalue, and argument of the system.

Equilibrium point		Eigenvalues of Jacobian matrix		Stability
0	0	–0.3354	1	Unstability
2.3800	0	–1	0.7421	Unstability
0.3131	1.4843	0.0138–0.4838i	0.0138+0.4838i	Unstability

| Show Table

DownLoad: CSV

It can be seen from Table 3 that all three equilibrium points are unstable points, and the system may generate chaos. Next, we introduce a high-precision numerical method for numerical simulation.

Table 3. Comparison of methods for problem (2.9) when $ \beta = 0.9 $ and $ t \in [0, 1]$.

Method	$ \alpha=0.4 $		$ \alpha=0.7 $		$ \alpha=1 $
Method	$ h=0.01 $	$ h=0.05 $	$ h=0.01 $	$ h=0.05 $	$ h=0.01 $	$ h=0.05 $
Closed-from solution ^[35]	3.2307e-02	5.4302e-02	4.5734e-02	2.5713e-01	2.0716e-01	4.5136e-01
Present method	8.9463e-05	5.6724e-04	9.1684e-05	7.3524e-04	2.1527e-04	4.5136e-03

| Show Table

DownLoad: CSV

2.2. Numerical comparison

In simulations of this section, let $ \alpha = 1.1, \tau = 4.4, \zeta = 0.8, G = 8, n = 1.75, \kappa = 5.3 $, and $ m = 2.5 $ and initial conditions $ u_0 = 1 $ and $ v_0 = 1 $. We set the time step size $ h = 0.01 $ and the time as $ T = 600 $ for simulation. Figure 1 shows numerical results at $ (\alpha_1, \alpha_2) = (0.895, 0.994) $.

Figure 1. Comparative numerical result of the two methods of the model $ \left(\alpha, \tau, \zeta, G, n, \kappa, m\right) = \left(1.1, 4.4, 0.8, 8, 1.75, 5.3, 2.5\right) $, $ \alpha_1 = 0.895 $, and $ \alpha_2 = 0.994 $.

DownLoad: Full-Size Img PowerPoint

We consider the following fractional order GL initial value problem

$\begin{equation} \left\{ \begin{aligned} D_t^\alpha y(t) & = t^\beta, \\ y(0) & = 0. \end{aligned} \right. \end{equation}$ $

(2.9)

Numerical simulations are performed for Eq (2.9). Compared with the closed-from solution, the simulation results of the high-precision numerical method are in good agreement with the simulation results of other methods, and the accuracy is higher. This verifies the effectiveness of the high-precision numerical method.

2.3. Brief description of the numerical methodology

We first consider the following GL differential equation:

$\begin{eqnarray} \frac{d^{\alpha}}{dt^{\alpha}} y\left( t \right) = f\left( t, y\left( t \right) \right). \end{eqnarray}$ $

(2.10)

According to ^[36], we can get the high-precision numerical method ^{[36, 37, 38, 39]}:

$\begin{eqnarray} y\left( t_k \right) = h^{\alpha}f\left( t_k, y\left( t_{k-1} \right) \right) -\sum\limits_{i = 1}^m{\vartheta_{j}^{\left( \alpha, p \right)} y\left( t_{k-j} \right)}, \end{eqnarray}$ $

(2.11)

where,

$\begin{eqnarray} \begin{array}{ll} \left\{ \begin{array}{l} \vartheta^{(\alpha, p)}_0 = \mathbf{g}_0, k = 0, \\ \vartheta^{(\alpha, p)}_k = -\frac{1}{\mathbf{g}_0}\sum\limits_{i = 0}^{k-1}(1-i\frac{1+\alpha}{k})\mathbf{g}_i \vartheta^{(\alpha, p)}_{k-i}, k = 1, 2, ...., p-1, \\ \vartheta^{(\alpha, p)}_k = -\frac{1}{\mathbf{g}_0}\sum\limits_{i = 0}^{p}(1-i\frac{1+\alpha}{k})\mathbf{g}_i \vartheta^{(\alpha, p)}_{k-i}, k = p, p+1, p+2, .... \end{array} \right.\end{array} \end{eqnarray}$ $

(2.12)

$\begin{eqnarray} \begin{array}{ll} \begin{pmatrix} \mathbf{g}_0 \\ \mathbf{g}_1\\ \mathbf{g}_2\\ \vdots\\ \mathbf{g}_p\\ \end{pmatrix} = - \begin{pmatrix} 1& 1 & 1& \cdots & 1 \\ 1 & 2 & 3 & \cdots & p+1\\ 1 & 2^2 & 3^2 & \cdots &(p+1)^2\\ \vdots & \vdots & \vdots & \cdots & \vdots\\ 1 & 2^p& 3^p & \cdots& (p+1)^p\\ \end{pmatrix}^{-1} \begin{pmatrix} 0 \\ 1\\ 2\\ \vdots\\ p\\ \end{pmatrix}. \end{array} \end{eqnarray}$ $

(2.13)

Where $ m = \left[t_k/h \right] +1 $, $ h $ is the step size. Therefore, a high-precision numerical method for the system (2.2) is the following

$\begin{align} \left\{ \begin{array}{l} \bar{u}_{k} = h^{\alpha_1} \left(\bar{u}_{k-1}\left(1 - \frac{\bar{u}_{k-1}}{\kappa}\right) - \frac{\bar{u}_{k-1} \bar{v}_{k-1}}{1 + \alpha \zeta + \bar{u}_{k-1}} \right) + \bar{u}_0 - \sum\limits_{j = 1}^{m} \vartheta_j^{(\alpha_1, p)} \bar{u}_{k-j}, \\ \bar{v}_{k} = h^{\alpha_2} \left( \frac{\tau (\bar{u}_{k} + \zeta) \bar{v}_{k-1}}{1 + \alpha \zeta + \bar{u}_{k}} - m \bar{v}_{k-1} - \mu \frac{\bar{u}_{k} \bar{v}_{k-1}}{G + \bar{u}_{k}} \right) + \bar{v}_0 - \sum\limits_{j = 1}^{m} \vartheta_j^{(\alpha_2, p)} \bar{v}_{k-j}. \end{array} \right. \end{align}$ $

(2.14)

According to ^{[36, 40]}, the least common order of the system is identified as $ 0.9819 $. Let $ \left(\alpha _1, \alpha _2 \right) = \left(1.05, 0.95 \right), \alpha = 0.55, \tau = 4.1, \zeta = 0.72, G = 8.43, \mu = 0.79, \kappa = 2.38 $, and $ m = 2.45 $, and the characteristic equation of the system is given by

$\begin{eqnarray} \lambda ^{40}-\frac{10667}{28823}\lambda ^{21}-\frac{24789}{90071}\lambda ^{19}+\frac{1039692}{4437077} = 0, \end{eqnarray}$ $

(2.15)

with characteristic roots $ \lambda _{1, 2} = 1.0358\pm 0.0439i $, as $ |arg\left(\lambda _{1, 2} \right) | = 0.0423 < \frac{\pi}{40} $. This indicates that the instability of the system may generate chaos. It can be seen from Figure 2 that the system generates chaos. The chaotic phenomenon of the system gradually disappears over time and tends to a stable state.

Figure 2. Numerical simulation of phase diagrams and time series plot of system (2.2) at $ (\alpha_{1}, \alpha_{2}) = (1.05, 0.95). $.

DownLoad: Full-Size Img PowerPoint

We set different orders $ \left(\alpha _1, \alpha _2 \right) = \left(1.1, 0.95 \right) $ to observe the influence of the order on the dynamic behavior of the system at parameters $ \alpha = 0.55, \tau = 4.1, \zeta = 0.72, G = 8.43, n = 0.79, \; \kappa = 2.38 $, and $ m = 2.45 $. The characteristic equation of the system is given by

$\begin{eqnarray} \lambda ^{41}-\frac{10667}{28823}\lambda ^{22}-\frac{24789}{90071}\lambda ^{19}+\frac{1039692}{4437077} = 0, \end{eqnarray}$ $

(2.16)

with characteristic roots $ \lambda _{1, 2} = 1.0021\pm 0.2687i $, as $ |arg\left(\lambda _{1, 2} \right) | = 0.2620 < \frac{\pi}{40} $. It indicates that the instability of the system may generate chaos. It can be seen from Figure 3 that the system is in a quasi-periodic state.

Figure 3. Numerical simulation of phase diagrams and time series plot of system (2.2) at $ (\alpha_{1}, \alpha_{2}) = (1.1, 0.95). $.

DownLoad: Full-Size Img PowerPoint

To gain a clearer view of the system's dynamic behavior, Figure 4 shows the phase diagram of the system under different $ \alpha $ orders.

Figure 4. Comparison of the phytoplankton-zooplankton phase diagram of the system (2.2) at different fractional-order derivatives.

DownLoad: Full-Size Img PowerPoint

2.4. Hopf bifurcation analysis

Next, we explore the potential for a Hopf bifurcation at $ E^{*} $. A Hopf bifurcation is critical as it marks a transition in the model's stability and the emergence of periodic solutions. To achieve this, the characteristic roots of system (2.2) must be purely complex.

The solution to system (2.2) is obtained as follows:

$\begin{eqnarray} \lambda _{1, 2} = \frac{Tr_0\pm \sqrt{Tr_{0}^{2}-4det_0}}{2}. \end{eqnarray}$ $

(2.17)

When $ \alpha = 1 $, system (2.2) undergoes a destabilizing Hopf bifurcation under the conditions $ tr_0 = 0 $ and $ det_0 > 0 $. Given that the stability of system (2.2) is influenced by the fractional derivative, this derivative can be considered a parameter in the Hopf bifurcation. We now determine the conditions for the Hopf bifurcation for system (2.2) around parameter $ E_* $, $ \alpha = \alpha_h $:

(1) At the equilibrium point $ E_* $, the Jacobian matrix possesses a pair of complex conjugate eigenvalues $ \lambda _{1, 2} = a_i+ib_i $, which transition to being purely imaginary at $ \alpha = \alpha_h $;

(2) $ m\left(\alpha _h \right) = 0 $ where $ m\left(\alpha \right) = \alpha \frac{\pi}{2}-\underset{1\le i\le 2}{\min}|arg\left(\lambda _i \right) | $;

(3) $ \left. \frac{\partial m\left(\alpha \right)}{\partial \alpha} \right|_{\alpha = \alpha _h}\ne 0 $.

We now demonstrate that $ E_* $ experiences a Hopf bifurcation as $ \alpha $ crosses the value $ \alpha_h $.

Proof. Given that $ tr_{0}^{2}-4\det _0 < 0 $ and $ tr_0 > 0 $, the eigenvalues form a complex conjugate pair with a positive real component. Thus,

$\begin{eqnarray} 0 < arg\left( \lambda _{12} \right) = \tan ^{-1}\left( \frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right) < \frac{\pi}{2}, \end{eqnarray}$ $

(2.18)

and $ \alpha \frac{\pi}{2} > \Big{|}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right) \Big{|} $ for some $ \alpha $. Let $ \alpha_h \frac{\pi}{2} = \Big{|}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right) \Big{|} $, get $ \alpha _h = \frac{2}{\pi}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right) $. Moreover, $ \left. \frac{\partial m\left(\alpha \right)}{\partial \alpha} \right|_{\alpha = \alpha _h} = \frac{\pi}{2}\ne 0 $. Therefore, all Hopf conditions satisfy. □

3. Weakly nonlinear analysis

Above, we analyze the dynamical behavior of fractional system (2.2). The expression of the solution of system (2.1) involves more complex special functions, such as the Mittag-Leffler function and other hypergeometric functions. Although the form of the solution exists, its expression is very large. Therefore, in this section, we reveal various spatiotemporal behaviors around the Turing bifurcation threshold $ d _2 $ and derive the amplitude equations associated with system (2.1) at $ \alpha_1 = \alpha_2 = 1 $ using multi-scale and weakly nonlinear analysis. The solution of system (2.1) is written in the following form:

$\begin{eqnarray} \left( \begin{array}{l} \bar{u}\\ \bar{v}\\ \end{array} \right) = \left( \begin{array}{c} \bar{u}^*\\ \bar{v}^*\\ \end{array} \right) +\sum\limits_{j = 1}^3{\left( \begin{array}{l} A_{j}^{\bar{u}^*}\\ A_{j}^{\bar{v}^*}\\ \end{array} \right)}e^{iq_j\cdot r}+\, \, \text{c.c}., \end{eqnarray}$ $

(3.1)

where $ \left(A_{j}^{\bar{u}^*}, A_{j}^{\bar{v}^*} \right) ^T $ represents the magnitude of the wave vector $ \boldsymbol{q}_j $, fulfilling the condition that $ \left|\boldsymbol{q}_j \right| = q_T $.

Adding a perturbation to the equilibrium point $ E^{*}(u^*, v^*) $, such that $ \bar{u} = \bar{u}^*+u $, $ \bar{v} = \bar{v}^*+v $, and substituting them into system (2.1), and performing a Taylor expansion, we can get the following form:

$\begin{eqnarray} \left\{ \begin{array}{l} \frac{\partial u}{\partial t} = a_{11}u+a_{12}v+k_{20}u^2+k_{02}v^2+k_{11}uv+k_{30}u^3+k_{03}v^3+k_{21}u^2v+k_{12}uv^2+d_1\nabla ^2u, \\ \frac{\partial v}{\partial t} = a_{11}u+a_{22}v+m_{20}u^2+m_{02}v^2+m_{11}uv+m_{30}u^3+m_{03}v^3+m_{21}u^2v+m_{12}uv^2+d_2\nabla ^2u+d_3\nabla ^2v, \end{array} \right. \end{eqnarray}$ $

(3.2)

where

$\begin{gather*} k_{02} = 0, \quad k_{03} = 0, \quad k_{12} = 0, \quad m_{02} = 0, \quad m_{03} = 0, \quad m_{12} = 0, \\ k_{11} = -\frac{1 + \alpha \zeta}{(1 + \alpha \zeta + \bar{u}^*)^2}, \quad m_{11} = \frac{(1 + \alpha \zeta - \zeta) \tau}{(1 + \alpha \zeta + \bar{u}^*)^2} - \frac{\mu G}{(G + \bar{u}^*)^2}, \\ k_{21} = -\frac{1 + \alpha \zeta}{(1 + \alpha \zeta + \bar{u}^*)^3}, \quad m_{21} = \frac{-\tau (1 + \alpha \zeta - \zeta)}{(1 + \alpha \zeta + \bar{u}^*)^3} + \frac{\mu G}{(G + \bar{u}^*)^3}, \\ k_{30} = -\frac{(1 + \alpha \zeta) \left(1 - \frac{\bar{u}^*}{\kappa}\right)}{(1 + \alpha \zeta + \bar{u}^*)^3}, \quad m_{30} = \frac{\tau \left(1 - \frac{\bar{u}^*}{\kappa}\right) (1 + \alpha \zeta - \zeta)}{(1 + \alpha \zeta + \bar{u}^*)^3} - \frac{\mu G \bar{v}^*}{(G + \bar{u}^*)^4}, \\ k_{20} = -\frac{1}{\kappa} + \frac{(1 + \alpha \zeta) \left(1 - \frac{\bar{u}^*}{\kappa}\right)}{(1 + \alpha \zeta + \bar{u}^*)^2}, \quad m_{20} = \frac{-\tau (1 + \alpha \zeta - \zeta) \left(1 - \frac{\bar{u}^*}{\kappa}\right)}{(1 + \alpha \zeta + \bar{u}^*)^2} + \frac{\mu G \bar{v}^*}{(G + \bar{u}^*)^3}. \end{gather*}$ $

Let $ U = \left(u, v \right) ^T $ and system (3.2) can be simplified as

$\begin{eqnarray} \frac{\partial U}{\partial t} = LU+N\left( U, U \right), \end{eqnarray}$ $

(3.3)

$\begin{eqnarray*} L = \left( \begin{matrix} a_{11}+d_1\nabla^{2}& a_{12}\\ a_{21}+d_2\nabla^{2}& a_{22}+d_3\nabla^{2}\\ \end{matrix} \right), \end{eqnarray*}$ $

where

$\begin{eqnarray*} N = \left(\begin{array}{l} k_{20} u^2 + k_{02} v^2 + k_{11} u v + k_{30} u^3 + k_{03} v^3 + k_{21} u^2 v + k_{12} uv^2 \\ m_{20} u^2 + m_{02} v^2 + m_{11} uc + m_{30} u^3 + m_{03} v^3 + m_{21} u^2 v + m_{12} uv^2 \end{array}\right). \end{eqnarray*}$ $

Here, $ L $ denotes a linear operator and $ N $ represents a nonlinear operator.

We expand parameter $ d _2 $ with a sufficiently small parameter $ \varepsilon $, so we can obtain:

$\begin{eqnarray} d^{T}_{2} = \varepsilon d^{(1)} _2+\varepsilon ^2d^{(2)} _2+\varepsilon ^3d^{(3)} _2+\mathcal{O}\left(\varepsilon ^3 \right). \end{eqnarray}$ $

(3.4)

We develop variable $ U $ and nonlinear term $ N $, respectively, the small parameter $ \varepsilon $:

$\begin{eqnarray} U = \left( \begin{array}{c} u\\ v\\ \end{array} \right) = \varepsilon \left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) +\varepsilon ^2\left( \begin{array}{c} u_2\\ v_2\\ \end{array} \right) +\varepsilon ^3\left( \begin{array}{c} u_3\\ v_3\\ \end{array} \right) +\mathcal{O}\left(\varepsilon ^3 \right), \end{eqnarray}$ $

(3.5)

$\begin{eqnarray} N = \varepsilon ^2N_2+\varepsilon ^3N_3+\mathcal{O}\left(\varepsilon ^4 \right), \end{eqnarray}$ $

(3.6)

where

$\begin{gather*} N_2 = \left( \begin{array}{c} k_{20}u_{1}^{2}+k_{11}u_1v_1\\ m_{20}u_{1}^{2}+m_{11}u_1v_1\\ \end{array} \right), \\ N_3 = \left( \begin{array}{c} 2k_{20}u_1u_2+k_{11}u_1v_2+k_{11}u_2v_1+k_{30}u_{1}^{3}+k_{21}u_{1}^{2}v_1\\ 2m_{20}u_1u_2+m_{11}u_1v_2+m_{11}u_2v_1+m_{30}u_{1}^{3}+m_{21}u_{1}^{2}v_1\\ \end{array} \right). \end{gather*}$ $

We decompose operator $ L $ into the following form:

$\begin{eqnarray} L = L_c+\left( d_{2} -d ^{T}_2 \right) M, \end{eqnarray}$ $

(3.7)

where

$\begin{gather*} L_c = \left( \begin{matrix} d_1\nabla ^2+a_{11}& a_{12}\\ d_2^{T}\nabla ^2+a_{21}& d_3\nabla ^2+a_{22}\\ \end{matrix} \right), M = \left( \begin{matrix} b_{11}& b_{12}\\ b_{21}& b_{22}\\ \end{matrix} \right) . \end{gather*}$ $

Using the multiple-scale approach, we differentiate the system's time scales into $ T_1 = \varepsilon t, T _2 = \varepsilon ^2t, T_3 = \varepsilon ^3t $, which are mutually independent. Consequently, the time derivative can be articulated as:

$\begin{eqnarray} \frac{\partial}{\partial t} = \varepsilon \frac{\partial}{\partial T_1}+\varepsilon ^2\frac{\partial}{\partial T_2}+\varepsilon ^3\frac{\partial}{\partial T_3}+\mathcal{O}\left(\varepsilon ^3 \right). \end{eqnarray}$ $

(3.8)

Substituting Eqs (3.4)–(3.8) into Eq (3.3), we can obtain the following equation:

$\begin{eqnarray} \varepsilon :L_c\left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) = 0, \end{eqnarray}$ $

(3.9)

$\begin{eqnarray} \varepsilon ^2:L_c\left( \begin{array}{c} u_2\\ v_2\\ \end{array} \right) = \frac{\partial}{\partial T_1}\left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) -d^{(1)} _2M\left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) -N_2, \end{eqnarray}$ $

(3.10)

$\begin{eqnarray} \varepsilon ^3:L_c\left( \begin{array}{c} u_3\\ v_3\\ \end{array} \right) = \frac{\partial}{\partial T_1}\left( \begin{array}{c} u_2\\ v_2\\ \end{array} \right) +\frac{\partial}{\partial T_2}\left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) -d^{(1)} _2M\left( \begin{array}{c} u_2\\ v_2\\ \end{array} \right) -d^{(2)} _2M\left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) -N_3. \end{eqnarray}$ $

(3.11)

The $ 1^{st} $ order of $ (\mathcal{O}\left(\varepsilon \right)) $:

$\begin{eqnarray} L_c\left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) = 0. \end{eqnarray}$ $

(3.12)

The solution of Eq (3.12) is

$\begin{eqnarray} \left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) = \left( \begin{array}{c} \varphi\\ 1\\ \end{array} \right) \sum\limits_{j = 1}^3{\left( W_je^{i\boldsymbol{q}_j\cdot r} \right)}+c.c., \end{eqnarray}$ $

(3.13)

where $ \varphi = \frac{-a_{12}}{a_{11}-{d}_{1}q_{c}^{2}}, \left| \boldsymbol{q}_{_j} \right| = q_c, q_c = q_T\left(d ^{T}_2 \right) $. Additionally, $ W_j $ denotes the amplitude of $ e^{i\boldsymbol{q}_j\cdot r} $.

The $ 2^{nd} $ order of $ (\mathcal{O}\left(\varepsilon ^2 \right)) $:

$\begin{eqnarray} L_c\left( \begin{array}{c} u_2\\ v_2\\ \end{array} \right) = \frac{\partial}{\partial T_1}\left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) -d^{(1)} _2M\left( \begin{array}{c} u_1\\ v_1\\ \end{array} \right) -N_2 = \left( \begin{array}{c} F_u\\ F_v\\ \end{array} \right), \end{eqnarray}$ $

(3.14)

$ F_u^j $ and $ F_v^j (j = 1, 2, 3) $ are the coefficients with respect to $ e^{i\boldsymbol{q}_j\cdot r} $ in $ F_u $ and $ F_v $, respectively.

Using the Fredholm solvability condition for Eq (3.14), the vector function on the right-hand side must be perpendicular to the zero eigenvalue of $ L_{c}^{+} $. The zero eigenvector of $ L_{c} $ is:

$\begin{eqnarray} \left( \begin{array}{c} 1\\ \psi\\ \end{array} \right) e^{-iq_jr}+c.c., \ j = 1, 2, 3, \end{eqnarray}$ $

(3.15)

with $ \psi = -\frac{a_{11}-{d}_{1}q_{c}^{2}}{a_{21}-{d}^{T}_{2}q_{c}^{2}} $. From the orthogonal condition:

$\begin{eqnarray} \left( 1, \psi \right) \left( \begin{array}{c} F_{u}^{j}\\ F_{v}^{j}\\ \end{array} \right) = 0, \, \, j = 1, 2, 3. \end{eqnarray}$ $

(3.16)

Utilizing the Fredholm solvability condition as stated in Eq (3.16), we arrive at the following:

$\begin{eqnarray} \left\{ \begin{array}{c} \left( \varphi +\psi \right) \frac{\partial W_1}{\partial T_1} = d^{(1)} _2\eta_0W_1+2\left( \eta_1+\psi \eta_2 \right) \overline{W}_2\overline{W}_3, \\ \left( \varphi +\psi \right) \frac{\partial W_2}{\partial T_1} = d^{(1)} _2\eta_0W_2+2\left( \eta_1+\psi \eta_2 \right) \overline{W}_1\overline{W}_3, \\ \left( \varphi +\psi \right) \frac{\partial W_3}{\partial T_1} = d^{(1)} _2\eta_0W_3+2\left( \eta_1+\psi \eta_2 \right) \overline{W}_1\overline{W}_2, \end{array} \right. \end{eqnarray}$ $

(3.17)

where

$\begin{eqnarray*} \left\{ \begin{array}{l} \eta_0 = \left( \varphi b_{11}+b_{22} \right) +\psi \left( \varphi b_{21}+b_{22} \right) , \\ \eta_1 = \frac{1}{2}k_{20}\varphi ^2+k_{11}\varphi , \\ \eta_2 = \frac{1}{2}m_{20}\varphi ^2+m_{11}\varphi . \end{array} \right. \end{eqnarray*}$ $

Next, higher-order disturbance terms are introduced:

$\begin{eqnarray} \left( \begin{array}{c} u_2\\ v_2\\ \end{array} \right) = \left( \begin{array}{c} U_0\\ V_0\\ \end{array} \right) +\sum\limits_{j = 1}^3{\left( \begin{array}{c} U_j\\ V_j\\ \end{array} \right)}e^{i\boldsymbol{q}_j\cdot r}+\sum\limits_{j = 1}^3{\left( \begin{array}{c} U_{jj}\\ V_{jj}\\ \end{array} \right)}e^{2i\boldsymbol{q}_j\cdot r}+\\ \left( \begin{array}{c} U_{12}\\ V_{12}\\ \end{array} \right) e^{i\left( \boldsymbol{q}_1-\boldsymbol{q}_2 \right) \cdot r}+\left( \begin{array}{c} U_{23}\\ V_{23}\\ \end{array} \right) e^{i\left( \boldsymbol{q}_2-\boldsymbol{q}_3 \right) \cdot r}+\left( \begin{array}{c} U_{31}\\ V_{31}\\ \end{array} \right) e^{i\left( \boldsymbol{q}_3-\boldsymbol{q}_1 \right) \cdot r}+c.c.. \end{eqnarray}$ $

(3.18)

We substitute formulas (3.13) and (3.18) into formula (3.10). By solving the equations corresponding to different modes, it is known that the coefficients have the following forms:

$\begin{eqnarray*} \begin{aligned} \left( \begin{array}{l} U_0\\ V_0\\ \end{array} \right) & = \left( \begin{array}{l} u_{00}\\ v_{00}\\ \end{array} \right) \left( |W_1|^2+|W_2|^2+|W_3|^2 \right) , \ U_j = \varphi Y_j, \\ \left( \begin{array}{l} U_{jj}\\ V_{jj}\\ \end{array} \right) & = \left( \begin{array}{l} u_{11}\\ v_{11}\\ \end{array} \right) W_{j}^{2}, \, \, \left( \begin{array}{l} U_{ij}\\ V_{ij}\\ \end{array} \right) = \left( \begin{array}{l} \bar{u}^{*}\\ \bar{v}^{*}\\ \end{array} \right) W_i\overline{W_j}, \\ \end{aligned} \end{eqnarray*}$ $

where

$\begin{eqnarray} \left(\begin{array}{c} \bar{u}_{00} \\ \bar{v}_{00} \end{array}\right) & = \frac{-2}{a_{11}a_{22} - a_{21} a_{12}} \left(\begin{array}{c} a_{22} \eta_1 - a_{12} \eta_2 \\ a_{11} \eta_2 - a_{21} \eta_1 \end{array}\right), \end{eqnarray}$ $

(3.19)

$\begin{eqnarray} \left(\begin{array}{c} \bar{u}_{11} \\ \bar{v}_{11} \end{array}\right) & = \left(\begin{array}{c} \frac{a_{12} \eta_2 - (a_{22} - 4 d_3 q_c^2) \eta_1}{(a_{11} - 4 d_1 q_c^2)(a_{22} - 4 d_3 q_c^2) - (a_{21} - 4 d ^{T}_2 q_c^2) a_{12}} \\ \frac{(a_{21} - 4 d^{T}_2 q_c^2) \eta_1 - (a_{11} - 4 d_1 q_c^2) \eta_2}{(a_{11} - 4 d_1 q_c^2)(a_{22} - 4 {d}_3 q_c^2) - (a_{21} - 4 d^{T}_2 q_c^2) a_{12}} \end{array}\right), \end{eqnarray}$ $

(3.20)

$\begin{eqnarray} \left(\begin{array}{c} \bar{u}_{*} \\ \bar{v}_{*} \end{array}\right) & = 2 \left(\begin{array}{c} \frac{a_{12} \eta_2 - (a_{22} - 3 d_3 q_c^2) \eta_1}{(a_{11} - 3 d_1 q_c^2)(a_{22} - 3 d_3 q_c^2) - (a_{21} - 3 d^{T}_2 q_c^2) a_{12}} \\ \frac{(a_{21} - 3 d^{T}_2 q_c^2) \eta_1 - (a_{11} - 3 d_1 q_c^2) \eta_2}{(a_{11} - 3 d_1 q_c^2)(a_{22} - 3 d_3 q_c^2) - (a_{21} - 3 d^{T}_2 q_c^2) a_{12}} \end{array}\right). \end{eqnarray}$ $

(3.21)

From the orthogonal condition, we obtain:

$\begin{eqnarray} \left\{ \begin{array}{l} \left( \varphi +\psi \right) \left( \frac{\partial Y_1}{\partial T_1}+\frac{\partial W_1}{\partial T_2} \right) = \eta_0\left( d^{(2)} _2W_1+d^{(1)} _2Y_1 \right) -\left[ \left( Q_1+\psi R_1 \right) \left| W_1 \right|^2+\left( Q_2+\psi R_2 \right) \left( \left| W_2 \right|^2+\left| W_3 \right|^2 \right) \right] W_1\\ \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, +2\left( \eta_1+\psi \eta_2 \right) \left( \overline{W_2}\overline{Y_3}+\overline{W_3}\overline{Y_2} \right) , \\ \left( \varphi +\psi \right) \left( \frac{\partial Y_2}{\partial T_1}+\frac{\partial W_2}{\partial T_2} \right) = \eta_0\left( d^{(2)} _2W_2+d^{(1)} _2Y_2 \right) -\left[ \left( Q_1+\psi R_1 \right) \left| W_2 \right|^2+\left( Q_2+\psi R_2 \right) \left( \left| W_1 \right|^2+\left| W_3 \right|^2 \right) \right] W_2\\ \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, +2\left( \eta_1+\psi \eta_2 \right) \left( \overline{W_1}\overline{Y_3}+\overline{W_3}\overline{Y_1} \right) , \\ \left( \varphi +\psi \right) \left( \frac{\partial Y_3}{\partial T_1}+\frac{\partial W_3}{\partial T_2} \right) = \eta_0\left( d^{(2)} _2W_3+d^{(1)} _2Y_3 \right) -\left[ \left( Q_1+\psi R_1 \right) \left| W_3 \right|^2+\left( Q_2+\psi R_2 \right) \left( \left| W_2 \right|^2+\left| W_1 \right|^2 \right) \right] W_3\\ \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, +2\left( \eta_1+\psi \eta_2 \right) \left( \overline{W_2}\overline{Y_1}+\overline{W_1}\overline{Y_2} \right) , \end{array} \right. \end{eqnarray}$ $

(3.22)

where

$\begin{eqnarray*} Q_1 = -\left( 2\varphi k_{20}+k_{11} \right) \left( u_{00}+u_{11} \right) -\varphi k_{11}\left( v_{00}+v_{11} \right) -3k_{30}\varphi ^3-3k_{21}\varphi ^2, \\ Q_2 = -\left( 2\varphi k_{20}+k_{11} \right) \left( u_{00}+u_* \right) -\varphi k_{11}\left( v_{00}+v_* \right) -6k_{30}\varphi ^3-6k_{21}\varphi ^2, \\ R_1 = -\left( 2\varphi m_{20}+m_{11} \right) \left( u_{00}+u_{11} \right) -\varphi m_{11}\left( v_{00}+v_{11} \right) -3m_{30}\varphi ^3-3m_{21}\varphi ^2, \\ R_2 = -\left( 2\varphi m_{20}+m_{11} \right) \left( u_{00}+u_* \right) -\varphi m_{11}\left( v_{00}+v_* \right) -6m_{30}\varphi ^3-6m_{21}\varphi ^2. \end{eqnarray*}$ $

The amplitude $ A_{j} $ is the coefficient of $ e^{i\boldsymbol{q}_j\cdot r} $ at each level, so

$\begin{eqnarray} A_{j} = \varepsilon W_j+\varepsilon ^2Y_j+\mathcal{O}\left(\varepsilon^3 \right) . \end{eqnarray}$ $

(3.23)

Substituting Eqs (3.17) and (3.22) into Eq (3.23), we get the following amplitude equations:

$\begin{eqnarray} \left\{ \begin{array}{c} \tau _0\frac{\partial Z_1}{\partial t} = \eta Z_1+\chi \bar{Z}_2\bar{Z}_3-\left[ g_1|Z_1|^2+g_2\left( |Z_2|^2+|Z_3|^2 \right) \right] Z_1, \\ \tau _0\frac{\partial Z_2}{\partial t} = \eta Z_2+\chi \bar{Z}_1\bar{Z}_3-\left[ g_1|Z_2|^2+g_2\left( |Z_1|^2+|Z_3|^2 \right) \right] Z_2, \\ \tau _0\frac{\partial Z_3}{\partial t} = \eta Z_3+\chi \bar{Z}_1\bar{Z}_2-\left[ g_1|Z_3|^2+g_2\left( |Z_1|^2+|Z_2|^2 \right) \right] Z_3, \\ \end{array} \right. \end{eqnarray}$ $

(3.24)

where

$\begin{gather*} \eta = \frac{d _2 -d ^{T}_2}{d ^{T}_2}, \\ \tau _0 = \frac{\varphi +\psi}{d ^{T}_2\left[ \left( \varphi b_{11}+b_{12} \right) +\psi \left( \varphi b_{21}+b_{22} \right) \right]}, \\ \chi = \frac{k_{20}\varphi ^2+2k_{11}\varphi +m_{20}\varphi ^2\psi +2m_{11}\varphi \psi}{d ^{T}_2\left[ \left( \varphi b_{11}+b_{12} \right) +\psi \left( \varphi b_{21}+b_{22} \right) \right]}, \\ g_1 = \frac{Q_1+\psi R_1}{d ^{T}_2\left[ \left( \varphi b_{11}+b_{12} \right) +\psi \left( \varphi b_{21}+b_{22} \right) \right]}, \\ g_2 = \frac{Q_2+\psi R_2}{d ^{T}_2\left[ \left( \varphi b_{11}+b_{12} \right) +\psi \left( \varphi b_{21}+b_{22} \right) \right]}. \end{gather*}$ $

Since each amplitude $ A_j = \omega _je^{i\theta _j}\left(j = 1, 2, 3 \right) $ in Eq (3.24) can be decomposed into mode $ \omega _j = |A_j| $ and phase angle $ \theta_j $, substituting $ A_{j} $ into Eq (3.24) to separate the real and imaginary parts yields the following equation:

$\begin{eqnarray} \left\{ \begin{array}{l} \tau _0\frac{\partial \theta}{\partial t} = -\chi \frac{\omega _{1}^{2}\omega _{2}^{2}+\omega _{1}^{2}\omega _{3}^{2}+\omega _{2}^{2}\omega _{3}^{2}}{\omega _1\omega _2\omega _3}\sin \theta, \\ \tau _0\frac{\partial \omega _1}{\partial t} = \eta \omega _1+\chi \omega _2\omega _3\cos \theta -\left[ g_1\omega _{1}^{3}+g_2\left( \omega _{2}^{2}+\omega _{3}^{2} \right) \omega _1 \right], \\ \tau _0\frac{\partial \omega _1}{\partial t} = \eta \omega _1+\chi \omega _2\omega _3\cos \theta -\left[ g_1\omega _{1}^{3}+g_2\left( \omega _{2}^{2}+\omega _{3}^{2} \right) \omega _1 \right] , \\ \tau _0\frac{\partial \omega _1}{\partial t} = \eta \omega _1+\chi \omega _2\omega _3\cos \theta -\left[ g_1\omega _{1}^{3}+g_2\left( \omega _{2}^{2}+\omega _{3}^{2} \right) \omega _1 \right], \\ \end{array} \right. \end{eqnarray}$ $

(3.25)

where, $ \theta = \theta _1+\theta _2+\theta _3 $. From Eq (3.25), it can be deduced that the solution is stable under the conditions $ \chi > 0, \psi = 0 $ and $ \chi < 0, \psi = \pi $. The solutions of Eq (3.25) are shown in Table 4.

Table 4. The relationship between pattern shape and steady-state solutions.

Conditions	Solution	Pattern shape
—	$ \omega_1=\omega_{2}=\omega_3=0 $	Stationary state
—	$ \omega_1=\sqrt{\dfrac{\eta}{g_1}}, \omega_2=\omega_3=0 $	Strip pattern
$ \eta > \eta_1=\dfrac{-\chi^2}{4(g_1+2g_2)} $	$ \omega_1=\omega_{2}=\omega_3=\dfrac{\|\chi\|\pm\sqrt{\chi^2+4(g_1+2g_2)\eta}}{2(g_1+2g_2)} $	Hexagon pattern
$ g_2 > g_1, \eta > g_1\omega^2_1 $	$ \omega_1=\dfrac{\|\chi\|}{g_2-g_1}, \omega_2=\omega_3=\sqrt{\dfrac{\eta-g_1\omega^2_1}{g_1+g_2}} $	Mixed state

| Show Table

DownLoad: CSV

4. Numerical simulation

In this section, we conduct numerical simulations to verify the theoretical analysis and observe its pattern dynamic behavior. In order to observe the dynamic behavior, we introduce the numerical method of the PZEM model (2.1). If $ 0 < \alpha_i < 1 (i = 1, 2) $, the time derivative is expressed by formula (2.11), and the space derivative is expressed by the following formula (4.2). If $ \alpha_i = 1(i = 1, 2) $, we employ the Euler discretization approach for conducting numerical simulations in a two-dimensional domain $ \varOmega = \left[0, L_x \right] \times \left[0, L_y \right] $. We select $ L_x = 250, L_y = 250 $, a time increment $ \varDelta t = 0.2 $, and a spatial increment $ \varDelta h = 0.69 $. We denote $ u_{pq}^{n} = u\left(x_p, y_q, n\varDelta t \right) $ and $ v_{pq}^{n} = v\left(x_p, y_q, n\varDelta t \right) $. The model (2.1) is discretized using the Euler method, as outlined below:

$\begin{eqnarray} \left\{ \begin{array}{l} \frac{\bar{u}_{pq}^{n+1} - \bar{u}_{pq}^{n}}{\Delta t} = \bar{u}_{pq}^{n} \left(1 - \frac{\bar{u}_{pq}^{n}}{\kappa}\right) - \frac{\bar{u}_{pq}^{n} \bar{v}_{pq}^{n}}{1 + \alpha \zeta + \bar{u}_{pq}^{n}} + d_{1} \nabla^{2} \bar{u}_{pq}^{n}, \\ \frac{\bar{v}_{pq}^{n+1} - \bar{v}_{pq}^{n}}{\Delta t} = \frac{\tau \left(\bar{u}_{pq}^{n} + \zeta\right) \bar{v}_{pq}^{n}}{1 + \alpha \zeta + \bar{u}_{pq}^{n}} - m \bar{v}_{pq}^{n} - \mu \frac{\bar{u}_{pq}^{n} \bar{v}_{pq}^{n}}{G + \bar{u}_{pq}^{n}} + d_{2} \nabla^{2} \bar{u}_{pq}^{n} + d_{3} \nabla^{2} \bar{v}_{pq}^{n}, \end{array} \right. \end{eqnarray}$ $

(4.1)

where

$\begin{eqnarray} \left\{ \begin{array}{l} \nabla^{2} \bar{u}_{pq} = \frac{\bar{u}_{p+1, q+1} + \bar{u}_{p-1, q-1} + \bar{u}_{p+1, q-1} + \bar{u}_{p-1, q+1} + 4 \left( \bar{u}_{p+1, q} + \bar{u}_{p-1, q} + \bar{u}_{p, q+1} + \bar{u}_{p, q-1} \right) - 20 \bar{u}_{pq}}{6h^{2}}, \\ \nabla^{2} \bar{v}_{pq} = \frac{\bar{v}_{p+1, q+1} + \bar{v}_{p-1, q-1} + \bar{v}_{p+1, q-1} + \bar{v}_{p-1, q+1} + 4 \left( \bar{v}_{p+1, q} + \bar{v}_{p-1, q} + \bar{v}_{p, q+1} + \bar{v}_{p, q-1} \right) - 20 \bar{v}_{pq}}{6h^{2}}. \end{array} \right. \end{eqnarray}$ $

(4.2)

Next, we select the parameters shown in Table 5 and use Eq (4.1) to conduct numerical simulations.

Table 5. Parameter values.

$ \alpha $	$ \tau $	$ \zeta $	$ G $	$ n $	$ \kappa $	$ m $	$ \kappa_{c} $	$ d_{1} $	$ d_{2} $	$ d_{3} $
$ 0.55 $	$ 4.1 $	$ 0.72 $	$ 8.43 $	$ 0.79 $	$ 2.38 $	$ 2.45 $	$ 1.91 $	$ 0.0036 $	$ 0.00158 $	$ 0.0141 $

| Show Table

DownLoad: CSV

The parameter values are shown in Table 5, and the calculation results are as follows:

$\begin{gather*} E^{*} = \left( 0.3131, 1.4843 \right) , \chi = -0.4312, \eta = 0.2461, \tau_{0} = 0.3270, g_1 = 1.2869, g_2 = 1.7464. \end{gather*}$ $

The initial conditions are specified as follows:

$\begin{eqnarray*} u\left( x, y, 0 \right) = u_*\left( 1+0.1\left( rand-0.5 \right) \right) , v\left( x, y, 0 \right) = v_*\left( 1+0.1\left( rand-0.5 \right) \right). \end{eqnarray*}$ $

The outcomes of the numerical simulation indicate that speckle and mixed-structure patterns emerge in the graphical representation with this particular set of parameters, as illustrated in Figure 5.

Figure 5. Density distribution pattern of system (2.1) at $ \alpha = 0.55, \tau = 4.1, \zeta = 0.72, G = 8.43, n = 0.79, \kappa = 2.38, m = 2.45, d_{1} = 0.0036, d_{2} = 0.00158, d_{3} = 0.0141 $, and $ \alpha_1 = \alpha_2 = 1 $.

DownLoad: Full-Size Img PowerPoint

Now, we use the parameters in Table 5 and symmetric initial conditions to perform numerical simulations by changing parameters $ d_{1} $, $ d_{2} $, and $ d_{3} $. Numerical simulation was carried out with other parameters remaining unchanged, and the results showed that there is a symmetric mixed structure solution. The subsequent figure illustrates the intricate pattern evolution of the symmetric mixed pattern across parameters. We select the following symmetric initial conditions:

$\begin{eqnarray*} u\left( x, y, 0 \right) = \left\{ \begin{array}{l} u^*, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, x, y\in \left( 80, 120 \right) , \\ u^*-0.001, \, \, other, \\ \end{array} \right. \\ v\left( x, y, 0 \right) = \left\{ \begin{array}{l} v^*, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, x, y\in \left( 80, 120 \right) , \\ v^*-0.001, \, \, other.\\ \end{array} \right. \end{eqnarray*}$ $

Figure 6 shows the spatial distribution patterns of phytoplankton and zooplankton at different diffusion coefficients, $ d_{2} $. The diffusion coefficient $ d $ reflects the influence of phytoplankton on the distribution of zooplankton. When $ d_{2} = 0.003 $, phytoplankton and zooplankton form highly clustered patterns, with zooplankton closely following the prey. When $ d_{2} = 0.005 $, phytoplankton and zooplankton are more dispersed, and the zooplankton pattern is more complex. This indicates that the diffusion ability of phytoplankton affects the hunting efficiency of zooplankton, and zooplankton tend to migrate to areas with a high density of phytoplankton. When $ d_{2} = 0.004 $, the distribution densities of phytoplankton and zooplankton are between the distribution densities of phytoplankton and zooplankton at $ d_{2} = 0.003 $ and $ d_{2} = 0.005 $.

Figure 6. Population density distribution pattern of phytoplankton and zootoplankton with with different parameters $ {d}_{1} = 0.002 $ and $ {d}_{3} = 0.018 $.

DownLoad: Full-Size Img PowerPoint

Let $ d_{1} = 0.002, d_{2} = 0.004 $, and $ d_{3} = 0.015 $ through the observation of the distribution pattern of plankton population density in Figure 7. At $ t = 6000 $, the system presents an incomplete pattern structure. Phytoplankton begins to form an irregular pattern distribution, while zooplankton shows an aggregation pattern following predation. At $ t = 8000 $, the pattern structure gradually forms a stable state. At $ t = 10000 $, the system reaches a stable state and forms a predator aggregation area that is misaligned with the plant pattern. We find that as time goes by, $ t $, the distribution pattern of plankton population density becomes more complex. This indicates that the system presents a complex competitive process over time. This reflects the process by which the ecosystem reaches a stable state through dynamic adjustment.

Figure 7. The distribution patterns of phytoplankton and zooplankton population densities at different time periods.

DownLoad: Full-Size Img PowerPoint

Let parameter $ d_{1}=0.002, d_{2}=0.004 $, and $ d_{3}=0.017 $, the patterns of phytoplankton and zooplankton at different times are shown in Figure 8$ a_{1}-a_{4} $ and Figure 8$ c_{1}-c_{4} $. For diffusion coefficients at $ d_{1}=0.002, d_{2}=0.006 $, and $ d_{3}=0.017 $, the patterns of phytoplankton and zooplankton at different times are shown in Figure 8$ b_{1}-b_{4} $ and Figure 8$ d_{1}-d_{4} $. In Figure 8, the density of phytoplankton gradually increases over time. It can be seen from the pattern that the difference between the high-density area and the low-density area significantly increases. The internal structure of the pattern becomes more complex, showing more details and layers. The distribution pattern of zooplankton shows a significantly misaligned predation aggregation area with that of plants, and the density of predation hotspots continues to increase over time. It can be seen from the pattern that the spatial distribution of phytoplankton and zooplankton shows complex dynamic behaviors.

Figure 8. Population density distribution patterns of phytoplankton and zootoplankton with parameters $ {d}_{1} = 0.002 $ and $ {d}_{3} = 0.017 $.

DownLoad: Full-Size Img PowerPoint

Through the analysis of the patterns of phytoplankton and zooplankton, we see that, as time goes by, the patterns of zooplankton evolve from simple patterns to complex ones, demonstrating the complexity of this ecosystem. The diffusion coefficient is highly sensitive to the generation of patterns. A smaller diffusion coefficient leads to more regular and symmetrical patterns, while a larger diffusion coefficient results in more complex and irregular patterns, significantly affecting the distribution pattern of plankton. The initial conditions and diffusion coefficient significantly impact the long-term ecological behavior of the model. Moreover, initial conditions affect the formation of the pattern, while the parameters affect the distribution structure of the ecosystem.

5. Conclusions

We combine theoretical analysis with numerical simulation to deeply explore the dynamic behavioral characteristics of the fractional-order PZEM. In theoretical analysis, we conduct stability, Turing instability, Hopf bifurcation, and nonlinear analysis for the PZEM. In numerical methods, we develop a high-precision numerical method for fractional-order PZEM without diffusion terms and verify the superiority of this method through comparison with other methods. Moreover, an effective discretization method is established for the model with diffusion terms. The numerical simulation results not only verify the correctness of the theoretical analysis, but also demonstrate complex dynamic behaviors at different diffusion coefficients. Moreover, the system presents significantly different spatial pattern evolution characteristics. The results reveal that the diffusion coefficient is crucial to the generation of patterns. A bigger diffusion coefficient will lead to more complex and irregular pattern structures, which directly affects the spatial distribution patterns of plankton populations. Furthermore, the initial conditions and the diffusion coefficient have a decisive influence on the long-term ecological behavior of the model: The initial conditions affect the formation process of the pattern, while the diffusion coefficient affects the dynamic evolution of the pattern structure. In future work, we will explore the response mechanism of the system at the coupling effect of multiple factors, as well as the effect of dynamic behavior for fractional-order.

Author contributions

Conceptualization, Methodology, Software, Data, Formal analysis and Funding acquisition, Writing-original draft and writing-review and editing: Shuai Zhang, Hao Lu Zhang, Yu Lan Wang and Zhi Yuan Li. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

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

Acknowledgments

This paper is supported by Doctoral research start-up Fund of Inner Mongolia University of Technology (DC2300001252) and the Natural Science Foundation of Inner Mongolia (2024LHMS06025).

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

References

[1]	Evenhuis C, Lenton A, Cantin NE, et al. (2015) Modelling coral calcification accounting for the impacts of coral bleaching and ocean acidification. Biogeosci 12: 2607-2630.
[2]	Guinottea JM, Fabry VJ (2008) Ocean Acidification and Its Potential Effects on Marine Ecosystems. Ann NY Acad Sci 1134: 320-342. doi: 10.1196/annals.1439.013
[3]	Cyronak T, Schulz KG, Jokiel PL (2015) The Omega myth: what really drives lower calcification rates in an acidifying ocean. ICES J Mar Sci 73: 558-562.
[4]	Suwa R, Nakamura M, Morita M, et al. (2010) Effects of acidified seawater on early life stages of scleractinian corals (Genus Acropora). Fish Sci 76: 93-99. doi: 10.1007/s12562-009-0189-7
[5]	Hii YS, Ambok Bolong AM, Yang TT, et al. (2009) Effect of elevated carbon dioxide on two Scleractinian corals: Porites cylindrica (Dana, 1846) and Galaxea fascicularis (Linnaeus, 1767). J Mar Biol 2009: 215196.
[6]	Kerrison P, Hall-Spencer JM, Suggett DJ, et al. (2011) Assessment of pH variability at a coastal CO2 vent for ocean acidification studies. Estuar Coast Shelf Sci 94: 129-137. doi: 10.1016/j.ecss.2011.05.025
[7]	Hall-Spencer JM, Rodolfo-Metalpa R, Martin S, et al. (2008) Volcanic carbon dioxide vents show ecosystem effects of ocean acidification. Nature 454: 96-99. doi: 10.1038/nature07051
[8]	Cigliano M, Gambi MC, Rodolfo-Metalpa R, et al. (2010) Effects of ocean acidification on invertebrate settlement at volcanic CO2 vents. Mar Biol 157: 2489-2502. doi: 10.1007/s00227-010-1513-6
[9]	Johnson VR, Brownlee C, Rickaby REM, et al. (2013) Responses of marine benthic microalgae to elevated CO2. Mar Biol 160: 1813-1824.
[10]	Inoue S, Kayanne H, Yamamoto S, et al. (2013) Spatial community shift from hard to soft corals in acidified water. Nat Clim Chang 3: 683-687. doi: 10.1038/nclimate1855
[11]	Gabay Y, Benayahu Y, Fine M (2013) Does elevated pCO2 affect reef octocorals? Ecol Evol 3: 465-473. doi: 10.1002/ece3.351
[12]	Gabay Y, Fine M, Barkay Z (2014) Octocoral Tissue Provides Protection from Declining Oceanic pH. PloS ONE 9: e91553. doi: 10.1371/journal.pone.0091553
[13]	Michalek-Wagner K, Bourne DJ, Bowden BF (2001) The effects of different strains of zooxanthellae on the secondary-metabolite chemistry and development of the soft-coral host Lobophytum compactum. Mar Biol 138: 753-760. doi: 10.1007/s002270000505
[14]	Changyun W, Haiyan L, Changlun S, et al. (2008) Chemical defensive substances of soft corals and gorgonians. Acta Ecol Sin 28: 2320-2328. doi: 10.1016/S1872-2032(08)60048-7
[15]	Sotka E, Forbey J, Horn M, et al. (2009) The emerging role of pharmacology in understanding consumer-prey interactions in marine and freshwater systems. Integr Comp Biol 49: 291-313. doi: 10.1093/icb/icp049
[16]	Lages BG, Fleury BG, Ferreira CE, et al. (2006) Chemical defense of an exotic coral as invasion strategy. J Exp Mar Biol Ecol 328: 127-135.
[17]	Kahng SE, Grigg RW (2005) Impact of an alien octocoral, Carijoa riisei, on black corals in Hawaii. Coral Reefs 24: 556-562. doi: 10.1007/s00338-005-0026-0
[18]	Aceret TL, Sammarco PW, Coll JC (1995) Toxic effects of alcyonacean diterpenes on scleractinian corals. J Exp Mar Biol Ecol 188: 63-78.
[19]	Sammarco PW, Coll JC, Barre SL (1995). Competitive strategies of soft coral (Coelenterata : Octocorallia), II, variable defensive responses and susceptibility to scleractinian corals. J Exp Mar Biol Ecol 91: 199-215.
[20]	Sammarco PW, Coll JC (1990) Lack of predictability in terpenoid function - multiple roles and integration with related adaptations in soft corals. J Chem Ecol 16: 273-289. doi: 10.1007/BF01021284
[21]	Yang B, Liu J, Wang J, et al. (2015) Cytotoxic Cembrane Diterpenoids. InHandbook of Anticancer Drugs from Marine Origin. Springer International Publishing, 649-672.
[22]	Liu X, Zhang J, Liu Q, et al. (2015) Bioactive Cembranoids from the South China Sea Soft Coral Sarcophyton elegans. Molecules 20: 13324-13335. doi: 10.3390/molecules200713324
[23]	Rocha J, Peixe L, Gomes N, et al. (2011) Cnidarians as a source of new marine bioactive compounds—An overview of the last decade and future steps for bioprospecting. Mar Drugs 9: 1860-1886. doi: 10.3390/md9101860
[24]	Fabricius KE, Langdon C, Uthicke S, et al. (2011) Losers and winners in coral reefs acclimatized to elevated carbon dioxide concentrations. Nat Clim Chang 1: 165-169. doi: 10.1038/nclimate1122
[25]	Pierrot DE, Lewis E, Wallace DWR (2006) MS Exel Program Developed for CO2 System Calculations. ORNL/CDIAC-105a. Oak Ridge, Tennessee, USA: Carbon Dioxide Information Analysis Centre, Oak Ridge National Laboratory, US Department of Energy.
[26]	Fabricius KE, Alderslade P (2001) Soft corals and sea fans: a comprehensive guide to the tropical shallow water genera of the central west Pacific, the Indian Ocean and the Red Sea. Australian Institute of Marine Science, 264.
[27]	Zachary I (2003) Determination of cell number, in: Cell proliferation and apoptosis. D. Hughes and H Mehmet (eds), Bios Scientific Publishers, 13-35.
[28]	Kohler KE, Gill SM (2006) Coral Point Count with Excel extensions (CPCe): A visual basic program for the determination of coral and substrate coverage using random point coral methodology,”. Comput Geosci 32: 1259-1269. doi: 10.1016/j.cageo.2005.11.009
[29]	Hammer O, Harper DAT, Ryan PD (2001) Past: Paleontological Statistics Software package for education and data analysis. Palaeontol Electron 4: 9.
[30]	Anthony KR, Kline DI, Diaz-Pulido G (2008) Acidification causes bleaching and productivity loss in coral reef builders. P Natl Acad Sci USA 105: 17442-17446. doi: 10.1073/pnas.0804478105
[31]	Crook ED, Potts D, Rebolledo-Vieyra M (2012) Calcifying coral abundance near low-pH springs: implications for future ocean acidification. Coral Reefs 31: 239-245. doi: 10.1007/s00338-011-0839-y
[32]	Edmunds PJ (2011) Zooplanktivory ameliorates the effects of ocean acidification on the reef coral Porites sp. Limnol Oceanogr 56: 2402-2410. doi: 10.4319/lo.2011.56.6.2402
[33]	Doney SC, Fabry VJ, Feely RA, et al. (2009) Ocean acidification: the other CO2 problem. Ann Rev Mar Sci 1: 169-192. doi: 10.1146/annurev.marine.010908.163834
[34]	Sammarco PJ, Coll JC (1992) Chemical adaptations in the Octocorallia: evolutionary considerations. Mar Ecol Prog Ser 88: 93-93.
[35]	Luter HM, Duckworth AR (2010) Influence of size and spatial competition on the bioactivity of coral reef sponges. Biochem Syst Ecol 38: 146-153.
[36]	Januar HI, Marraskuranto E, Patantis G, et al. (2012) LC-MS Metabolomic Analysis of Environmental Stressors Impacts to the Metabolites Diversity in Nephthea sp.. Chron Young Sci 2: 57-62.
[37]	Januar HI, Pratitis A, Bramandito A (2015) Will the increasing of anthropogenic pressures reduce the biopotential value of sponges? Scientifica 2015: 734385.
[38]	Januar HI, Chasanah E, Tapiolas DM, et al. (2015) Influence of anthropogenic pressures on the bioactivity potential of sponges and soft corals in the coral reef environment. Squallen Bull Mar Fish Postharvest Biotech 10: 51-59.
[39]	Arnold T, Mealey C, Leahey H, et al. (2012) Ocean Acidification and the Loss of Phenolic Substances in Marine Plants. PLoS ONE 7: e35107.
[40]	Suggett DJ, Hall-Spencer J, Rodofo-Metalpa R, et al. (2012) Sea anemones may thrive in a high CO2 world. Global Chang Biol 18: 3015-3025. doi: 10.1111/j.1365-2486.2012.02767.x

This article has been cited by:

1.	Yangyan Zeng, Yidong Zhou, Wenzhi Cao, Dongbin Hu, Yueping Luo, Haiting Pan, Big data analysis of water quality monitoring results from the Xiang River and an impact analysis of pollution management policies, 2023, 20, 1551-0018, 9443, 10.3934/mbe.2023415
2.	Xiangzheng Fu, Yifan Chen, Sha Tian, DlncRNALoc: A discrete wavelet transform-based model for predicting lncRNA subcellular localization, 2023, 20, 1551-0018, 20648, 10.3934/mbe.2023913
3.	Dongliang Yang, Zhichao Ren, Ge Zheng, The impact of pension insurance types on the health of older adults in China: a study based on the 2018 CHARLS data, 2023, 11, 2296-2565, 10.3389/fpubh.2023.1180024

Reader Comments

Your name:*

Email:*
© 2016 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

AIMS Environmental Science

1.2 2.7

Metrics

Article views(6034) PDF downloads(1223) Cited by(7)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(3) / Tables(2)

AIMS Environmental Science

Changes in soft coral Sarcophyton sp. abundance and cytotoxicity at volcanic CO₂ seeps in Indonesia

Related Papers:

Abstract

1. Introduction

2. Stability analysis

2.1. Stability analysis

2.2. Numerical comparison

2.3. Brief description of the numerical methodology

2.4. Hopf bifurcation analysis

3. Weakly nonlinear analysis

4. Numerical simulation

5. Conclusions

Author contributions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

AIMS Environmental Science

Changes in soft coral Sarcophyton sp. abundance and cytotoxicity at volcanic CO2 seeps in Indonesia

Related Papers:

Abstract

1. Introduction

2. Stability analysis

2.1. Stability analysis

2.2. Numerical comparison

2.3. Brief description of the numerical methodology

2.4. Hopf bifurcation analysis

3. Weakly nonlinear analysis

4. Numerical simulation

5. Conclusions

Author contributions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog

Changes in soft coral Sarcophyton sp. abundance and cytotoxicity at volcanic CO₂ seeps in Indonesia