Qualitative features of a NPZ-Model

— Qualitative study of higher order non linear dynamical systems is a rewarding experi-ence and a great challenge. This reﬂective paper is an attempt to deeply analyze interaction features between nutrients, phytoplanktons and zooplanktons by building a so-called NPZ-Model . We used classical methods (of Lyapunov, Hopf, etc.) to examine existence, positivity, boundedness and stability of solutions. Our main contribution is the implementation of a meaningful space parameter that simultaneously guarantees instability of equilibria at the border and stability of the internal equilibrium. In the case of internal equilibrium instability, we observed the emergence of limit cycle which means the existence of periodical solutions.


I. INTRODUCTION
Experiments and sporadic observations of front zones (intense upwelling) show significant primary production (i.e., phytoplankton and zoo-plankton) in these zones. This overproduction due to nutrients influx from ocean floor to the Ekman layer (surface layer) promotes the growth of animal populations living in this layer, sardines for example.
Dynamics of living species in aquatic environments involve biological phenomena such as predation or interspecies competition. There are also physical phenomena such as flowing water, temperatures changes and range of salinity. These phenomena are closely linked and are subject to important multidisciplinary scientific investigationsthat is, oceanography, ecology, mathematics and physics. Several deterministic models governing species dynamics in the food chain have emerged recently. In particular the interaction between nutrients, phytoplankton and zooplankton ( [1], [2], [3], [4], [5], [6]) .
These studies can be roughly broken down into two categories. The first is dynamical systems where observables are densities (or number) of living species. In continuous case, the temporal evolution is described by an ODE 1 . The sec-ond category is that of reaction-diffusion systems where evolution of observables is spatio-temporal, one uses PDEs 2 in the continuous case. Sometimes one can add transport effect to take into account water velocity. Note that algebraic systems are handled in discrete cases.
Among these systems there is a relevant model governing dynamics of three trophic chain species mentioned above. It was initiated by Steele and Henderson in 1981 ( [7], [8] ) and then fleshed out by Baretta and Ebenhöh in 1997 [9]). Although they have given a perfect description of numerical modeling and simulations,their studies failed to grasp the description of space parameters and properties of solutions. We suppose that this is due to a plethora of parameters and the strong nonlinearity of equations. This is why we derive a model which will be called NPZ-Model, It's a 3-D dynamical system. Our study is first articulated on the space of admissible parameters. It's a set of conditions guaranteeing existence, positivity and boundedness of solutions. Then we rigorously analyzed the equilibrium stability by applying Lyapunov method and LaSalle principle of invariance. In the case of instability, we examined bifurcation in which a supercritical Hopf bifurcation has been discovered, this gives a stable limit cycle.
The model describes interaction of three species in food chain: Nutrients n(t), Phytoplankton p(t) and Zooplanktons z(t). Molecular concentration can be affected over time in an aquatic environment (ocean, lake, etc.). We have the following 2 System of Partial Differential Equation 3 All chlorophyllian aquatic organisms from plankton, some microscopic, others large 4 as nitrate, phosphate or silicate 5 All animal species that are part of plankton system of equations: The first equation of system (1) expresses nutrients evolution. The combined and vertical mixture results from an optimized and recycled nutrients supply for bacteria and phytoplanktons. Vertical mixing transports nutrients from the deep layer of water to the mixing layer (Ekman layer), this transport is expressed by the flux where the function S designates design the front zone intensity. It depends essentially on horizontal position (x, y) and water velocity v from the upwelling phenomenon. In this paper, we consider S(x, y, v) constant. This simplification eliminates the spatial aspect and does not take into account water velocity. Nutrients are consumed by phytoplankton with a characteristic saturation described by a Holling type II functional response 6 − βnp kn+n . Their recycling by bacteria is modeled by the last three terms in brackets of the first equation in system (1). A part of all organic waste and exudation of zoo-plankton are recycled by bacteria. However, bacteria dynamics is not included in the model.
Phytoplankton's dynamics is considered in the second equation of system (1). Their concentration depends on nutrients which constitute a food source βnp k n + n . It is also reduced by natural mortality (−µ p p) and by zoo-plankton's grazing (predation) − αηp 2 α + ηp 2 z which is Holling type III functional response. 6 Which takes into account nutrient consumption time. Therefore, catch rate decreases with increasing nutrient density. Grazing is also considered as zooplankton growth term with the factor γ αηp 2 α+ηp 2 z − µ z z 2 (third equation of (1) . This factor is there to take into account that only a part of nutrients is converted into zoo-plankton biomass. The Other part (1 − γ) is recycled. Mortality of zoo-plankton is assumed to be quadratic.This hypothesis implies the existence of a super predator acting on zoo-plankton whose dynamics are not explicitly considered [12]. It may also be an interspecies competition. Parameters (Table I) are described as in the work of Pasquero & al [13]. Their numerical values will be discussed in the paper through the space of admissible parameters (III.1 ).
As the Holling type-II functional response is one of the most realistic forms to represent predation rate of predators on their preys [14], we change the Holling type-III functional to a Holling type-II like αp 1+αhp .
We then assume that there is neither implicit predation on the zoo-plankton nor interspecies competition. We assumed that the steady constancy of nutrient recruitment mid saturation k n is greater than the constant concentration of nutrients under the mixing layer n 0 ( i.e., k n > n 0 ). Then system (1) becomes what we will call NPZ-Model: To reduce occurrences of α in the NPZ-Model, we denote η ≡ α.h. Then we set n = µ n − 1, Let's define or, equivalently, where Φ(t) = (n(t), p(t), z(t)) T , Denoting by F i,x the partial derivative of F i with respect to x, we have: Thus, Jacobian matrix of field F is written: Thus, we have the main notations that we will use for qualitative study of NPZ-Model. Note that intermediate ratings will be defined as needed.

A. Solutions Properties
Proposition III.1. (Existence, uniqueness, positivity & boundedness) System (4) associated with initial condition has a unique positive and bounded solution.
Proof: (Proposition III.1) Elements of Jacobian matrix of field F (8 ) are all continuous, then the function F is locally Lipschitzian. Thus by Cauchy-Lipschitz theorem [15], system (4) has a unique local solution.
To show the positively invariance of R 3 + , let's consider these elements of boundaries: Thus, by choosing a function L(n, p, z) = −n, we see that F, ∇L ≤ 0. From the foregoing, we can deduce by the barrier theorem ( theorem 16.9 in [16] ) that the vector field point into the domain along Γ 1 and by analogy along Γ 2 and Γ 3 . Hence, R 3 + is a positively invariant domain, i.e., for any positive initial condition, solutions remain positive.
To prove boundedness, let's define with a > 0 such that Thuṡ We have Thuṡ Let us set By applying Gronwall's inequality, we obtain: Accordingly, for a large t , Θ(t) < e t0 Θ(t 0 ) + aSn 0 e −t0 , and the solutions are bounded.
Proof: (Proposition III.2) To calculate equilibrium points, we replace the first equation by the sum of three others and we obtain the following system : (10) From the third equation of system (10), we have z = 0 where p = µz αγ−ηµz . By treating the case where z = 0 , we have two equilibrium points :  with relations: and the fact that A > 0 because R 2 < 1, we have : Hence We have also Therefore Ultimately Definition III.1 (Space of admissible parameters).
We define the space of admissible parameters P ad by Remark III.1. The space of admissible parameters defined above (III.1) P ad is not empty . Table  II gives examples of admissible parameters.
Proof: (Theorem III.1) We will directly show the global stability before checking condition (i) .
As ( * n, * p, * z) is an equilibrium, then from the first equation of system ( 4) we have : By substituting this equality intoV 1 , we get: Let's consider : Thus, .
We have : .
We have also So, we can rewriteV 1 like : Let's look aṫ From the second equation of system (4) we have: By substituting µ p inV 2 we obtain: .
Thus, V satisfy conditions of Lyapunov stability theorem, [15, page 194], then e 3 is asymptotically stable in G. Now let's show condition (i): We thus obtain If we choose ω 0 such that We have also Let's set the condition And let δ 1 and δ 2 the two real roots of the polynomial Aδ 2 + Bδ + C such that δ 1 < δ 2 . Then, there exist admissible parameters (ω 0 = 1.79 and ∆ = 0.0077, See Fig 1 for the rest of parameters ) verifying the condition (24). So it's sufficient to choose δ < δ 1 or δ > δ 2 and we will have Aδ 2 + Bδ + C > 0, which implies : As we have obviously ζ( * n, * p, * z) = 0 ≤ 0 , we conclude that e 3 ∈ G.

C. Bifurcation around equilibrium e 3
We are interested in system behavior according to upwelling intensity, the parameter S. We set a range of admissible parameters where equilibrium e 3 is stable. We denote S 0 the parameter resulting from this choice and we will call it bifurcation parameter. By varying S in [S min , S max ] = [0.1, 5.85], we look for the bifurcation value denoted by S b , this will be the point where equilibrium e 3 loses its stability. Let P (X) be the characteristic polynomial of system linearized around By fixing S 0 = 0.85, we observe numerically that H 2 (3.38) = H 3 (3.38) = 0 (Fig 2). Then S b = 3.38 is a bifurcation value. It is a Hopf bifurcation according to the Liu criterion [18] (III.1).
Lemma III.1 (Liu criterion). Suppose that the following conditions are satisfied Then there is a simple Hopf bifurcation.
Looking at the typology of eigenvalues (Fig 2) 7 , we have : (i) For S ∈ [S min , S b [ , a negative real eigenvalue, and two conjugate complexes of strictly negative real part. The equilibrium e 3 is therefore locally stable. (ii) For S ∈]S b , S max ], we have a negative real eigenvalue, and two conjugate complexes with a strictly positive real part. The equilibrium e 3 is therefore unstable. We conclude from the Poincaré-Andronov-Hopf theorem [19], that there exists a limit cycle for S ∈]S b , S max ] =]3.38; 5.85].

IV. RESULTS & DISCUSSION
After the definition of NPZ-Model, we are interested in qualitative study of system. We made sure of the existence, uniqueness, positivity and boundedness of solutions. Then, we calculated the stationary points and established sufficient conditions for their positivity. We have discovered that the two equilibria on the border of R 3 + are unconditionally unstable, which is interpreted as a correlation of dependence between the three species. We have shown thanks to Lyapunov's theorem [15] , that the internal equilibrium is stable in a subset G ∈ R 3 + . This set is not unique and depends strongly on parameters of system.
We have written a python program (NPZ-Prog.py) 8 whose main task is to select a range of admissible parameters ( i.e., checking proposition III.1 ). Then, we used these parameters to evaluate respectively equilibria and the associated eigenvalues.Thus, by drawing 10 5 ranges of admissible parameters, it was observed that *  (Fig 4).
Then comes the integration of the NPZ system with the odein function of python which is based on a Runge-Kutta method of order 4 with adaptive steps.
As shown in Theorem III.1, there is a range of parameters where the equilibrium e 3 is stable. In this case, simulations confirm that eigenvalues are of negative real part and two of them are com- 8 Attached with additional files plex conjugate. Orbits spiral and converge towards equilibrium e 3 ( Figure 6). By varying upwelling intensity S, we have shown that a super-critic Hopf bifurcation occurs, that is to say, there is a value of S where the equilibrium e 3 loses its stability and a limit cycle appears. We have not rigorously studied stability of this cycle but the numerical simulations conjecture its stability. This phenomenon of periodic maintained oscillations was observed by Edwards & Bees [10] where they considered zooplankton predation of Holling type-III. This shows that our choice of Holling type-II does not remove this fundamental aspect of the original model (Fig 5).
In the case where we have these periodic solutions, we are interested in primary production which is defined as the growth term in phytoplankton's dynamics and defined by The quantity P P is periodic because n and p are, looking at evolution of its average with respect to S ∈ [3.4, 5.85] (after bifurcation ), we see that the average primary production is constant ( P P =  0.29 ) with a growing standard deviation through S (Figure 7).

V. CONCLUSION
To study the dynamics of nutrients, phytoplanktons and zoo-planktons, we focused on a system of ODE (4). There is a space of admissible parameters where solutions are well defined. We also have an interior equilibrium point. In the case where the equilibrium is unstable, there is a stable limit cycle which means periodical solutions. We had algebraically determined four equilibria, of which e 4 / ∈ R 3 + . This is due to the sufficient conditions stated on positivity of three others (9). It would be interesting to look for positivity conditions for the 3. four equilibria and to study if necessary, the bistability of e 3 and e 4 . A biological interpretation of the parameter and time spaces would strengthen understanding of the model. Finally, it would be enriching to consider transport-diffusion effect and to couple the NPZ-model with an upwelling equation.

APPENDIX
Tables III, IV and V show a list of some functions used to perform this paper. There are three python codes 9 : NPZ _ Sympy_ Calcul.py , NPZ_ Prog.py, Hurwitz_ Prog.py . Tables below give functions and their outputs