The Role of Photosynthesis Processes in Dynamics of Plant Community

1. Introduction

Phytoplankton is a key link in the global carbon cycle. The study of the mechanisms of primary production formation is paramount. Primary production is the result of the interaction of carbon dioxide and minerals when exposed to solar energy. Methods to calculate primary production using the assimilation number are known [1,2,3] This scheme is simplified because it does not take into account the variability of photosynthetic substances in phytoplankton cells. Chlorophyll and other plant cell photosynthetic substances play an important role in the primary production processes [4,5,6]. In what follows, we will refer to chlorophyll «a» as the main photosynthetic substance in cells. It actively absorbs solar energy. The absorbed photons are used to start the flow of electrons across the thylakoid membrane. The energy thus converted is used in the chemical reaction for the synthesis of organic [7,8]. Chlorophyll grows together with biomass by metabolism. The ratio of chlorophyll to phytoplankton is usually called the chlorophyll quota [8]. It is measurable but quite variable [9,10]. According to experts, the proportion can vary from 0.001 to 0.097 [11]. The ratio of chlorophyll to biomass varies depending on the state of cells and environmental conditions [12,13,14,15]. Such variations significantly affect the productivity of phytoplankton.

A substance of ATP (adenosine triphosphate) is a universal energy source for all organic biochemical processes. There is a traditional way to assess the functional indicators of primary production through the concentration of chlorophyll per unit of biomass [16]. An alternative method based on the amount of ATP is considered in a number of works [17,18,19,20]. Such works became relevant again after the methods for measuring ATP levels in cells were invented and improved [21]. In the study [19] a statistical modeling method for estimating the rate of biomass growth is presented. The method is built both on the basis of measurements of the proportion of chlorophyll per unit of biomass (chlorophyll quota), and on the basis of the estimated proportion of ATP.

We developed mathematical models that are based on the chlorophyll quota and on the concentration of energy-intensive substances. Just like many researchers [22,23,24,25,26], we based our constructions on the Droop model [27,28]. Traditionally, the functioning of phytoplankton in the aquatic environment is described by the Monod dependence [29,30]. In same times, microorganism growth continues for some time after the nutrient resources are depleted. This phenomenon is important for photosynthesis process. A detailed analysis of these mechanisms is discussed in [23]. Monod's kinetics are unable to imitate this phenomenon. Droop's concept makes it possible to separate the nutrient uptake rate from the growth rate. Many authors use it for a detailed description of nutrientically limited reproduction of phytoplankton [31,32,33]. Our models describe an important impact of energy-intensive substance concentration and chlorophyll quote on biomass dynamics.

2. Droop’s model

The Droop model for describing the growth of phytoplankton in a chemostat [34] is chosen as the base one. This model is the system of differential equations. For an open system, it has the form:

$$\dot{p}=\left(\mu \left(q\right)-D\right)p,$$

(2.1a)

$$\dot{s}=D\left(s_{in}-s\right)-\nu \left(s, q\right)p,$$

(2.1b)

$$\dot{q}=\nu \left(s, q\right)-\mu \left(q\right)q.$$

(2.1c)

The model has next denotations:

• s is the concentration of the nutrient substance;
• p is the phytoplankton biomass in volume unity;
• ν is the pace of nutrient assimilation;
• q is the cell quota;
• μ is the rate of biomass growth;
• D is the rate of water flow;
• s_in is the concentration of the nutrient substance in the inner flow.

Cell quota is the concentration of a mineral substance in a phytoplankton cell.

The pace of nutrient assimilation is described as a function

$$\nu \left(s,q\right)={\nu }_0\frac{s}{s+s_0}\left(1-\frac{q}{Q_1}\right).$$

(2.2)

The rate of biomass growth is described as a function

$$\mu \left(q\right)={\mu }_0\left(1-\frac{Q_0}{q}\right).$$

(2.3)

A function μ(q) describes the velocity of phytoplankton growth, and ${\mu }_0$ is the maximum of this velocity. Parameter $Q_0$ and $Q_1$ are the minimum and maximum of cell quota q respectively. At $q<Q_0$ biomass growth stops. This form for the function μ(q) is the traditional form for the Droop model [35]. We propose that the next conditions are true: all constant parameters have positive values, and $Q_0<Q_1$ .

The properties of Droop model were studied very well [28,31,33]. The equilibria are found in the system of equations:

$$\left(\mu \left(q\right)-D\right)p=0, D\left(s_{in}-s\right)-\nu \left(s, q\right)p=0, \nu \left(s, q\right)-\mu \left(q\right)q=0.$$

The single trivial equilibrium $\left(p_e=0, {s_e=s}_{in}, q_e\right)$ exists always. The single positive equilibrium $(p_{*}>0, s_{*}, q_{*})$ exist at next conditions:

$${D\le \mu }_0\left(1-\frac{Q_0}{Q_1}\right) , D{q}_{*}\le \nu \left(s_{in},q_{*}\right) \mathrm{with} q_{*}={\mu }^{-1}\left(D\right).$$

The Droop model has an attractor. A common content of the nutrients will be denoted as $x\left(t\right)=p\left(t\right)·q\left(t\right)+s\left(t\right)$. From model (2.1) follows a property $\underset{t\to \infty }{\mathit{lim}}x\left(t\right)=s_{in}$. The proof follows from equation $\dot{x}=D{s}_{in}-Dx$ as a corollary of the system (2.1). We denote this attractor as set

$$A=\left\{\left(p,s, q\right)\in R_{+}^3: p·q+s=s_{in}, q\in \left[Q_0,Q_1\right] \right\}.$$

(2.4)

Two possibility equilibria are contained in the set A.

3. Droop’s model modification with chlorophyll quota (Model C)

Next our models are extensivity of Droop model. O. Bernard proposed a modification of the Droop model [23]. The effect of lighting on phytoplankton growth is described using a restricted function: $\mu \left(q\right)={\mu }_0\left(1-\frac{Q_0}{q}\right)\cdot \chi \left(\cdot \right)$.

The function $\chi \left(\cdot \right)$ is discussion. In [8] a transition was made from the rate of biomass growth to the rate of carbon uptake per unit of chlorophyll. Thus, the dynamic process is transferred from the mass balance to the category of photosynthetic transformations. At this stage, the authors introduce the concept of chlorophyll quota as the ratio of the total concentration of chlorophyll in the system to the number of cells containing chlorophyll. Many researchers have a mean that the chlorophyll quota is somewhat dependent on the cell quota [36]. Our modification of the Droop’s model uses this proposition.

The Model C modification has three equations (2.1) from the Droop model, to which the chlorophyll quota dynamics relation is added:

$$\dot{c}=c_m\gamma \left(q,c\right)-\rho c.$$

(3.1)

Chlorophyll quota ${c\in [0,c}_m]$. The normalized chlorophyll quota seeks to normalized cell quota. Chlorophyll consumption by photosynthesis process [37] is taken into account too.

The system of equations (2.1)+(3.1) is Model C.

The formula (2.2) is true here. The growth rate of phytoplankton biomass has the form instead (2.3):

$$\mu \left(q,c,I,T\right)={\mu }_0\left(1-\frac{Q_0}{q}\right)·(\alpha \phi \left(I\right)·\frac{c}{c_m}+{\alpha }_0)·\psi \left(T\right).$$

(3.2)

The parameters $\alpha$, ${\alpha }_0$ are associated with the rate of biomass growth during the light and dark stages of photosynthesis, respectively.

In addition to the mineral nutrition constantly entering the system, PAR and water temperature are among the influencing factors. The growth rate depends on PAR as function $\phi \left(I\right)=\frac{I}{I+I_0}$ and from water temperature as function $\psi \left(T\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left(\frac{T-T_{opt}}{T_1-T_0}\right)}^2\right)$.

The influence of cell quota to chlorophyll quota is regulated of the function

$$\gamma \left(q,c\right)={\gamma }_0\left(\frac{q}{Q_1}-\frac{c}{c_m}\right).$$

(3.3)

3.1. Existence and stability of equilibrium at the Model C for constant environmental

We propose that the parameters $I, T$ of the environment are constants in this subsection. We lose these parameters from function in this subsection.

We have for define of possible equilibrium a next system:

$$\left(\mu \left(q,c\right)-D\right)p=0,$$

(3.4a)

$$D\left(s_{in}-s\right)-v\left(s, q\right)p=0,$$

(3.4b)

$$v\left(s, q\right)-\mu \left(q,c\right)q=0,$$

(3.4c)

$$c_m\gamma \left(q,c\right)-\rho c=0.$$

(3.4d)

Theorem 3.1. The single trivial equilibrium $(p_e=0,s_e=s_{in}, q_e, c_e)$ exists always. The single positive equilibrium $(p_{*}>0, s_{*}, q_{*}, c_{*})$ exists under inequations:

$$\mu \left(Q_1,c\left(Q_1\right)\right)-D\ge 0,$$

$$\nu \left(s_{in}, q_{*}\right)-\mu \left(q_{*},c_{*}\right)q_{*}\ge 0.$$

The values $q_{*}, c_{*}$ is defined unambiguously in system of equations:

$$\mu \left(q,c\right)-D=0,$$

$$c_m\gamma \left(q,c\right)-\rho c=0.$$

Proof.

1) Trivial equilibrium. The linear function $c\left(q\right)$ is defined at the equation (3.4d) of system (3.4). The unique value $q_e$ is defined from equation (3.4c) of system (3.4) for $s_e=s_{in}$. The left part of this equation has properties of strongly decreases about q and different signs in interval $q\in [Q_0,Q_1]$.

2) The positive equilibrium $(p_{*}>0, s_{*}, q_{*}, c_{*})$ is defined at equations system (3.4) for $p>0$. The linear function $c\left(q\right)$ is defined at the equation (3.4d) of system (3.4). It is increasing function. The unique $q_{*}$ is defined at the equation (3.4a) if left part of this equation has nonnegative value for $q=Q_1$

$$\mu \left(Q_1,c\left(Q_1\right)\right)-D\ge 0.$$

The components $q_{*}, c_{*}=c\left(q_{*}\right)$ must do an inequation

$$\nu \left(s_{in}, q_{*}\right)-\mu \left(q_{*},c_{*}\right)q_{*}\ge 0.$$

We have then a unique $s_{*}\in [0, s_{in}]$ from equation (3.4c). The component $p_{*}$ is calculated in the equation (3.4b). To end we have equilibrium $(p_{*}, s_{*}, q_{*}, c_{*})$. □

Two possibility equilibria belong to set A (2.4). There are the one trivial equilibrium $p=0$ and no more than one positive equilibrium $p>0$ for positive constant parameters in the Model C. It is possible to study the stability of these solutions in the range of changes in key parameters using a linearized system by the Lyapunov method [38,39]. The illustration of these dynamical properties is on Figure 1 in projection on plane $(\alpha ,Q_0)$.

The area of instability of trivial equilibrium and asymptotical stability positive equilibrium has blue color. The «green» area of parameter variation corresponds to the existence of both equilibria. Of these, the trivial equilibrium is stable, and the positive one is unstable. A stable trivial equilibrium exists in «red» set.

Depending on the area in which the model parameters fall, the phase portrait of the main components $(p, s, q)$ of the model system changes. The phase portraits of the system (2.1)+(3.1) show on Figure 2 for the points indicated by numbers 1 - 3 on Figure 1.

The solution of the Model C is calculated from different points of the area of the initial approximation. The system parameters are set based on hitting each of the designated stability areas. The parameter $\alpha$ is constant for all solutions. The parameter $Q_0$ of minimal сell quota is variable. In all cases, the solution asymptotically tends to the attractor (Section 2). For the "blue" region, solutions from different initial states tend to a positive equilibrium. The "attraction" of the equilibrium point has a global impact on the dynamics of the solution. In the "green" region, solutions from different initial states tend to different points belonging to the attractor, and then "wander" around the attractor. And in the "red" region, the "attraction" of the trivial equilibrium is global in nature. Thus, the minimal сell quota has a big role for quality characteristics of solutions.

4. Model with separation of light and dark stages of photosynthesis (Model E)

A light stage of photosynthesis includes an accumulation of energy-intensive compounds. These substances are then involved in the synthesis of organic matter at the dark stage of photosynthesis. The model has the three equations of Droop model (2.1) and the next equations:

$$\dot{c}=c_m\gamma \left(q,c\right)-\delta (q,I)·c,$$

(4.1a)

$$\dot{r}=\epsilon \delta (q,I)·c-\mu \left(q,r,T\right)r.$$

(4.1b)

The parameter $r\ge r_0$ has a sense of growth pool. This is an amount of the energy-intensive substances (ATP) in biomass. The parameter $r_0$denotes a value of the energy pool below which the fusion reaction does not occur. We have the coefficient of growth biomass in the form:

$$\delta \left(q,I\right)={\delta }_0·\left(1-\frac{Q_0}{q}\right)·\phi \left(I\right),$$

(4.2a)

$$\nu \left(s,q\right)={\nu }_0\frac{s}{s+s_0}\left(1-\frac{q}{Q_1}\right),$$

(4.2b)

$$\gamma \left(q,c\right)={\gamma }_0\left(\frac{q}{Q_1}-\frac{c}{c_m}\right),$$

(4.2c)

$$\mu \left(q,r,T\right)={\mu }_0·\left(1-\frac{Q_0}{q}\right)·\beta (1-\frac{r_0}{r})·\psi \left(T\right).$$

(4.2d)

The function $\beta (1-\frac{r_0}{r})$ describes an influence of $r$ to biomass growth similar to the influence of the cell quota.

The system of equations (2.1)+(4.1) is Model E.

4.1. Existence and stability of equilibrium at the Model E for constant environmental

We propose that the parameters $I, T$ of the environment are constants in this subsection like in subsection 3.1.

To define of possible equilibrium we have the next system:

$$\left(\mu \left(q,r,T\right)-D\right)p=0,$$

(4.3a)

$$D\left(s_{in}-s\right)-\nu \left(s, q\right)p=0,$$

(4.3b)

$$\nu \left(s, q\right)-\mu \left(q,r,T\right)q=0,$$

(4.3c)

$$c_m\gamma \left(q,c\right)-\delta \left(q,I\right)c=0,$$

(4.3d)

$$\epsilon \delta \left(q,I\right)·c-\mu \left(q,r,T\right)r=0.$$

(4.3e)

Theorem 4.1. The Model E has no more than one trivial equilibrium$(p_e=0,s_e=s_{in}, q_e, c_e,r_e)$. This model has no more than one positive equilibrium$\left(p_{*}>0, s_{*}, q_{*}, c_{*},r_{*}\right)$.

Proof.

1) Parameters $\left(q_e, c_e,r_e\right)$ of the trivial equilibrium is defined from equations (4.3c), (4.3d), (4.3e) of system (4.3)

$$\nu \left(s_{in}, q\right)-\mu \left(q,r,T\right)q=0,$$

(4.4a)

$$c_m\gamma \left(q,c\right)-\frac{1}{\epsilon q}\nu \left(s_{in}, q\right)r=0,$$

(4.4b)

$$\epsilon \delta \left(q,I\right)·c-\mu \left(q,r,T\right)r=0.$$

(4.4c)

The strong decrease function $r=r\left(q\right)$ is denoted at equation (4.4a) of this system. Further there is a strong decrease function $c=c\left(q\right)$ at equation (4.4c) with help of function $r=r\left(q\right)$.These functions are put in equation (4.4b).

$$c_m\gamma \left(q,c\left(q\right)\right)-\frac{1}{\epsilon q}\nu \left(s_{in}, q\right)r\left(q\right)=0.$$

The function on the left side is strictly decreasing. This equation has no more than one solution $q_e$. A no more than one $c_e,r_e$ is defined from this fact.

2) The positive equilibrium $(p_{*}>0, s_{*}, q_{*}, c_{*},r_{*})$ is defined at equations system

$$\mu \left(q,r,T\right)-D=0,$$

(4.5a)

$$\left(s_{in}-s\right)-qp=0,$$

(4.5b)

$$v\left(s, q\right)-Dq=0$$

(4.5c)

$$c_1\gamma \left(q,c\right)-\frac{1}{\epsilon }Dr=0,$$

(4.5d)

$$\epsilon \delta \left(q,I\right)·c-Dr=0$$

(4.5e)

This system is received from the system (4.3). The scheme of equilibrium definition is next. The strong increase function $q=q\left(s\right)$ is denoted at equation (4.5c) of this system. Further the strong decrease functions $r=r\left(s\right)$, $p=p\left(s\right)$, $c=c\left(s\right)$ are denoted at the equations (4.5a), (4.5b), (4.5e). These functions are put to equation (4.5d). The function in left part is a strong increase function in the variable s. Therefore, the equation (4.5d) has no more than one solution. And at most one positive equilibrium is defined. □

Two possibility equilibria belong to set A (2.4). One of them is trivial equilibrium $p=0$ and other one positive equilibrium $p>0$.

In addition to the minimum cell quota, numerical experiments have shown a significant effect on the stability of the solutions of the parameter β. The domains of existence and stability of both solutions are in the plane of the ratios of these parameters (Figure 3).

Solutions are defined of a trivial or of a positive equilibrium. The tendency to attraction by some equilibrium is the more obvious, the farther the matching parameters point is located from the region of instability. The phase portraits of the main components of the Model E illustrate different behaviors of the solution trajectories depending on the area from which the model parameters are selected (Figure 4).

The dynamics of the solutions in Model C and Model E have similar properties. The Droop’s submodel makes some general basis in these models. A zone of instability is smaller for Model E. It represents a very narrow boundary region in which the Jacobian eigenvalues vanish. The complex dynamics can be in this “green” set on plane $(\beta ,Q_0)$ parameters. The attractor largely affects the dynamics of solutions. These properties connect with existence and stability of the equilibrium solutions. In general, the solutions dynamics are subject to attraction by the attractor of the Droop's submodel.

Let us illustrate the behavior of the solution of both models with consistent coefficients near the attractor.

The Figure 5 shows the dynamics of biomass along the attractor under constant external conditions calculated using the Model C (left) and the Model E (right). As can be seen, according to Model C, the quasi-stationary state of biomass is reached almost four times longer. The dynamics of other indicators on the attractor is shown in Figure 6.

Figure 6 show that Model E makes it possible to reach the equilibrium state of the system in a shorter period. Obviously, the reason is in the structure of the equations of the Model E. The chlorophyll quota directly affects the growth of biomass. The Model E configuration allows cellular structures to quickly "adjust" to the biomass changes that occurred in the previous time step. Sequential (one after another) changes in the main components of the system allow you to approach equilibrium faster.

5. Model E: the optimal solution on the attractor

Studying the behavior of solutions on the attractor A (2.4) gives a more complete picture of the properties of the Model E. Let us consider the trajectories under the additional condition of optimality according to the criterion of the largest increase in phytoplankton biomass p for a fixed time.

First, we transform the Droop model for solutions on the attractor A.

Lemma 5.1.The Droop model on the attractor A is transformed to the form:

$$\dot{p}=(\mu \left(q,r,T\right)-D)p,$$

(5.1a)

$$\dot{q}=\nu \left(p, q\right)-\mu \left(q,r,T\right)q,$$

(5.1b)

where

$$\tilde{\nu }\left(p, q\right)=\nu \left(s_{in}-pq, q\right).$$

(5.2)

If the initial condition for solving the system of equations (5.1) is positive and admissible, then the solution satisfies the condition $p\left(t\right)q\left(t\right)\in (0, s_{in}]$ for any t.

Proof. According to the definition of the attractor A (2.4), the formula $pq+s=s_{in}$ is valid. From here we substitute s into $\nu \left(s, q\right)$ in the equation (2.1c) of the Droop model (2.1). Then the equation (2.1b) can be excluded from this system. Equations (5.1) are obtained.

Due to the validity of formula (5.2), the solution of system (5.1) at any time t satisfies the condition ${\partial }_t\left(pq\right)=\left(\nu \left(p, q\right)-Dq\right)p$ and $\tilde{\nu }\left(p, q\right)=0$ for $s_{in}-pq=0$. This implies the second part of the Lemma assertion. □

5.1. Formulation of the optimal control problem

Assume the dynamics of phytoplankton biomass is modeled on the interval $[0, t_e]$. We will proceed from the condition of phytoplankton biomass maximization at the final time $p\left(t_e\right)\to max$. This is equivalent to the condition ${\int }_0^{t_e}\dot{p}\left(t\right)dt={\int }_0^{t_e}\left(\mu \left(q,r,T\right)-D\right)pdt\to max$. Let us take into account equations (4.1) for c, r in the form of the comparative smallness of $\dot{r}$ in the equation (4.1e). Then we get the following optimization criterion:

$${\int }_0^{t_e}\{\left(\mu \left(q,r,T\right)-D\right)p-\omega {\left(\epsilon \delta \right(q,I)·c-\mu \left(q,r,T\right)r)}^2\}dt\underset{c\in \left[0,c_m\right], r\ge r_0}{\to }max,$$

(5.3)

where ω is a relatively small positive number.

The relevant restrictions are the equations (5.1) under given initial conditions for p and q.

Problem (5.1)+(5.3) is a Lagrange problem [40,41] and is solved on the basis of the Pontryagin maximum principle [42] according to the Lagrange approach.

5.2. Solution of the optimal control problem

Compose the Lagrange function:

$$L\left(p,q,c,r,\lambda \right)={(\lambda }_1-1)\left(\mu \left(q,r,T\right)-D\right)p+\omega {\left(\epsilon \delta \right(q,I)·c-\mu \left(q,r,T\right)r)}^2+{\lambda }_2[\tilde{\nu }\left(p, q\right)-\mu \left(q,r,T\right)q]$$

Then the Euler-Lagrange equation is:

$${\dot{\lambda }}_1=\left({\lambda }_1-1\right)\left(\mu \left(q,r,T\right)-D\right)+{\lambda }_2{\partial }_p\tilde{\nu }\left(p, q\right)$$

$${\dot{\lambda }}_2=\left({\lambda }_1-1\right){\partial }_q\mu \left(q,r,T\right)p+2\omega \left(\epsilon \delta \left(q,I\right)c-\mu \left(q,r,T\right)r\right)\left(\epsilon {\partial }_q\delta \left(q,I\right)c-{\partial }_q\mu \left(q,r,T\right)r\right)+{\lambda }_2[{\partial }_q\tilde{\nu }\left(p, q\right)-{\partial }_q(\mu \left(q,r,T\right)q\left)\right]$$

under conditions

$${\lambda }_1\left(t_e\right)=0, {\lambda }_2\left(t_e\right)=0.$$

The functions c, r are determined from the condition of minimizing L with respect to these variables.

5.3. The optimal control problem modification

The minimizing c is calculated from the equation $\epsilon \delta \left(q,I\right)·c-\mu \left(q,r,T\right)r=0$ for $q>Q_0$ and $\frac{\mu \left(q,r,T\right)r}{\epsilon \delta \left(q,I\right)}\le c_m$. Denote by $r_m$ the value r at which the last inequality becomes an equality.

The minimum condition for L with respect to r turns into a condition for minimizing with respect to r the expression:

$$L\left(p,q,c,r,\lambda \right)={(\lambda }_1-1)\left(\mu \left(q,r,T\right)-D\right)p+{\lambda }_2[\tilde{\nu }\left(p, q\right)-\mu \left(q,r,T\right)q]$$

It boils down to minimizing over r the expression:

$$[{(\lambda }_1-1)p-{\lambda }_{2}q]\mu \left(q,r,T\right)$$

Denote $\xi =[{(\lambda }_1-1)p-{\lambda }_{2}q]$. Then the following condition is satisfied:

$$r=\left\{\begin{array}{c}{r}_m, \xi <0 \\ r_0, \xi \ge 0\end{array}\right..$$

Then the equations for ${ \lambda }_{\mathrm{1,2}}$ are simplified:

$${\dot{\lambda }}_1=\left({\lambda }_1-1\right)\left(\mu \left(q,r,T\right)-D\right)+{\lambda }_2{\partial }_p\tilde{\nu }\left(p, q\right)$$

$${\dot{\lambda }}_2=\left({\lambda }_1-1\right){\partial }_q\mu \left(q,r,T\right)p+{\lambda }_2[{\partial }_q\tilde{\nu }\left(p, q\right)-{\partial }_q(\mu \left(q,r,T\right)q\left)\right].$$

Finally, problem (5.1)+(5.3) acquires the following form.

Theorem 5.1. To solve problem (5.1)+(5.3) it is necessary to solve the following boundary value problem

$$\dot{p}=(\mu \left(q,r,T\right)-D)p,$$

(5.4a)

$$\dot{q}=\tilde{\nu }\left(p, q\right)-\mu \left(q,r,T\right)q,$$

(5.4b)

$${\dot{\lambda }}_1=\left({\lambda }_1-1\right)\left(\mu \left(q,r,T\right)-D\right)+{\lambda }_2{\partial }_p\tilde{\nu }\left(p, q\right),$$

(5.4c)

$${\dot{\lambda }}_2=\left({\lambda }_1-1\right){\partial }_q\mu \left(q,r,T\right)p+{\lambda }_2[{\partial }_q\tilde{\nu }\left(p, q\right)-{\partial }_q(\mu \left(q,r,T\right)q\left)\right]$$

(5.4d)

for a given function $r_{\left(t\right)}=\left\{\begin{array}{c}{r}_m, \xi \left(t\right)<0 \\ r_0, \xi \left(t\right)\ge 0\end{array}\right.$ with.

The boundary conditions are:

The functions included in the equations correspond to the Model Е (formula (4.2)) with a change in (5.2).

Proof. The proof is given before the statement of the theorem.

Test calculations for problem (5.4) were carried out under the following conditions: the left boundary conditions were chosen so that the initial point lay on the attractor. The right boundary conditions for ${ \lambda }_{\mathrm{1,2}}$ are set equal to zero. The values of other parameters are presented in Table. 1 (Section 6). The calculations were made under constant external conditions. The ϕ, ψ values are constant and correspond to the average daily values of light irradiance and temperature. Experimental examples showed that in this case the solution to the boundary value problem can be found using standard numerical procedures. Some initial values may require some regularization.

More details about the selected values of environmental indicators and their changes during the daily and annual cycle are described in Section 6, devoted to numerical experiments.

6. Calculation experiments

The comparative analysis of solutions’ properties of Models C, E is a basic content of computer experiments. A data of research is defined in Table 1 from references [43,44,45].

In calculations, qualitative properties of dynamics of the main model variables are analyzed depending on significant parameters. These are α, β and Q₀. All parameters are involved in the function of specific biomass growth rate. The parameter $\alpha$ is "responsible" for the influence of illumination in the first model, and the parameter $\beta$ characterizes the influence of ATP in the second model. Preliminary numerical experiments showed the greatest significance of the parameters α и β (apart from Q₀) for the dynamics of the described processes.

6.1. Dynamics in a variable environment. Influence of Light and Darkness

The behavior of models on an attractor under constant external conditions is considered (Section 4). It is interesting to know how sensitive they are to changes in external conditions. We simulated the conditions similar to a laboratory experiment [20]. At some initial concentrations, the dynamics of the system was calculated sequentially over a time interval corresponding to a day and for 100 days. The illumination function was modeled as follows: for 12 hours the light level was constant (60 Ein⋅m^-3day^-1). This was followed by a twelve-hour period of complete darkness. In this case, the temperature of the medium was assumed to be constant (7° C).

The numerical experiment demonstrates the behavior of the model components corresponding to the description of the results of a laboratory experiment under similar conditions. The daily cycle is divided into equal periods of presence/absence of light. During the light period (Figure 7, left), a decrease in the proportion of chlorophyll is observed. This is in complete agreement with the observations described in [20]. In the biological sense, such process is associated with an increase in biomass and the dispersion of the proportion of chlorophyll in it. The model behavior of the chlorophyll quota is due to the fact that it tends to correspond to the intracellular quota of nutrients. The cellular pool of nutrients is consumed during the progressive growth of biomass. The value of the chlorophyll quota decreases following the cell quota to its minimum value, below which the biomass growth stops. Next, the consumption of chlorophyll for the formation of ATP stops. In a 100-day period, after a certain surge, there is a weak trend towards a decrease in the chlorophyll quota in combination with daily fluctuations (Figure 7, right).

The proportion of ATP changes in antiphase with respect to the chlorophyll quota. This is also in line with the findings presented in [20]. Biologically, this means that during the light period, the energy to start the enzymatic reaction comes from the capture of photons by the reaction centers of the chlorophyll molecules. Energy-intensive substances, therefore, acquire the ability to accumulate. During the dark period, energy for chemical transformations can only be obtained by releasing it as a result of the breakdown of ATP. Accordingly, their pool is reduced (Figure 8, left).

In the long-term light/dark alternation, the ATP pool tends to decrease (Figure 8, right). With a long-term progressive increase in biomass (Figure 10, left), the trend towards a decrease in both indicators seems to be adequate (Figure 7, 8, right).

Other numerical experiments have shown that phytoplankton acquires a tendency to extinction with a significant deterioration in external conditions, namely, a change in temperature relative to the optimum, a reduction in the period of illumination and a decrease in the level of illumination. The dynamics of plankton with a threefold reduction in the above indicators is shown in Figure 10 (right).

At the same time, the chlorophyll quota and the proportion of ATP reveal an interesting trend (Figure 9). At first, the indicators go in the anti-phase. Under sufficient initial concentrations of phytoplankton and nutrients, at the first stage, an increase in biomass is observed. This is followed by an increase in the chlorophyll quota. A rapid consumption of ATP is due to the low level of light energy coming from outside. Then the trend changes. Both indicators show weak growth. But a slight accumulation of ATP and chlorophyll do not cope with the lack of energy and unfavorable temperature. Phytoplankton biomass is decreasing (Figure 10, right).

Figure 9. Dynamics of the chlorophyll quota and ATP under periodic lighting changes (one day on the left, 100 days on the right). Daily period of light is 4 hours.

Figure 9. Dynamics of the chlorophyll quota and ATP under periodic lighting changes (one day on the left, 100 days on the right). Daily period of light is 4 hours.

Figure 10. Dynamics of biomass (g/m³) under periodic lighting changes during 100 days (daily period of light 12 hours on the left, daily period of light 4 hours on the right).

Figure 10. Dynamics of biomass (g/m³) under periodic lighting changes during 100 days (daily period of light 12 hours on the left, daily period of light 4 hours on the right).

6.2. Annual cycle

6.2.1. Solution of the initial problem on the annual period

Under changing environmental conditions, the model systems of equations cease to be autonomous. The annual dynamics calculated for both models undergo two characteristic bursts, spring and autumn.

The annual cycle in Figure 11 shows the dynamics of phytoplankton biomass, typical for the Far Eastern Seas of Russia. The simulation of the annual cycle is provided by seasonal changes in ambient temperature and illumination. The daily change in illumination is considered.

It is generally accepted that the spring rise in biomass is higher than the autumn one [46,47]. In the Far East region, the temperature regime has some peculiarities. In the Far Eastern seas, cold-loving species (mainly diatoms) [48] predominate, and a certain proportion is occupied by more thermophilic (dinophyte) species [49]. The most favorable temperature for the entire plant complex forms in late summer-early autumn up to its middle. The spring increase in biomass is associated with the consumption of nutrients stored during the winter. The numerical experiment simulates the natural conditions of the region. The external illumination function is a parabolic curve with a maximum in mid-July. Temperature is a similar smooth function with a maximum in the second half of August and relatively constant values throughout the year except for summer. In accordance with the winter accumulation of mineral nutrition, their initial concentrations are high compared to average observed concentrations.

The temperature and light functions are based on real data provided by the Center for Multifunctional Satellite Environmental Monitoring, by the Institute of Automation and Control Processes, the Far Eastern Branch, Russian Academy of Sciences (IACP FEB RAS), Vladivostok [50]. Based on remote sensing data, initial approximations for the numerical solution of model systems are formed. The species composition of phytoplankton in the Far Eastern seas was formed under the influence of regional natural conditions and is tolerant to low temperatures. Biomass growth is quite intensive throughout the year. This is the reason for the intensive consumption of nutrients. Both models support this. The concentration of nutrients throughout the entire vegetative cycle approaches its value in the input flow (Figure 11, Part 1).

The value of the cell quota fluctuates between the minimum and maximum values, and these fluctuations occur in antiphase to the volume of biomass, especially during periods of its rise. With intensive growth of phytoplankton, the intracellular pool of nutrients is dispersed over an increasing number of cells. Therefore, the value of the cell quota decreases. Daily fluctuations in the chlorophyll quota (Figure 11, Part 2) become more pronounced. As already mentioned [51], this is consistent with experimental observations.

The dynamics of the share of ATP in phytoplankton repeats the annual dynamics of the chlorophyll quota with some time lag. Fluctuations around daily averages are much less noticeable in comparison with the variability of the chlorophyll quota. This is also confirmed by the conclusions contained in [20].

6.2.2. Solution of the problem of optimal control in the annual period.

Annual fluctuations in biomass can also be considered as a solution to the optimal control problem (see Section 5). By modeling the external conditions in the manner described above, it is possible to obtain the biomass dynamics presented in Figure 12.

The obtained solution of the boundary value problem (5.4) quite closely repeats the solution of the initial Cauchy problem for models C and E with a characteristic spring and autumn burst of biomass. Moreover, the maximum values of seasonal concentrations are either comparable in absolute values (spring outbreak) or equal (autumn outbreak).

Thus, under the condition that the initial-boundary conditions belong to the attractor of system (4.3), the numerical solution of the optimal control problem (5.4) can be used to predict phytoplankton dynamics. In this case, the dimension of the space of the system components decreases. At the same time, one should not ignore the fact that the numerical solution of boundary value problems is associated with problems of convergence of computational procedures. This requires additional regularization algorithms, which cannot always be effectively applied.

6.3. Climate fluctuations

The illumination and temperature functions correspond to the climatic features of the Far Eastern seas. In this case, both models demonstrate adequate behavior of the phytoplankton biomass. The solution of Model E is characterized by smoother dynamics compared to Model C. This is due to the structure of the model equations in Model E. The coefficient ρ at the chlorophyll quota ceases to be constant and becomes a function $\delta \left(q,I\right)$ depending on external parameters. Thus, a large proportion of the dispersion due to high gradients of the illumination function is accounted for by chlorophyll. It is partially repaid due to the balance correlation with the share of ATP. This provides more smoothness for the biomass growth function. Thus, Model E seems to be more preferable for understanding the trend of biomass changes. When modeling anomalous scenarios, the emerging complex dynamics of phytoplankton is explained by changes in external parameters, and not by the properties of the model itself. On the basis of Model E, a number of calculations were undertaken to simulate climatic fluctuations during a given season. Figure 13 shows the results of changes in phytoplankton dynamics, corresponding to a 2° C cooling in winter (left) and a 2° C warming in summer. The general conclusion is that the deviation of average temperatures from the optimum within the tolerance interval has a negative effect on the growth rate. But, in general, such temperature fluctuations do not have a significant effect on the average values of biomass during the year. Thus, Model E adequately reflects the regulatory role of the temperature factor in the phytoplankton growth process [52]. This fact is a qualitative conclusion of most experimental observations [53,54,55].

More sensitive Model E is in relation to the parameter β. Even with a slight increase in its value (above 10^-3), that is, with an increase in the amount of ATP consumed for the synthesis of organic matter, the specific biomass increases markedly (Figure 14). At the same time, the temperature still has an insignificant effect on the annual dynamics.

6.4. Long-term dynamics

Figure 15 shows changes in biomass dynamics over many years. Scenarios corresponding to different values of β in the presence/absence of climate change are considered.

The results of model calculations are as follows: if the value β is small and the replenishment of the nutrient resource is limited, then there is a gradual extinction of the phytoplankton biomass. Numerical experiments have shown that the smaller the value of the parameter β, the higher the negative gradient of the growth function. With a slight increase in the proportion of ATP that stimulates the synthesis reaction, the downtrend gradually breaks and turns into an exponential one. With a larger increase in β, the dynamics becomes stable, including a short-term surge, and then a transition to a quasi-stationary regime. In this case, the step in β affects the maximum biomass value at the burst stage and the average quasi-stationary value, but does not affect the dynamic contour itself (Figure 15). Among the endogenous parameters of the model, the value of the minimum cell quota Q₀ has a comparable impact. The influence of this parameter does not need to be considered separately, since it has a clear biological meaning. The larger the minimum cellular quota of nutrients required to start growth, the slower it is and vice versa. The cumulative effect of the parameters β and Q₀ on the equilibrium values of the specific biomass is shown in the Figure 16.

7. Conclusions

The presented models are based on Droop's concept of the relation [28] between growth rate and intracellular nutrient quota. This model makes it possible to describe the growth of biomass after the cessation of the supply of mineral nutrition to the system.

In addition to the Droop model, both proposed models are based on the concept of chlorophyll quota [56,57]. Experimental studies [37] give grounds to suspect a relationship between the proportion of chlorophyll in phytoplankton and its cell quota. Simplified, it can be a linear relationship between indicators with a certain proportionality coefficient. The presented models allow to avoid simplification and describe the relationship between chlorophyll and cell quota, taking into account the mechanics of primary production processes. Furthermore, Model E describes the stage of accumulation of the internal energy reserve necessary to start the enzymatic reaction. One of the model systems includes the percentage of ATP in biomass. The form of Droop dependence was borrowed to describe the growth of biomass in the dark stage of photosynthesis (Model E), when illumination as a necessary resource dries up. Then the only source of energy needed to start the enzymatic reaction is the energy-intensive substances stored in the light stage of photosynthesis [58,59]. Modern research methods make it possible to measure this quantity directly [17,18]. As a result of the enzymatic reaction, organic substances are formed. They are the building blocks of the living cell [60,61].

The Model C has complex dynamics when the stabile equilibria don’t exist. These dynamics are submitted to the attractor of the Droop model (see Lemma 3.2).

The behavior of the chlorophyll quota and the proportion of ATP in the biomass, calculated according to Model E, is in full agreement with the results of laboratory experiments presented [20]. In the absence of light, phytoplankton cells use ATP as an energy source. Chlorophyll begins to be consumed only when a sufficient amount of energy is accumulated. The reaction of organic matter synthesis is slower. Therefore, chlorophyll is also consumed slowly. Consequently, its highest concentration in phytoplankton cells is recorded in absolute darkness [62] (Figure 7, left). The ATP concentration, on the contrary, drops to a minimum (Figure 8, left), since the consumption of energy-intensive substances increases in the dark due to the absence of other energy sources [63].

Modeling of the annual cycle shows adequate dynamics in both models (Figure 13). Seasonal habitat changes lead to spring and autumn bursts of biomass productivity. This agrees with known field observations [46,47]. Given the availability of data from remote observations and the results of laboratory experiments, the presented models can be used to assess the biological productivity of aquatic ecosystems.

Model E can be transformed by solving the optimal control problem. The application of optimal management methods is relevant if the plant community is considered as a certain population striving for the maximum possible increase in the total biomass. The described and investigated approach (Section 5) allows us to reduce the space of the original components of the model. This can greatly simplify the collection of experimental data in the practical use of models.