Bifurcation and pattern dynamics in the nutrient-plankton network

: This paper used a Holling-IV nutrient-plankton model with a network to describe algae’s spatial and temporal distribution and variation in a specific sea area. The stability and bifurcation of the nonlinear dynamic model of harmful algal blooms (HABs) were analyzed using the nonlinear dynamic theory and de-eutrophication’s e ff ect on algae’s nonlinear dynamic behavior. The conditions for equilibrium points (local and global), saddle-node, transcritical, Hopf-Andronov and Bogdanov-Takens (B-T) bifurcation were obtained. The stability of the limit cycle was then judged and the rich and complex phenomenon was obtained by numerical simulations, which revealed the robustness of the nutrient-plankton system by switching between nodes. Also, these results show the relationship between HABs and bifurcation, which has important guiding significance for solving the environmental problems of HABs caused by the abnormal increase of phytoplankton.


Introduction
Differential equation models, such as infectious disease models [1][2][3][4][5][6][7] and neurodynamic models [8,9], are common tools for describing natural and social phenomena changes.They are powerful tools for understanding, explaining, and predicting the behavior of complex systems in the real world.The defensive behavior of a population plays a crucial role in its survival and reproduction of the whole population.It is necessary to further study the intrinsic nature of predation systems with defensive behavior (Holling-IV functional response function).The Holling-type IV predation model has been a popular topic in biomathematics and the research mainly includes bifurcation, equilibria and patterns.Some scholars performed bifurcation analysis in a Holling-type IV predation model and showed that the system undergoes a Codimension-2 Bogdanov-Takens (B-T) bifurcation [10], B-T singularity of Codimension-4 and also multiple other nonhyperbolic and degenerate equilibria, bistability, tristability or even tetrastability [11], as well as the Hopf bifurcation by choosing the level of fear as a bifurcation parameter [12], the Hopf bifurcation of coexisting steady-state for both the delayed and non-delayed systems [13] or supercritical Hopf bifurcation independently of the growth rate of the prey [14].While other researchers obtained the sufficient conditions for the existence of a periodic solution [15], discussed the existence, uniqueness, non-negativity, and boundedness of the solutions [16], or carried out the existence of equilibrium points and stability analysis for correlation model [17], or considered predation model with combined continuous white noise and discontinuous Lvy noise, and proved the global positive solutions' existence and uniqueness [18].Some scholars derived the Turing instability condition for the spatial predation system [19].
Harmful algal blooms (HABs) are ecological anomalies caused by the explosive proliferation or high aggregation of many plankton in seawater under certain environmental conditions.With the rapid development of modern industrial and agricultural production and the increase in population in coastal areas, a lot of industrial and agricultural wastewater and domestic sewage are discharged into the ocean, resulting in an increasing degree of eutrophication in offshore and harbor areas [20].At the same time, the increase in the degree of coastal development and the expansion of mariculture has brought about the pollution of the marine ecological environment and mariculture itself [21].On the other hand, the development of the maritime industry has led to the introduction of exotic and HABs species; global warming has also led to the frequent occurrence of HABs [22].The events of HABs cause significant economic losses and threaten human health [23].
The paper gives the systematic-dynamical analysis of equilibria in section two.In section three, we analyze conditions for several bifurcations (such as transcritical, saddle-node, Hopf and Bogdanov-Takens bifurcation).Finally, we have some discussions in section four.Some interesting results are obtained, which help us better learn the dynamics in the real world.
To better explain the causes of HABs, we next briefly consider the formation mechanism of HABs from the differential equations perspective.We first consider the following equation where x is the concentration of nutrient at time s, y is the density of plankton at time s, r is the conversion rate of inorganic nutrients into biomass, k is the carrying capacity of the nutrient, m > 0 is the mortality of plankton (natural death, predation, etc.), µ is the maximum plankton rate, x a+bx+x 2 is the Holling-IV functional response, a is the half-saturation constant of nutrient concentration growth, b is the denominator of the functional response [24] and c is the effect changes in the nutrient ratio of seawater.Changes in the nutrient ratio of seawater directly change the composition and structure of marine planktonic ecosystems [25].For example, a low nitrogen-to-phosphorus concentration ratio is conducive to Cyanobacteria blooms [26].When the nutrient proportion of silicate in seawater is high, the growth of plankton can be promoted.Diffusion and convection dominate the transfer of HABs from one area to another.HABs, in many parts of the world, are due to the influence of ocean currents, such as G. breve off the coast of West Florida and Diphysis sp off the coast of Sweden [27,28].In addition, the results of population dynamics, especially the nonlinear population dynamics, show that the complex dynamics of the food chain and food web structure can occur in plankton systems.HABs are likely to be a reflection of such complex dynamics.We researched the system with the diffusion network as follows: (1. 2) The description of the network Laplacian matrix L = {L i j } can be found in [29].
Proof.For convenience, we define a function Ψ by Differentiating Ψ with respect to t, together with (1.2), one has Via the theory of differential equations and inequalities, we obtain It is easy to know that the righthand limit of the above equation is Φ m as t → ∞, from which Ψ(t) is bounded.It can be further deduced that the solution of the system is bounded in R 2 + .The system (1.2) has a trivial equilibrium, an axial equilibrium and, at most, two internal equilibrium, recorded as E 0 = (0, 0), , respectively.The trivial equilibrium E 0 always exists.The actual value of the function of nutrients should be positive, so the axial equilibrium E 1 exists when c r > −1, x * 1 and x * 2 are determined by and ).If ∆ < 0, then there is no internal equilibrium; especially if ∆ = 0, then there is one internal equilibrium E * .
We now analyze the equilibrium points' stability by using the stability theory, and the Jacobian matrix J iE (i = 1, ..., n) of the system (1.2) is If c + r < 0, the trivial equilibrium point E 0 is globally asymptotically stable; E 0 is saddle point when c + r > 0.
from which one has tr(J iE * 1 ) = If tr(J iE * 1 ) < 0, then eigenvalues λ 1 and λ 2 are negative and E * 1 is locally asymptotically stable; if tr(J iE * 1 ) > 0, then eigenvalues λ 1 and λ 2 are positive and E * 1 is unstable.It's easy to calculate that det(J iE * 2 ) < 0, then in the existence range of E * 2 , it is always a saddle point.To sum it up theoretically, the stability of the system solution is conditional and the conditions for stability are different between different equilibrium states.There is a phenomenon of overlapping between different equilibrium states.The stability of the equilibrium state E 0 is straightforward.However, for E 1 , although there is a possibility of stability in theory, it is impossible in practice.

Bifurcation of the System
The linear part of (1.2) is ẋi ẏi = a 11 a 12 a 21 a 22 where We expressed the general solution as in [30] where n j L i j ϕ k j = Λ i ϕ k i .Plug (3.2) into (3.1) and the Jacobian matrix B i (i = 1, ..., n) reads ) denotes the eigenvalues of the matrix L.
Proof.When a a T c and d 2 = 0, the Jacobian matrix J iE 1 (i = 1, ..., n) of the system (1.2) is Clearly, det(J iE 1 ) = 0 and J iE 1 has an eigenvalue λ = 0. Let its corresponding eigenvector be p iT and the eigenvector (λ = 0) of the transpose matrix J T iE 1 is q iT , then Therefore, the transcritical bifurcation condition is satisfied.

Saddle-node bifurcation
In section two, we researched the interior equilibrium point's existence when ∆ > 0. The system (1.2) has two distinct internal equilibria when ∆ = 0, i.e. a = (µ−mb) 2 4m 2 a * and two equilibrium points merged into one, which is called a degenerate equilibrium, denoted by E * .Theorem 7 * .When a = a * and d 2 = 0, the system (1.2) exists a saddle-node bifurcation.

Hopf bifurcation
Theorem 8 * .For the system (1.2), we assume that the parameters satisfy the existence condition of c * be such that tr(J E * 1 )| c=c * = 0, then the stability of E * 1 changes and the system (1.2) exists a Hopf bifurcation.
Proof.When c = c * , the Jacobian matrix J iE * and it is easy to verify that We then compute the first Lyapunov coefficient [31].Let c = c * , where the coordinates of the internal equilibrium point are (x,y).Change variables 1 , the Taylor expansion is performed and the system becomes , where p 10 + q 01 = 0, p 10 q 01 − p 01 q 10 > 0, and , q 01 = q 02 = q 12 = q 03 = 0.The expression of the first Lyapunov coefficient is x * 2 1 ) 2 ))}, if lc 1 < 0 (resp.> 0), the limit cycle is stable(resp.unstable).

Bogdanov-Takens bifurcation
Now, we will discuss the Bogdanov-Takens bifurcation of (1.2).Let d 2 = 0 (meaning Plankton at each node do not spread out due to human intervention), then the Jacobian matrix about E * is then J iE * has two zero characteristic roots.Now we will use the method given by Xiao and Ruan in [32] to verify that E * is a cusp of codimensional two and that there exists B-T bifurcation in the neighborhood of the equilibrium point E * .
We select (a, c) as the bifurcation parameters to study the bifurcation of the system (1.2). (a, c) changes in a small neighborhood near (a * , c * ), where (a * , c * ) is the B-T bifurcation point.We consider the following system where (ε a , ε c ) are small parameters.Through the following transformation where 3  , Then, in the small neighborhood of (0, 0), we make the following C ∞ transformation and system (3.4) becomes a 01 , P 3 (X i2 , Y i2 ) is a smooth function and Introducing a new time variable τ that satisfies ds = (1 − p 02 X i2 ), and τ is still recorded as s, then system (3.5)becomes The system (3.6)becomes where P 4 (X i3 , Y i3 ) is a smooth function and ).Note that the sign of p 20 − 2p 02 p 10 + p 00 p 2 02 is difficult to determine.We consider the following two scenarios.Case I: If p 20 − 2p 02 p 10 + p 00 p 2 02 > 0, then we make the following variable transformation The system (4.1)becomes where r 00 = q 00 − q 2 10 4 , r 01 = q 01 − q 10 q 11 2 , r 11 = q 11 , P 6 (X 5 , Y 5 ) is a smooth function and , from which we are able to get the universal fold of the system (1.2) ), where ξ 00 = r 00 r 4  11 , ξ 01 = r 01 r 11 , P 7 (X i6 , Y i6 ) is a smooth function and Case II: If p 20 − 2p 02 p 10 + p 00 p 2 02 < 0, then we make the following variable transformation: and τ is still recorded as s, then system (3.7)becomes , The system (10') becomes where r 00 = q 00 − q 2 10 4 , r 01 = q 01 − q 10 q 11 2 , r 11 = q 11 , P 6 (X i5 , Y i5 ) is a smooth function and 0, and we make another variable transformation as follows: from which one can get the universal fold of the system (1.2) where ξ 00 = −r 00 r 4  11 , ξ 01 = −r 01 r 11 , P 7 (X i6 , Y i6 ) is smooth function and

Numerical simulation
To further discuss the stability of the solution of the model equation, considering the actual situation of HABs occurrence and the critical influence of human behavior on the food chain model, taking a and c as control variables, the stability of the equilibrium state of the model equation was studied by using the bifurcation theory.We consider the following system (4.1) To get the model dynamics in the neighborhood of the Transcritical, Saddle-node, Hopf, Bogdanov-Takens bifurcations, we give a series of simulations to illustrate the emergence of complex patterns that arise for the system (4.1).We choose some parameter values for which (4.1) attains these bifurcations (see Figures 1a,2a,4a,7a).To observe the peculiarities of patterns that arise in the neighborhood of these bifurcations, we choose two sets of parameter values (see Figures 1b,2b) and four sets of network connection probability values (see Figures 5 and 8).

Transcritical bifurcation
When transcritical bifurcation occurs, it can be seen from the comparison that the equilibrium point's stability has changed in Figure 1 (left).In general, the distribution of plankton is uneven and will vary due to regional differences, as shown in Figure 1a (right).When a < a T c , the equilibrium point E 1 of (1.2) is in a stable state in Figure 1b (left), and the nodes of the network decrease first and then increase successively in Figure 1b (right).When a > a T c , the equilibrium point E 1 of (1.2) becomes unstable in Figure 1c (left), and the nodes of the network tend to another stable equilibrium point in different ways in Figure 1c (right).
Regarding biological mechanism, the physical meaning is the idealized initial level for the equilibrium point E 1 ; that is, the population density of plankton is zero.The half-saturation constant a can be used as the critical concentration to maintain the practical form of nutrients in the water for the average growth of algal cells.It can also be used to compare the ability of different phytoplankton to absorb nutrients.When light intensity, water temperature and other conditions were suitable but nutrient content was low (a < a T c ), phytoplankton with a smaller a value were likelier to become the dominant species.The growth of phytoplankton with a higher a value was limited due to the need for more nutrients.However, when nutrients are too abundant (a > a T c ), the community structure of phytoplankton will change significantly, which may lead to the rapid reproduction of some harmful phytoplankton; that is, transcritical bifurcation occurs (transition from one equilibrium point to another equilibrium point).It's also how organisms choose profitable ways to adapt to their environment.

Saddle-node bifurcation
For the equilibrium point E * , the cell density of plankton is not zero; E * > 0 is necessary for HABs to occur.As can be seen from Figure 2 (left), c = −1.95< 0 means de-eutrophication of effluent wastewater when the parameter crosses the critical value a * (other parameters were fixed) and the equilibrium points number goes from zero to one.Then, it splits into two, and their stability changes accordingly.With a = a * , there is one stable positive equilibrium in Figure 2a (left).Due to the network's action, each node region's stable point is different.That is, turing instability occurs in Figure 2a (right).As a < a * , there is one stable positive equilibrium point and one unstable positive equilibrium point in Figure 2b (left).Due to the effect of the network, almost all nodes simultaneously tend to be stable in an exponential decay way in Figure 2b (right).With a > a * , the positive equilibrium point disappears in Figure 2c (left).In this case, plankton tends to be extinct in this environment, resulting in too abundant nutrients.However, nutrients will not grow indefinitely.Still, it will reach a stable state, and under the action of the network, this stabilization is synchronous in Figure 2c (right).Similarly, the simulation showed that the network connection probability still could not affect the pattern.

Hopf bifurcation
In the following, a is selected as Hopf bifurcation parameters.Using the MATCONT toolbox for analysis, it is observed that the system has a Branch point (labelled as BP), Limit point (fold) bifurcation (labelled as LP) and two Hopf bifurcations (labelled as H), as shown in Figure 3  First Lyapunov coefficient = −1.289975× 10 −4 < 0, so the Hopf bifurcation is supercritical.When the bifurcation value a = 5.5015, the stability of the system changes abruptly, a limit cycle "pops out" around it and the system oscillates violently and continuously, as in Figure 4a.It can be seen from Figure 4 that if HABs break out at this time, the explosive growth of HABs must be related to the Hopf bifurcation behavior of the system.The periodic outbreak period at each node is different in the Nutrient-Plankton network due to the different environments in each region.However, the system at each node is periodic in Figure 4b for the whole network.With the increase of the network connection probability(p = 0.01, 0.05), it is evident that the outbreak period of each node gradually tends to be consistent in Figure 5.  Due to the limitation of environmental factors and physiological characteristics of algae, algae growth shows periodicity.When the surrounding environment is suitable, a specific population density of algae simultaneously begins to multiply and consume nutrients from the surrounding environment.The first HABs occur when the thickness of algae is above a particular value.When algae concentrations first peak (when the effects of HABs are most severe), they reach their limits.Due to nutrient restriction (Figure 4, nutrient concentrations decreased to low values during the same period), algal densities declined sharply until they reached deficient levels.Nutrient concentrations gradually recover due to the dramatic decrease in algae density, the degradation of algae metabolism and dead algae and the diffusion effect of ocean currents.Under these conditions, algae enters a period of rapid growth, forming the second HABs (the appearance of the second peak).

Bogdanov-Takens bifurcation
In this section, (a, c) are selected as Bogdanov-Takens bifurcation parameters.Using the MATCONT toolbox for analysis, we can get that (1.2) has a Bogdanov-Takens bifurcation point (labelled as BT) and a Cusp bifurcation point (labelled as CP), as shown in Figure 6 According to the previous analysis on the existence of the internal equilibrium point (Theorem 2), we know that when ∆ = 0, the two positive equilibrium points of (1.2) merge into one equilibrium point.The only positive equilibrium point E * is a cusp of codimension two, and the system state is shown in Figure 7.The positive equilibrium point E * is unstable in Figure 7a.Moreover, when the connection probability is small, as can be seen from the figure Figure 7b, the orbits near the equilibrium point far away from the equilibrium point are almost out of sync, which means eutrophication of effluent wastewater at each node in the network almost does not affect each other.However, with the increase of the connection probability (p = 0.035, 0.04), this lack of mutual influence has almost disappeared in Figure 8.Ultimately, we conclude the relationship between harmful algal blooms and bifurcation.First, marine eutrophication is the essential condition leading to the outbreak of HABs; when the nutrient concentration is low, there will be no HABs (nutrient concentration is below the threshold).Second, the explosion of HABs also depends on the point of algae concentration.On one hand, the location of bifurcation changes with the change of parameters.If the bifurcation location reaches the threshold of HABs occurrence, it may trigger HABs.On the other hand, it is often determined that the equilibrium point of the system loses stability and bifurcations occur (such as transcritical bifurcation, saddle-node bifurcation, Hopf bifurcation, B-T bifurcation), which causes a significant change in the concentration of algae.Once the threshold of the concentration of algae required for the occurrence of HABs is reached at some time, it will lead to an outbreak of HABs.Therefore, using the nonlinear dynamics theory to discuss the stability and bifurcation behavior of the model near the equilibrium point helps us study the HABs generation mechanism and explore the HABs outbreak's internal mechanism and the dependence on dynamic parameters and environmental factors.The study of the stability and bifurcation behavior of the plankton ecosystem plays an essential role in the prediction and early warning of HABs.

Conclusions
The paper deals with equilibria (local and global), Saddle-node, Transcritical, Hopf-Andronov and B-T bifurcation of the Nutrient-Plankton system on networks.In the numerical simulations, by controlling the change of parameters and observing the transformation process of the pattern, we observed several types of dynamics and the evolution process of the population of each node, which will bring great reference value to scientific research.When the system produces transcritical or saddle-node bifurcation, the change of parameters a and c will lead to changes in the number of equilibrium points or stability of (1.2), and the corresponding pattern will also change greatly.Moreover, in the simulation process, it was found that the change in network connection probability will not affect the changing trend of the pattern.The pattern does not change much before, after or on the critical value of the transcritical bifurcation (saddle-node bifurcation).The simulation results also show that when Hopf or B-T bifurcation was generated in (1.2), the network connection probability greatly influenced the pattern.
Regarding biological mechanism, the half-saturation constant a and de-eutrophication of effluent wastewater c < 0 directly affect the occurrence of HABs.When transcritical bifurcation (saddle-node bifurcation, B-T bifurcation) occurs, the system tends to be stable.The nutrient was relatively deficient and a limiting environmental factor at this time.Therefore, the growth rate of plankton is slow, the mortality rate and environmental loss rate are more significant than the growth rate and the concentration of plankton decreases or even tends to extinction.When Hopf bifurcation occurs, the