Mathematical Modeling of the Evolutionary Dynamics of a Planktonic Community Using a Discrete-Time Model

Abstract

This study proposes a discrete-time eco-genetic model of a planktonic community that includes zooplankton and two competing phytoplankton haplotypes with and without a toxicity trait. The Holling type II response function describes predator consumption. We use the Ricker model to consider density limitation and regulation. The model is analytically and numerically studied. The loss of stability of fixed points occurs via the Neimark–Sacker scenario and a cascade of period-doubling bifurcations. The model reveals bistability and multistability. Therefore, the initial conditions can determine which of the coexisting dynamic modes will be attracted. If the competition of haplotypes is weaker than their self-regulation, then the variation in the current densities of community components can shift the observed dynamics, while the evolution direction remains unchanged. The ratio of haplotype fitnesses and predator pressure generally determines the asymptotic genetic composition of phytoplankton. If competition of haplotypes is higher than their self-regulation, then the bistability of monomorphic fixed points occurs when the displacement of one haplotype by another depends on initial conditions. The presence of predators can maintain the genetic polymorphism of the prey. This system shows dynamic modes similar to experimental dynamics: oscillation with delay, long-period antiphase fluctuations, and cryptic cycles emerging due to rapid evolution.

1. Introduction

Plankton are the basis of the marine food chain, supporting life from invertebrates and crustaceans to fish and large mammals. At the same time, phytoplankton provides more than half of the world’s oxygen needs and absorbs carbon dioxide, potentially mitigating one of the causes of global warming [1]. Understanding the dynamics of plankton communities is extremely important to predict fisheries [2,3,4,5] and the consequences of global climate change [6,7,8,9]. Using mathematical models in studies of phytoplankton dynamics is a desired tool for managing aquatic ecosystems and resources by better understanding the processes. Mathematical models allow virtual experiments and studies that may be impossible or undesirable in natural ecological systems. Direct estimates of phytoplankton biomass and tens of thousands of species of herbivorous zooplankton that make up the plankton community are traditionally complex and expensive [10]. Recently, biomass estimates of phytoplankton based on satellite observations have become increasingly popular [11,12,13,14]. The identification and interpretation of the part of the spectrum emitted by plankton based on satellite monitoring data remains a non-simple challenge [15,16,17,18]. In this context, theoretical modeling of plankton abundance dynamics remains relevant, as it allows us to understand the main mechanisms and physical processes that influence the dynamics of plankton and the consequences of its changes for the oceans and atmosphere.

Three-component models (NPZ models) that consider the interconnected dynamics of phytoplankton (P), zooplankton (Z), and critical nutrients (N) are the most popular in modeling plankton ecosystems. Analogous to plants, photosynthesis and nutrient availability determine the growth of phytoplankton. NPZ models, as a rule, consider zooplankton as an herbivorous species that feeds on phytoplankton [14,19]. Models of phytoplankton and zooplankton dynamics without explicit consideration of nutrients (the Rosenzweig–MacArthur models) are equally popular [20]. Zooplankton consists of herbivorous and predatory species; self-limitation allows one to explicitly or implicitly describe it during modeling [21].

Researchers pay considerable attention to modeling phytoplankton “blooms” characterized by a sharp growth of biomass, while toxic blooms negatively affect aquaculture, coastal tourism, and human health [22,23]. The effect of phytoplankton toxicity can vary. Some species produce toxins that poison organisms, inhibiting their growth or leading to death. Others dramatically increase biomass, leading to oxygen starvation and the death of marine inhabitants [24]. Delay equations are often applied to model the blooming effect [25,26]. Recurrent equations, widely used in modeling populations and communities [27] and allowing us to describe delay effects naturally, are virtually unused in this subject area.

In mathematical models of a predator–prey plankton community without separating phytoplankton into toxic and non-toxic species, the harmful properties of phytoplankton contribute to the zooplankton mortality function [25,26,28,29]. Zooplankton, which consume phytoplankton, can obtain toxins that can lead to their death. Modeling the contribution of consuming toxic phytoplankton to zooplankton development uses various functions. The difference between the contribution to birth and mortality rates is not infrequently expressed only in coefficients. An additional equation in the model to separate phytoplankton into toxic and non-toxic species allows us to consider that some phytoplankton species do not release toxins unless damaged [30]. At the same time, zooplankton can reduce the level of predation in response to phytoplankton-emitted toxins, thus avoiding the adverse effects of these chemicals. For example, Copepoda are selective in their food choices and avoid consuming phytoplankton that produce toxins [31].

Within the eco-evolutionary modeling approach that has recently gained popularity [32,33], the toxic properties of phytoplankton are its adaptive defense traits against predator grazing. Phytoplankton can evolve rapidly in response to zooplankton predation. Comparing experimental data on the development of the Brachionus calyciflorus rotifer community and the unicellular green alga Chlorella vulgaris in a chemostat and model trajectories confirmed the ability of algae to evolve rapidly by increasing the level of defense against grazing of rotifers and decreasing, as compensation, its birth rate [34]. The presence of a trade-off between the ability of algae to defend themselves from predators (or defense) and a decrease in their reproduction rate was usually a key factor in the evolution of the community. In particular, genetic analysis showed the presence of different algal genotypes in the chemostats, experimentally confirming the relationship between algal size, population density, and nutrition. However, a decrease in the reproduction rate did not always accompany the defense. In this context, the toxicity naturally transitions to the food value of algae, since this value of algae with maximum protection against rotifer grazing is zero. This fact is due to the refusal of zooplankton to consume toxic phytoplankton. The unicellular alga Chlamydomonas reinhardtii demonstrates another interesting example of defense. In the presence of rotifers, its protected genotypes form palmelloid cell clusters that are too large for the predator to ingest [35,36]. Therefore, the algal food value for zooplankton of the protected phytoplankton genotype is zero. This type of phytoplankton implicitly impacts the community, since competing for resources with non-toxic phytoplankton regulates zooplankton food availability.

Clonal selection models [32] naturally describe the asexual reproduction of algae (including C. vulgaris) and are the most adequate for modeling algae evolution. This is consistent with the results of genetic analysis that reveal multiple clones in laboratory populations [37]. In eco-evolutionary clonal selection models, the number of equations describing phytoplankton dynamics equals the number of algal genotypes. Without a detailed study of the eco-evolutionary NPZ model with clonal selection of prey and continuous time, J.F. Fussmann et al. [32] have shown its ability to capture the main features of experimental dynamics. The model reveals the long-period out-of-phase cycles of predator and prey during prey evolution in chemostats with many different algal genotypes. Classical predator–prey non-evolutionary dynamics also emerge [38], which are oscillations with a quarter-period lag of predator dynamics behind prey in chemostats with a single algal genotype. In addition to out-of-phase cycles, so-called cryptic cycles [39] were discovered in which the density of prey remains almost unchanged while the density of predators fluctuates. Such dynamics of the time series give the impression that there is no interspecific interaction between rotifers and algae. However, the rapidly changing genetic composition of prey species can support such cryptic cycles when the defense is not accompanied by a sharp decrease in the growth rate of algae. Thus, the masking effect of rapid evolution is that it can conceal trophic species interactions. Similar dynamics have previously been experimentally discovered in bacteriophage communities [40], which may indicate the presence of hidden cycles in various natural systems with rapid evolution.

In this work, we propose a model with a minimal number of discrete-time equations that considers the basic features of the evolutionary dynamics of a plankton community based on the classical ideology for this field [28,29,41]. Without consideration of the phytoplankton phenotypes located in the interval of the trade-off curve between the defense and growth rate of phytoplankton, we focus only on the extreme values of the trait, assuming that two haplotypes with zero food value (toxic) and unity food value (non-toxic) represent phytoplankton. The present study proposes a three-component discrete-time model of the planktonic community, where zooplankton feed on non-toxic phytoplankton and two phytoplankton haplotypes compete with each other for resources. Accordingly, in the community under consideration, “predator–prey” relationships are realized between zooplankton and the non-toxic phytoplankton haplotype, and competition between the two genetic groups of phytoplankton.

2. A Model of the Dynamics of the Plankton Community Considering the Genetic Composition of Phytoplankton

To describe the dynamics of toxic and non-toxic phytoplankton, we use a discrete-time Ricker’s competition model, which considers both interspecies and intraspecies competition [42,43]. Note that an argument in favor of applying such models to the analysis of phytoplankton community dynamics is the daily rhythm of such a community since many processes occurring in the phytoplankton community synchronize with circadian rhythms–cyclic fluctuations in the intensity of various biological processes associated with the alternation of day and night. Additionally, the Ricker model correctly describes changes in phytoplankton biomass because they occur through cell division. Additionally, in such a community, the model considers the consumption of phytoplankton by zooplankton. We propose toxicity as an adaptive trait of phytoplankton defense against zooplankton grazing, and such an approach has recently gained popularity [32,33]. On the basis of this, we represent phytoplankton by two genetically different groups (or haplotypes): toxic and non-toxic.

$$\{\begin{array}{l}{X}_{1,n+1}=A_1{X}_{1,n}\exp (-a_1{X}_{1,n}-b_1{X}_{2,n}\cdot (1-u))\\ X_{2,n+1}=A_2^{}\cdot (1-u)\cdot X_{2,n}\exp (-b_2{X}_{1,n}-a_2{X}_{2,n}\cdot (1-u))\end{array},$$

(1)

where X₁ and X₂ are the abundances of the toxic and non-toxic phytoplankton haplotypes, n is the day number, and A₁ and A₂ are the growth rates of the first and second phytoplankton haplotypes, respectively. Coefficients a₁ and a₂ are the self-limitation, and b₁ and b₂ are parameters characterizing the intensity of competitive interactions between haplotypes. u is the proportion of the abundance of non-toxic phytoplankton consumed by zooplankton, which depends on the biomass of prey and predator and will be formalized later.

Furthermore, let us denote X_n = X_1,n + X_2,n as the total abundance of phytoplankton, and then p_n = X_2,n/X_n is the proportion of the non-toxic phytoplankton haplotype, which corresponds to the allele frequency without the toxicity trait in the population. The concentration of the toxic haplotype (or frequency of the other allele) is denoted by q_n = X_1,n/X_n. Then, based on X_1,n = q_nX_n and X_2,n = p_nX_n, we obtain the dynamics of haplotype frequencies in neighboring generations from Equation (1):

$$\{\begin{array}{l}{q}_{n+1}=A_1\cdot q_n{X}_n\exp \left(-x_n\left(a_1{q}_n+b_1{p}_n\cdot (1-u)\right)\right)/X_{n+1}\\ p_{n+1}=A_2\cdot (1-u)\cdot p_n{X}_n\exp \left(-x_n\left(a_2{p}_n\cdot (1-u)+b_2{q}_n)\right)\right)/X_{n+1}\end{array}.$$

(2)

Now, based on the equality p + q = 1, the dynamics equation of phytoplankton abundance has the form:

$$X_{n+1}=X_n\left(A_1\cdot q_n\cdot e^{-X_n\left(a_1{q}_n+b_1{p}_n\cdot (1-u)\right)}+A_2\cdot (1-u)\cdot p_n\cdot e^{-X_n\left(a_2{p}_n\cdot (1-u)+b_2{q}_n)\right)}\right).$$

(3)

The average fitness of phytoplankton equals W_n = X_{n + 1}/X_n, then Equation (3) can be rewritten as $W_n=W_1\cdot q_n+W_2\cdot p_n$, where W₁ and W₂ are the fitnesses of the toxic and non-toxic haplotypes, respectively; as seen from Equation (3), they are determined by the formulas

$$W_1=A_1\cdot e^{-X_n\left(a_1{q}_n+b_1{p}_n\cdot (1-u)\right)}, W_2=A_2\cdot (1-u)\cdot e^{-X_n\left(a_2{p}_n\cdot (1-u)+b_2{q}_n)\right)}.$$

Given all the above, we obtain the following system:

$$\{\begin{array}{l}{p}_{n+1}=W_2\cdot p_n/W_n\\ X_{n+1}=X_n\cdot W_n\end{array}, \mathrm{where} q_n=1-p_n$$

We also use the Ricker model to describe the dynamics of zooplankton abundance:

$$Z_{n+1}=R\cdot Z_n\exp (-\omega \cdot Z_n),$$

(4)

where Z is the abundance of zooplankton, ω is its self-limitation coefficient, R = rα(X_2,n) is the growth rate of the predator, and r is its maximum possible value. To describe the feeding process of zooplankton α(X_2,n) = α(X_np_n), we use the Holling type II trophic function with predator saturation: α(X_np_n) = p_nX_n/(X* + p_nX_n), where X* is the abundance of non-toxic phytoplankton at which the reproductive potential of the predator will be equal to half of the maximum possible.

Predator feeding changes the abundance and genetic composition of prey:

$$p_n{X}_n\to p_n{X}_n-{\alpha }_0\cdot z_n\cdot \alpha (p_n{X}_n)=p_n{X}_n\left(1-\frac{{\alpha }_0\cdot z_n\cdot \alpha (p_n{X}_n)}{p_n{X}_n}\right)=p_n{X}_n(1-u),$$

where $u=\frac{{\alpha }_0\cdot z_n\cdot \alpha (p_n{X}_n)}{p_n{X}_n}=\frac{{\alpha }_0\cdot Z_n}{X*+p_n{X}_n}$ is the proportion of non-toxic phytoplankton consumed by zooplankton, and α₀ is the biomass conversion coefficient (α₀ >> 1). Thus, predators influence prey reproduction, competition, and survival. As a result, the community model takes the following form:

$$\{\begin{array}{l}{p}_{n+1}=W_2\cdot p_n/W_n\\ X_{n+1}=X_n\cdot W_n\\ Z_{n+1}=\frac{r\cdot X_n}{X*/p_n+X_n}{Z}_n\exp (-\omega \cdot Z_n)\end{array},$$

(5)

where $W_n=W_1\cdot q_n+W_2\cdot p_n$, $W_1=A_1\cdot e^{-X_n\left(a_1{q}_n+b_1{p}_n\cdot (1-\frac{{\alpha }_0\cdot Z_n}{X*+p_n{X}_n})\right)}$, $W_2=A_2\cdot \left(1-\frac{{\alpha }_0\cdot Z_n}{X*+p_n{X}_n}\right)\cdot e^{-X_n\left(a_2{p}_n\cdot \left(1-\frac{{\alpha }_0\cdot Z_n}{X*+p_n{X}_n}\right)+b_2{q}_n)\right)}$, q_n = 1 − p_n.

The following substitutions of variables a₂X→x, ωZ→z and coefficients ρ = b₁/a₂, x* = a₂X*, a₂α₀/ω = α, φ = b₂/a₂, a = a₁/a₂ transform model (5) to a simpler form

$$\{\begin{array}{l}{p}_{n+1}=w_2\cdot p_n/w_n\\ x_{n+1}=x_n\cdot w_n\\ z_{n+1}=\frac{r\cdot x_n}{x*/p_n+x_n}{z}_n\exp (-z_n)\end{array},$$

(6)

where $w_1=A_1\cdot e^{-x_n\left(a{q}_n+\rho p_n\cdot \left(1-\frac{\alpha \cdot z_n}{x*+p_n{x}_n}\right)\right)}$, $w_2=A_2\cdot \left(1-\frac{\alpha \cdot z_n}{x*+p_n{x}_n}\right)\cdot e^{-x_n\left(p_n\cdot \left(1-\frac{\alpha \cdot z_n}{x*+p_n{x}_n}\right)+\varphi q_n)\right)}$, $w_n=w_1\cdot q_n+w_2\cdot p_n$, q_n = 1 − p_n.

Here, x and z describe the relative abundances or densities of phytoplankton and zooplankton, respectively; α is the average number of prey (in relative units) consumed by a relative unit of predator; and x* characterizes the half-saturation constant of the predator. In this case, the predator’s half-saturation constant equals x*/p_n, and is a dynamic variable since it depends on the density of the non-toxic phytoplankton haplotype. The coefficients ρ, a, and φ characterize the competitive relationships between toxic and non-toxic phytoplankton and self-limiting processes. In general, the coefficients ρ and φ are coupling coefficients; their equality to zero reduces model (1) to a system of two noninteracting groups of individuals between which there is no competition. In the case of a = 1, the contribution to interspecies competition for limiting resources from each haplotype is the same.

Applying model (6) to describe plankton community dynamics in the case of high predation pressure can lead to negative values of phytoplankton density. To ensure the nonnegativity and boundedness of solutions, we must consider additional conditions for the parameters and initial conditions of the model. In particular, if the expression $\frac{\alpha \cdot z_n}{x*+p_n{x}_n}$ exceeds unity, then, due to overconsumption of prey by the predator, the non-toxic phytoplankton haplotype and zooplankton become extinct. Note that the growth rate and competitiveness of toxic haplotypes determine the success of their development.

3. Fixed Points of Model (6) and Scenarios of Community Development

One of the scenarios for community development whose dynamics can be described by model (6) is extinction, and it corresponds to the fixed point with coordinates:

$$\overline{x}=0 \mathrm{and} \overline{z}=0,$$

(7)

i.e., the stationary densities of the species equal zero, while the genetic composition can be either monomorphic $\overline{p}=0$, $\overline{p}=1$, or polymorphic $0<\overline{p}<1$ at $A_1=A_2$. Note that this scenario is possible at $A_{1,2}<1$ and follows from the nature of Equation (1). The Ricker model demonstrates degeneration when the reproductive potential is less than 1. Figure 1 shows the model (6) trajectories that demonstrate the extinction of the community with different genetic compositions of the phytoplankton population. Figure 1a shows that the community goes extinct at $\overline{p}=1$, when the monomorphic phytoplankton population is represented only by the non-toxic haplotype. Here, initially, the prey population is polymorphic, that is, both genetic groups are present. Then the more adapted haplotype, which is primarily determined by its growth rate $A_1<A_2$, displaces the less adapted competitor and then, due to the influence of factors of different natures, goes extinct. At $A_1>A_2$, other monomorphic solution $\overline{p}=0$, $\overline{x}=0$, $\overline{z}=0$ becomes attractive, where $\overline{p}=0$ corresponds to the community with the toxic phytoplankton haplotype (Figure 1b). With $A_1=A_2$ and $A_{1,2}<1$, Figure 1c shows that the phytoplankton population is polymorphic: variable p stabilizes at some intermediate level different from unity. Obviously, this solution is unstable, and the condition $A_1=A_2$ separates the attraction domains of solutions $\overline{p}=0$, $\overline{x}=0$, $\overline{z}=0$, and $\overline{p}=1$, $\overline{x}=0$, $\overline{z}=0$.

Another scenario of the development of such a community is the extinction of zooplankton $\overline{z}=0$, while phytoplankton exists. Here, at least one of the parameters A₁ and A₂ is greater than 1. The zooplankton extinction can be caused both by its low growth rate r in combination with limitation processes and by insufficient food (or non-toxic phytoplankton haplotype density). In this regard, in each iteration step of model (6), it is necessary to check the condition $1-\frac{\alpha \cdot z_n}{x*+p_n{x}_n}>0$. If this is not true, then there has been a complete removal of the non-toxic phytoplankton haplotype, which, in a closed system, will lead to zooplankton extinction. In general, the community without zooplankton can be represented by a monomorphic or polymorphic phytoplankton population.

The coordinates of the fixed point corresponding to the toxic haplotype monomorphism are as follows:

$$\overline{p}=0, \overline{x}={\overline{x}}_1=\frac{\ln A_1}{a}, \overline{z}=0.$$

(8)

The fixed point corresponding to the non-toxic phytoplankton haplotype monomorphism is:

$$\overline{p}=1, \overline{x}={\overline{x}}_2=\ln A_2, \overline{z}=0.$$

(9)

The solution corresponding to polymorphism is:

$$\overline{p}=\frac{a\ln A_2-\varphi \ln A_1}{\ln A_2(a-\rho )+\ln A_1(1-\varphi )}, \overline{x}={\overline{x}}_1+{\overline{x}}_2=\frac{\ln A_1-\rho \ln A_2}{a-\rho \varphi }+\frac{a\ln A_2-\varphi \ln A_1}{a-\rho \varphi }, \overline{z}=0.$$

(10)

Fixed points (8) and (9) are monomorphic solutions when only one of the phytoplankton haplotypes survives. Solution (10) is polymorphic and corresponds to the coexistence of two competing genetic groups of phytoplankton in the absence of zooplankton.

The results of model (1) can be extrapolated to model (6). In particular, in the works [44,45,46,47] for model (1) it has been shown that if the competition between species is higher than their self-limitation, i.e., a₁a₂ − b₁b₂ < 0, then the stable coexistence of two competing populations is impossible. Naturally, this statement is also true for model (6), i.e., if φρ > a, then the more adapted haplotype displaces the less adapted one, or there are parametric regions where the current ratio of genetic group densities determines the outcome of evolution, that is, the so-called “bistability trap” arises, when the more adapted haplotype can be displaced by the less adapted one.

Another group of community development scenarios characterizes zooplankton’s presence in the ecosystem. Here, the following cases are possible.

• The reduced community consists of a predator and a monomorphic prey population, i.e., when the toxic phytoplankton haplotype is not competitive and goes extinct. In this case, the nontrivial monomorphic fixed point of model (6) is the solution of the following system:

  $$\{\begin{array}{l}\overline{p}=1\\ A_2\cdot \left(1-\frac{\alpha \cdot \overline{z}}{x*+\overline{x}}\right)\cdot e^{-\overline{x}\left(1-\frac{\alpha \cdot \overline{z}}{x*+\overline{x}}\right)}\\ \overline{z}=\ln \frac{r\cdot \overline{x}}{x*+\overline{x}}\end{array}=1;$$

  (11)

2.

A complete community consists of a predator and a polymorphic prey population, where stationary densities are the solution of the system of equations:

$$\{\begin{array}{l}{w}_1=1\\ w_2=1\\ \frac{r\cdot \overline{x}}{x*/\overline{p}+\overline{x}}\exp (-\overline{z})=1\end{array}.$$

(12)

Note that systems (11) and (12) are nonlinear and, depending on parameter values, have from one to three roots corresponding to the number of nontrivial solutions of model (6) (Figure 2). The stability domains of fixed points are separated from each other by lines of transcritical bifurcations, at which stability exchange occurs between trivial and semitrivial, semitrivial and nontrivial solutions. System (11) must have three roots to provide the existence of three fixed points in system (12) for the complete community.

The graphical solution of Equation (11) for the reduced community demonstrates that as the value of the parameter characterizing the growth rate of the non-toxic phytoplankton haplotype increases, the number of roots changes (Figure 2a–d). Initially, there is one solution, which we denote as 1₁ since the curves intersect at a single point. Furthermore, with an increase in the values of the coefficient A₂, a saddle-node bifurcation occurs, and the unstable point 1₂ appears; it corresponds to the point of tangency (Figure 2b), which with further growth of the bifurcation parameter splits into stable point 1₃ and unstable 1₂; thus, the system has three nontrivial solutions. Then, the first fixed point 1₁ merges with unstable 1₂, and the only attractive solution 1₃ remains, born due to the saddle-node bifurcation (Figure 2d). The presence of a saddle-node bifurcation in the reduced system describing the dynamics of the “non-toxic phytoplankton haplotype—zooplankton” community provides the emergence of this bifurcation in the system corresponding to the complete community; it depends on the values of the parameters x*, r, α, and A₂ (Figure 2e). For the other coefficients, one can always select values at which tangent bifurcation will occur. Note that this bifurcation can emerge, but the attractive point will be outside the region of feasible values. For a complete community, Figure 2f shows the occurrence of tangent bifurcation and the variable p values that are greater than 1, and as seen, point 1₃ is outside the domain of feasible values. On the other hand, if we compare Figure 2d,g, we can see that in the first case, solution 1₃ remains after the collapse of the merging of points 1₁ and 1₂, while solution 1₁ in Figure 2g is attractive, and a variation in one of the parameters, for example, decreasing A₁ or φ, will lead to the occurrence of a saddle-node bifurcation. Thus, for some values of the coefficients, the system can have two fixed points, one of which is attractive—the point of intersection of two curves, and the other is the point of their tangency (Figure 2b). Moreover, the curves corresponding to the equations of system (12) can intersect at two points, and this is associated with a change in the curve shape due to variation in the values of parameters characterizing the limitation and competition in the phytoplankton population. In particular, at small a, the system can have two nontrivial fixed points, one attractive and the other not (Figure 2h).

4. Stability Regions of the Fixed Points of Model (6): Evolutionary Scenarios of Community Development

To study the stability of model (6) solutions, we write its Jacobian:

$$J=\Vert \begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\Vert =\Vert \begin{array}{ccc}F{p^{\prime }}_p& F{p^{\prime }}_x& F{p^{\prime }}_z\\ F{x^{\prime }}_p& F{x^{\prime }}_x& F{x^{\prime }}_z\\ F{z^{\prime }}_p& F{z^{\prime }}_x& F{z^{\prime }}_z\end{array}\Vert ,$$

(13)

where $Fp=\frac{w_2\cdot p}{w_1\cdot (1-p)+w_2\cdot p}$, $Fx=x(w_1\cdot (1-p)+w_2\cdot p)$, and $Fz=\frac{r\cdot x\cdot z\cdot \exp (-z)}{x*/p+x}$ at $w_1=A_1\cdot e^{-x\left(a(1-p)+\rho p\cdot \left(1-\frac{\alpha \cdot z}{x*+px}\right)\right)}$ and $w_2=A_2\cdot \left(1-\frac{\alpha \cdot z}{x*+px}\right)\cdot e^{-x\left(p\cdot \left(1-\frac{\alpha \cdot z}{x*+px}\right)+\varphi (1-p)\right)}$.

The traditional method for finding the stability area of fixed points is based on the fact that the roots of the characteristic equation of system (6) belong to the circle |λ| < 1:

λ³ − S λ ² + H λ − K = 0,

(14)

where S, H, and K are three invariants of the Jacobian matrix (13). The conditions for codimension one bifurcations of system (6) are found based on coefficients of its characteristic polynomial (14) and have the form [48]:

• Transcritical bifurcation, T, λ = 1 (the dominant root is real and is equal to 1): H = S + K − 1. At the transcritical bifurcation line, the stationary solutions of the system exchange their stability. However, at λ = 1, a saddle-node bifurcation (SN) can also occur when, in addition to the existing fixed points of the system, an unstable equilibrium arises, which then splits into stable and unstable attractors.
• Period-doubling bifurcation, PD, λ = −1 (the dominant root is a real and equal to −1): H = −S − K − 1.
• Neimark–Sacker bifurcation, NS, λ = e ^±iφ (the real root modulus is less than 1, and the absolute values of the complex conjugate roots are equal to 1): H = SK − K² + 1.

4.1. Boundaries of the Stability Regions of the Fixed Points of System (6)

Based on the methodology presented above, the stability region of fixed point (7) is shown to be bounded by the transcritical bifurcation lines A₁ = 1 and A₂ = 1, at which the trivial (7) and semitrivial (8) and (9) solutions exchange their stabilities. Consequently, in the parameter plane (A₁, A₂), the stability domain of the trivial fixed point is a unit square. Obviously, the extinction of both phytoplankton haplotypes with the following zooplankton elimination is observed for any values of the other parameters if $A_{1,2}<1$. Dividing the square along the bisector of the first coordinate angle, one can separate the parametric region where phytoplankton is represented only by the toxic haplotype from the one where only by the non-toxic haplotype by the time of phytoplankton elimination. The ratio of the parameters A₁ and A₂ determines this.

The stability domain of monomorphic fixed point $\overline{p}=0$, $\overline{x}={\overline{x}}_1=\left(\ln A_1\right)/a$, $\overline{z}=0$ (8) is formed by transcritical bifurcation lines A₁ = 1 and A₂ = A₁^φ/a and boundary A₁ = e², whose crossing leads to period-doubling bifurcation. For the other monomorphic solution $\overline{p}=1$, $\overline{x}={\overline{x}}_2=\ln A_2$, $\overline{z}=0$ (9), the boundaries of its stability area are set by the T curves: A₂ = 1, A₁ = A₂^ρ and A₂ = exp(x*/(r − 1)) at r > 1 + x*/2. At r < 1 + x*/2, the stability region is bounded by the transcritical bifurcation lines A₂ = 1 and A₁ = A₂^ρ and the period-doubling bifurcation line PD: A₂ = e². At r > 0, the stability domain of solution (9) is formed only by lines of transcritical bifurcations, and the larger the r value, the narrower the stability area. In particular, at very high r values, the boundary given by A₂ = exp(x*/(r − 1)) tends to A₂ = 1.

For polymorphic point (10), where both phytoplankton haplotypes coexist without zooplankton (z), one of the roots of the characteristic equation is ${\lambda }_1=r\left(\varphi \ln A_1-a\ln A_2\right)/\left(x*(\varphi \rho -a)+\varphi \ln A_1-a\ln A_2\right)$. In addition to A₂ = A₁^φ/a and A₁ = A₂^ρ, we can obtain another transcritical bifurcation line $\ln A_2=\frac{\varphi }{a}\ln A_1+\frac{x*(\varphi \rho -a)}{a(1-r)}$ that divides the stability regions of solutions (10) and (12). The Neimark–Sacker bifurcation does not occur in this case. Solution (10) loses its stability through a period-doubling bifurcation when the condition does not hold:$\ln A_1\ln A_2(a-\varphi \rho )+a\ln A_2(2-2\rho -\rho )+\ln A_1(\varphi \ln A_1+2a-2\varphi )+4\varphi \rho -4a\le 0$.

From the point of view of evolutionary processes, the main interest is in the stability and mechanisms of its loss by nontrivial fixed points (11) and (12), which correspond to the coexistence of phytoplankton and zooplankton; in the first solution, the phytoplankton population is monomorphic, and in the second solution it is polymorphic. The Jacobian (13) at point (11) takes the form:

$$J=\left(\begin{array}{ccc}{w}_1/w_2& 0& 0\\ F{x^{\prime }}_p& w_2+\overline{x}\cdot {(w_2{)}^{\prime }}_x& \overline{x}\cdot {(w_2{)}^{\prime }}_z\\ F{z^{\prime }}_p& \frac{x*\cdot \overline{z}}{(x*+\overline{x})\overline{x}}& 1-\overline{z}\end{array}\right).$$

It corresponds to the following characteristic equation:

$\left(\lambda -\frac{w_1}{w_2}\right)\left(\left(\frac{\alpha \overline{z}\overline{x}(x^{*}(\alpha \overline{x}-\alpha \overline{z}+\overline{x}+x^{*})+\alpha {\overline{x}}^2)}{(\overline{x}+x^{*}-\alpha \overline{z}){(\overline{x}+x^{*})}^2}-\overline{x}\right)(1-\overline{z}-\lambda )-\frac{\alpha x^{*}\overline{x}\overline{z}(\alpha (x^{*}+\overline{x}-\overline{z})-\overline{x}-x^{*}}{(\overline{x}+x^{*}-\alpha \overline{z}){(\overline{x}+x^{*})}^2}\right)=0,$ provided that $A_2\cdot \left(1-\frac{\alpha \cdot \overline{z}}{x*+\overline{x}}\right)\cdot e^{-\overline{x}\left(1-\frac{\alpha \cdot \overline{z}}{x*+\overline{x}}\right)}=1$. One of the eigenvalues of the polynomial is ${\lambda }_1=\frac{w_1}{w_2}$, and the other two are the roots of the quadratic equation. Based on system (11), we obtain the line of transcritical bifurcation, separating fixed points (11) and (12), which in parametric form is the following:

$$\mathrm{T}\text{:} \{\begin{array}{l}{A}_1=e^{\frac{t\rho \cdot \left(x*+t-\alpha \ln \frac{r\cdot t}{x*+t}\right)}{x*+t}}\\ A_2=\frac{x*+t}{x*+t-\alpha \ln \frac{r\cdot x}{x*+t}}{e}^{\frac{t\cdot \left(x*+t-\alpha \ln \frac{r\cdot t}{x*+t}\right)}{x*+t}}\end{array}.$$

(15)

To determine the nature of stability loss by solutions (11) and (12), it is reasonable to consider the behavior of bifurcation lines and the stationary solution $\overline{z}$ in the plane $\left(\begin{array}{cc}{A}_2,& z\end{array}\right)$ provided that $\overline{p}=1$, $\overline{x}=\frac{x*}{r\cdot \exp (-\overline{z})-1}$, and $\overline{p}=\frac{x*}{(r\cdot \exp (-\overline{z})-1)\overline{x}}$, $\overline{x}=\frac{\alpha \rho \overline{z}\exp (\overline{z})+r\ln A_1}{ar}-\frac{x*r\exp (\overline{z})(a-\rho )}{ar(\exp (\overline{z})-r)}$, respectively. Note that such a graph is essentially a bifurcation diagram showing the values of the fixed point $\overline{z}$ with corresponding bifurcation transitions (Figure 3).

Figure 3 shows that, depending on the value of the parameter A₂, the number of fixed points corresponding to both the reduced community represented by the non-toxic phytoplankton haplotype and zooplankton (Figure 3a) and the complete community (Figure 3b,c) changes: a specific value of A₂ can correspond to one to three values of $\overline{z}$. Note that the stability loss of one of the roots is always realized through a cascade of period-doubling bifurcations: the solution line $\overline{z}$ intersects the PD curve at A₂ = A₂^PD. In the range where the system has three solutions, one can see that a tangent bifurcation occurs; the SN curve intersects the solution line, while the third root lies in the instability region. The increase in A₂ values leads to the emergence of both direct (at A₂ = A₂^NS+) and inverse (at A₂ = A₂^NS−) Neimark–Sacker bifurcations. First, the fixed point loses stability by forming an invariant curve, which collapses with higher values of the bifurcation parameter, and only the attracting equilibrium remains. If the tangent bifurcation occurs after the collapse of the invariant curve or before the Neimark–Sacker bifurcation emergence, then bistability is observed; there are two stable nontrivial equilibria in the system, and which of them will be attractive is determined by the initial conditions (Figure 3a,b). In particular, Figure 3a,b show that two fragments of the solution curve lie in the stability area. If tangent bifurcation occurs before the collapse of the invariant curve, then an unstable fixed point undergoes a saddle-node bifurcation; two of the three nontrivial roots are located outside the stability region (Figure 3c).

4.2. Stability Regions of Fixed Points of System (6)

The stability analysis presented above can be generalized by constructing stability domains for all fixed points of system (6) in the parameter plane (A₁, A₂) (Figure 4). We numerically found the period-doubling and the Neimark–Sacker bifurcation lines for solutions (11) and (12).

Figure 4a shows a parametric portrait of the community without zooplankton; there is a large area of parameter values where phytoplankton is polymorphic (both haplotypes coexist). The picture is symmetrical, as the coefficients ρ and φ, characterizing the competitive relationships between different genetic groups of phytoplankton, are the same, and the parameter a, corresponding to self-limitation, equals 1. At r = 0, the stability region coincides with the results published in [43] and demonstrates the possible scenarios of the evolution of competing phytoplankton haplotypes without zooplankton. Furthermore, by increasing r, we can analyze the changes occurring in the community with the presence of another trophic level. In particular, the growth of r narrows the attraction areas of solutions (9) and (10). Note that the stability loss in Figure 4a–c,e,f occurs through a period-doubling bifurcation. However, at higher values of r, solutions (11) and (12) can lose stability according to the Neimark–Sacker scenario; Figure 4d shows characteristic discontinuities within the stability region: these boundaries were found numerically. Here, because of the bistability associated with several attracting solutions, the initial condition can influence the type of parametric portrait. The increase in prey consumption by the predator expands the range of parameter values where solution (12) is stable, while increasing the half-saturation constant leads to reverse changes due to the emergence of a vast zone where zooplankton goes extinct and only phytoplankton remains, maintaining its genetic variety.

At φρ < a, the stability regions of solutions (8) and (9) do not intersect (Figure 4, domains 2 and 3). However, at φρ = a, they have a common border A₁ = A₂^ρ (Figure 5a), and at φρ > a, their stability domains overlap each other (Figure 5b), which leads to bistability and multistability. Figure 5 shows the arrangement features of parametric areas corresponding to different community development scenarios depending on zooplankton reproduction and competitive relationships between different phytoplankton haplotypes. The lines of transcritical bifurcations: A₁ = 1, A₂ = A₁^φ/a, A₂ = 1, A₁ = A₂^ρ, A₂ = exp(x*/(r − 1)) and curve (15) form the areas.

If φρ ≥ a, polymorphism in the phytoplankton population does not arise in the community without a predator, i.e., the more adapted haplotype always displaces the less adapted one, or the bistability of monomorphic states occurs when the current ratio of genetic group densities determines the outcome of evolution. The presence of a predator leads to the possibility of the successful coexistence of competing haplotypes in some regions of the parameter space (Figure 5c). In general, in the community with a predator, the scenario, corresponding to one haplotype displacing the other, presents, but here, the displacement conditions are associated with both intraspecies competition between phytoplankton haplotypes and predator–prey interactions. At the same time, changing the values of the parameter ρ affects the type of parametric portraits of the system either without a predator or with small r. The higher the reproductive potential of the predator r, the wider the area with the polymorphic prey population. Note that the bistability of monomorphic states corresponding to different genetic compositions in the phytoplankton population also occurs in the community with a predator. However, if the phytoplankton haplotype with the toxicity trait displaces the other haplotype, the zooplankton will go extinct in the long run. If the haplotype without the toxicity trait displaces its competitor, then the coexistence of predator and prey will continue. Thus, variations in the current genetic composition of phytoplankton, for example, due to the influence of environmental factors, can change the direction of evolution, entailing a change in the community development scenario.

4.3. Evolutionary Scenarios for Community Development

The results obtained for stability analysis allow us to determine the evolutionary scenarios for the planktonic community depending on the values of parameters characterizing the fitnesses of phytoplankton haplotypes. The generalization of possible evolutionary development scenarios is presented in

Table 1. Evolutionary development scenarios for the plankton community, depending on the parameter values.

Parameter Value Ranges		Evolutionary Scenarios of Community Development
$r\le 1+x*/2$
$\rho \varphi \le a$	$A_1<A_2{}^{\rho }$$, A_2>1$	Monomorphism 1, zooplankton becomes extinct
	$A_2<A_1{}^{\varphi /a}$$, A_1>1$	Monomorphism 2, zooplankton becomes extinct
	$A_1>A_2{}^{\rho }$$, A_2>A_1{}^{\varphi /a}$$, A_1>1$$, A_2>1$	Polymorphism 3, zooplankton becomes extinct
$\rho \varphi >a$	$A_2>A_1{}^{\varphi /a}$$, A_2>1$	Monomorphism 1, zooplankton becomes extinct
	$A_1>A_2{}^{\rho }$$, A_1>1$	Monomorphism 2, zooplankton becomes extinct
	$A_1<A_2{}^{\rho }$$, A_2<A_1{}^{\varphi /a}$$, A_1>1$$, A_2>1$	Bistability of monomorphic states: displacement of one haplotype by another depends on initial conditions
$r>1+x*/2$
$\rho \varphi \le a$	$A_2>1$$,A_1<A_2{}^{\rho }$$,\ln A_2<x*/(r-1)$	Monomorphism 1, zooplankton becomes extinct
	$A_1>1$$,A_1>A_2{}^{a/\varphi }$	Monomorphism 2, zooplankton becomes extinct
	$A_1>A_2{}^{\rho }$$,A_1<A_2{}^{a/\varphi }$$,\ln A_2<(\varphi \ln A_1+x*(\varphi \rho -a)/(1-r))/a$	Polymorphism 3, zooplankton becomes extinct
	$\ln A_2>\frac{x*}{r-1}$$, A_2>A_2*$$\mathrm{at} A_1=A_1*$ 4	Monomorphism 1, zooplankton survives and develops
	$\ln A_2>(\varphi \ln A_1+x*(\varphi \rho -a)/(1-r))/a$$,A_2<A_2*$$\mathrm{at} A_1=A_1*$ 4	Polymorphism 3, complete “phytoplankton–zooplankton” community
$\rho \varphi =a$	$A_2>1$$, A_1<A_2{}^{\rho }$$, \ln A_2<x*/(r-1)$	Monomorphism 1, zooplankton becomes extinct
	$\ln A_2>\frac{x*}{r-1}$$, A_2>A_2*$$\mathrm{at} A_1=A_1*$ 4	Monomorphism1, zooplankton survives and develops
	$A_1<A_2{}^{\rho }$$, A_2<A_2*$$\mathrm{at} A_1=A_1*$ 4	Polymorphism 3, complete “phytoplankton–zooplankton” community
	$A_1>1$$, A_1>A_2{}^{\rho }$	Monomorphism 2, zooplankton becomes extinct
$\rho \varphi >a$	$A_2>1$$, A_2>A_1{}^{\varphi /a}$$, \ln A_2<x*/(r-1)$	Polymorphism 3, zooplankton becomes extinct
	$A_1<A_2{}^{\rho }$$,A_2<A_1{}^{\varphi /a}$$, \ln A_2<x*/(r-1)$	Bistability of monomorphic states: displacement of one haplotype by another depends on initial conditions, zooplankton becomes extinct
	$\ln A_2<x*/(r-1)$$, A_2>A_1{}^{\varphi /a}$, $A_2>A_2*$$\mathrm{at} A_1=A_1*$ 4	Bistability of monomorphic states: displacement of haplotype x2 leads to zooplankton extinction, otherwise zooplankton survives and develops
	$\ln A_2<x*/(r-1),$$A_2>A_2*$$\mathrm{at} A_1=A_1*$ 4	Monomorphism 1, zooplankton survives and develops
	$A_1>A_2{}^{\rho },$$A_2>A_1{}^{\varphi /a},$ $A_2>A_2*$$\mathrm{at} A_1=A_1*$ 4	Monomorphism 2, zooplankton becomes extinct
	$A_2<A_1{}^{\varphi /a}$$,\mathrm{at}{A}_1=A_1*$	Polymorphism 3, complete “phytoplankton–zooplankton” community

Monomorphism ¹: phytoplankton haplotype without toxicity trait x₂ displaces haplotype x₁. Monomorphism ²: phytoplankton haplotype with toxicity trait x₁ displaces haplotype x₂. Polymorphism ³: haplotypes x₁ and x₂ coexist.

.

⁴$A_1*=e^{\frac{t\rho \cdot \left(x*+t-\alpha \ln \frac{r\cdot t}{x*+t}\right)}{x*+t}}$, $A_2*=\frac{x*+t}{x*+t-\alpha \ln \frac{r\cdot x}{x*+t}}{e}^{\frac{t\cdot \left(x*+t-\alpha \ln \frac{r\cdot t}{x*+t}\right)}{x*+t}}$, and t is a parameter.

5. Bistability and Multistability in System (6)

At φρ < a and $\overline{z}\ne 0$, both monomorphic (11) and polymorphic (12) solutions demonstrate bistability (Figure 3). Let us examine this bistability in more detail by supplementing the graphs in Figure 3 with bifurcation diagrams constructed for specific initial conditions (Figure 6).

The bifurcation diagrams show sharp changes in dynamic modes, which is the identification of bistability and/or multistability. In particular, the intersection of line SN with the solution curve corresponds to a saddle-node bifurcation; as a result, we observe a transition from one dynamic mode to another. The dynamic behavior picture for communities represented by zooplankton with monomorphic phytoplankton without toxicity traits (A₁ = 0) and zooplankton with polymorphic phytoplankton (A₁ ≠ 0) is similar (Figure 6).

The evolution of the dynamic modes of model (6) with increasing values of parameter A₂ can be described as follows: stable solution 1₁ exists and then loses its stability via the Neimark–Sacker bifurcation after the invariant curve collapses, and solution 1₁ becomes stable (Figure 6a,c). Against the background of this fixed point, a tangent bifurcation results in the appearance of an unstable solution 1₂ and a stable solution 1₃, that loses stability through a cascade of period-doubling bifurcations with increasing A₂. Note that in the vicinity of the tangent bifurcation at higher A₂, solution 1₁ continues to exist and attract until the merger of 1₁ and 1₂. As a result, only solution 1₃ remains. If stable solution 1₁ does not yet lose stability via the Neimark–Sacker bifurcation, but attractive solution 1₃ appears due to saddle-node bifurcation, then the bistability of nontrivial fixed points is observed (Figure 6a–c).

On the other hand, the tangent bifurcation leading to the appearance of 1₃ may occur before the inverse Neimark–Sacker bifurcation, resulting in the invariant curve collapsing and only the attracting fixed point 1₁ remaining. In this case, the stable solution 1₃ coexists with quasi-periodic oscillations Q₁ (Figure 6d). However, the stable nontrivial fixed point 1₁ can also coexist with irregular oscillations. In particular, Figure 6b–d show that the oscillations denoted Q* arise before losing stability of point 1₁ via the Neimark–Sacker bifurcation: irregular dynamic modes occur before the bifurcation value A₂^NS+. As a result, 1₁ and Q* can coexist (Figure 7). We constructed basins of attraction in the space of initial conditions (p₀, z₀) to study the influence of genetic composition on the scenarios of community evolution.

At low values of the growth rate of non-toxic phytoplankton, for any initial conditions from the feasible region, the system’s dynamics stabilize. Then, a cycle with period 11, which is the solution of the 11-time iterated system (6), occurs (Figure 7a). At higher reproductive potential, along with a stable state and cycle with period 11, extensive regions of irregular Q* dynamics arise (Figure 7b). Note that we counted 15 thousand iterations of mapping (6) to check whether the quasi-periodic dynamics is not a stabilized 11-cycle. At a higher growth rate of the phytoplankton haplotype without toxicity, the basins of attraction demonstrate the coexistence of quasi-periodic oscillations Q* and the single fixed point 1₁.

A further increase in A₂ reduces the attraction basin of fixed point 1₁ (Figure 7d), and it loses stability via the Neimark–Sacker bifurcation (Figure 7e). As a result, the system demonstrates quasi-periodic oscillations with different phases for the various initial conditions (Figure 7c). To distinguish different phases of quasi-periodic oscillations, we selected one time series/trajectory and then calculated the correlation coefficients for it and others in each point of the parametric plane. It allowed the identification of areas with similar and distinct dynamics in region Q* (Figure 7c). Note that the loss of stability of point 1₁ leads to the appearance of an unstable curve Q₁, merging with Q*. This process is presented in phase portrait 7f using an initial condition near point 1₁. Consequently, subcritical Neimark–Sacker bifurcation occurs. In other words, the attraction basin of Q₁ collapses, and Q* captures the entire phase space (Figure 7c).

At the bifurcation value A₂ = A₂^SN, a saddle-node bifurcation appears, leading to the emergence of a new stable point 1₃, which is attractive in a significant part of the initial condition space (Figure 7e). By increasing A₂, we observe the bistability of two nontrivial fixed points, one of which is the result of a saddle-node (tangential) bifurcation, and its basin of attraction is substantially larger than the other (Figure 7e). Note that with further growth of parameter A₂ values, this stable point captures the entire space of initial conditions and then loses stability through a cascade of period-doubling bifurcations.

To investigate the multistability depicted in Figure 7, we constructed dynamic mode maps in the ranges of parameter values r and A₂ where quasi-periodic dynamics emerge (Figure 8). The maps were generated as follows: at each point (corresponding to a single pixel), 5000 iterations of the mapping were calculated. Based on the results of the last 500 steps, the period of the limit cycle was determined, and the point was colored according to the period obtained.

Figure 8 illustrates abrupt transitions between different dynamic modes, namely, multistability caused by the coexistence of a stable fixed point 1₁ and quasi-periodic oscillations Q* or 1₁/Q₁ with 1₃ (Figure 6b–d). Consequently, a variation in initial conditions can change the area with quasi-periodic and periodic dynamics. As observed, the region of quasi-periodic dynamics located within the stability domain of the nontrivial fixed point is bounded from above by two intersecting curves, forming an angle, which is also indicative of multistability, as the system “abruptly” transitions from Q₁ oscillations to stable point 1₃, as depicted in the bifurcation diagram shown in Figure 6d. In other words, a variation in initial conditions can both expand and contract the area of quasiperiodic dynamics. Accordingly, the “upper” boundary of the quasi-periodic oscillation region, generated as the stability loss by fixed point 1₁, must be higher than that represented in Figure 8, smoothing out the angle. Conversely, the lower boundary can be higher, as point 1₁ loses stability later than Q* emerges.

Thus, when φρ < a in the case of polymorphism, the variation in the current composition of the community can alter the observed dynamics mode, which changes the densities of different phytoplankton haplotypes. However, the evolution direction remains the same. At φρ < α in the case of monomorphism with a non-toxic haplotype of phytoplankton, variations in the current densities of algae and zooplankton can also lead to a shift in the dynamics mode.

At φρ ≥ a, system (6) reveals the bistability of monomorphic (8) and (9) states (Figure 5b,d, Table 1), as well as monomorphic (8) and polymorphic (11) solutions (Figure 5d–f, Table 1). We consider this bistability using genetic composition maps of the prey population with various initial conditions for the parametric portrait shown in Figure 5e (Figure 9). As seen, variations in initial conditions with specific values of population parameters can change the evolution direction and, consequently, the dynamics mode (Figure 9e,f). Such multistability is caused by the system having several attracting solutions simultaneously (Figure 9c–h). In the case of the stability of two monomorphic fixed points, the “bistability trap” occurs, where a more promising genetic form cannot displace a weaker competitor (Figure 9c) [49]. However, a drop in population size due to an external factor effect known as a “population bottleneck” [50] can lead to random changes in genetic composition and the transition of the population from one monomorphic state to another more promising but destabilizing dynamic.

Variation in current population sizes in the complete community can lead to the displacement of a phytoplankton haplotype with toxicity traits (Figure 9d,h) due to the bistability of monomorphic and polymorphic solutions. Specifically, Figure 9g shows that in the case of zooplankton presence, system (6) has three fixed points, two of which are attractive (Figure 9h). Figure 9 also demonstrates that the presence of predators contributes to the maintenance of genetic polymorphism in the prey population, particularly the fixation of the haplotype with toxicity traits. In this context, our findings are in agreement with the hypothesis proposed in the study by X. Irigoien et al. (2005) [51]. This hypothesis is based on natural observations and suggests that blooming species are capable of avoiding predation by microzooplankton during the early stages of blooms.

6. Dynamic Modes of Model (6)

System (6) comprises a pair of Ricker models coupled by the Holling type-II function. The dynamics of the Ricker model are well known: as the parameter characterizing the population growth rate increases, a cascade of period-doubling bifurcations leading to chaos is observed [42]. However, in the proposed model, the Neimark–Sacker scenario of stability loss arises. Note that the “predator–prey” interaction gives quasi-periodic oscillations. The changes in the community development scenarios during evolution, depending on the values of system (6) parameters, are shown in Figure 10.

Figure 10 illustrates that a nontrivial fixed point can lose stability through the Neimark–Sacker scenario and a cascade of period-doubling bifurcations. If the system has a single nontrivial solution and the value of predator growth rate is such that the Neimark–Sacker bifurcation has not yet occurred, an increase in the birth rate of zooplankton results in the appearance of a two-cycle, which evolves toward chaos through a cascade of period-doubling bifurcations. This scenario is possible due to the loss of stability of the fixed point when the growth of parameter r leads to crossing the bifurcation boundary on which one of the eigenvalues of the Jacobian matrix of system (6) becomes equal to (−1). If the Neimark–Sacker bifurcation occurs, then with higher values of the bifurcation parameter, the only nonzero solution loses stability through this scenario, forming an invariant curve that then “collapses”, and the nontrivial fixed point becomes stable again. Furthermore, a saddle-node bifurcation occurs, and a new solution arises, losing stability through a cascade of period-doubling bifurcations. This dynamic mode evolution is depicted in Figure 6a. The presence of inverse Neimark–Sacker bifurcation suggests that an increase in the prey growth rate can stabilize irregular predator dynamics. With higher values of r, when the system has three nontrivial fixed points due to saddle-node bifurcation, the dynamics are similar but complicated by multistability (Figure 10c).

Analyzing the impact of parameter values on the system’s dynamic behavior shows that an increase in the coefficient characterizing the predator’s pressure on the prey population leads to the emergence of long-periodic delayed oscillations at a lower growth rate of the predator population (Figure 10e,f). Higher values of the half-saturation constant significantly expand the stability domain of the nontrivial fixed point, resulting in quasi-periodic oscillations at higher values of the predator’s reproductive potential. Note that the non-smooth Neimark–Sacker bifurcation line suggests the maintenance of multistability due to the coexistence of several nontrivial fixed points in the system.

Maps constructed in the parameter space (r, A₂) allow us to analyze the changes in community dynamics caused by predator–prey interactions (Figure 10d–f). With a zero reproductive potential of the predator (r = 0), the abscissa axis corresponds to the dynamic modes of a community, represented by competing genetic groups of phytoplankton without the predator. If A₁ = 0, this community consists of a non-toxic haplotype, and the Ricker model describes its dynamics. At A₂ < 1, the non-toxic phytoplankton becomes extinct; at 1 < A₂ < e², the population density stabilizes; and at A₂ > e², the dynamics become complex via a cascade of period-doubling bifurcations. The presence of a predator with a low growth rate in the community of competitors does not change the phytoplankton dynamics because the predator goes extinct due to insufficient reproduction. Furthermore, by increasing the predator’s reproductive rate (r), we can compare the system’s dynamics with and without the predator. In particular, the emergence of quasi-periodic oscillations is the result of predator–prey interactions, as the nontrivial fixed point of two coupled Ricker models does not lose stability according to the Neimark–Sacker scenario [43].

Figure 10d–f illustrate that an increase in the values of coefficient r does not always change the nature of the community dynamics that coincides with the dynamic mode of the prey population in the absence of predators. In other words, predator dynamics follow prey one. At high values of the coefficient α, characterizing the predator’s pressure on the prey, the two-cycle is replaced by a stable state due to prey consumption by the predator. In the case of cycles with higher periods, an increase in predator reproductive potential results in a cascade of period-doubling bifurcations: the higher the growth rate of the phytoplankton haplotype without toxic traits, the lower the parameter r value giving chaotic dynamics. This phenomenon is associated with the superposition of oscillation periods in predator and prey populations and corresponds to situations when prey population dynamics are determined by predator one. In the case of the emergence of quasiperiodic oscillations due to predator–prey interactions, prey dynamics adapt to predator one.

7. Comparison of Model Community Dynamics with Experimental Observations

Here, we analyze the correspondence between the dynamic patterns of the model (6) and the qualitative characteristics of the experimental dynamics of the phyto-zoo plankton community described in the literature. Considering the population cycles of the predator (zooplankton) and its prey (the non-toxic haplotype of phytoplankton) separately, one can observe dynamics typical for the predator–prey model without evolution, as described in [38] for the chemostats with rotifers and genetically homogeneous algae: the predator’s dynamics lag behind the prey by approximately one-quarter of a period. Model (6) demonstrates such behavior (Figure 11a). Note that Figure 11a depicts the dynamics of phytoplankton represented only by the haplotype without toxicity trait because the prey is monomorphic. Dynamics with time lag are also observed in the case of polymorphism (Figure 11b,c). In particular, the cycle with period 4 is anti-phase oscillations of the predator and prey, represented by the non-toxic haplotype of phytoplankton, because the predator density peak coincides with its prey minimum, and the predator’s dynamics lag behind the prey by one time step, i.e., one-quarter of a cycle. Figure 11b shows a similar lag of variable p_n that corresponds to the density of the non-toxic haplotype in the total phytoplankton biomass.

Under conditions of high predator consumption rates, the dynamics of the community show long-period antiphase oscillations between predators and prey (Figure 11d), reminiscent of those observed during the evolution of prey in chemostats with different algal genotypes [38]. The trajectories in Figure 11b,d represent long-period oscillations with a time delay, similar to self-oscillations in the classical Lotka–Volterra model. Equation (6) also reveal dynamic patterns resembling batch dynamics, as seen in continuous-time systems (Figure 11f).

In Figure 11e, the quasi-periodic dynamics demonstrate antiphase oscillations in phytoplankton density. Considering the heterogeneity of phytoplankton, represented by two genetic forms, these oscillations can offset each other in the overall phytoplankton biomass, creating the appearance of a nearly steady state condition. Such dynamics, characterized by constant prey density and fluctuating predator one, are referred to as “cryptic cycles” and have been observed in experiments [39]. The example depicted in Figure 11e illustrates a potential mechanism for the emergence of cryptic cycles in communities with heterogeneous prey, where rapid evolution conceals trophic interactions among species. Note that in the case of stability loss of a fixed point via the Neimark–Sacker bifurcation, leading to an unstable invariant curve, the transient dynamics represent predator fluctuations with nearly constant prey density for an extended period. However, over time, the dynamics of the prey adjust to that of the predator, resulting in oscillations between the interacting species.

8. Conclusions and Discussion

Discrete-time equations naturally describe time lag and are an alternative to differential equations with delay and piecewise continuous functions in modifications of the classical phytoplankton–zooplankton interaction model. In line with the traditional ideology in this field, this study introduces a discrete-time eco-genetic model of a phytoplankton-zooplankton community, considering the genetic diversity of phytoplankton. The community is assumed to include zooplankton and two competing phytoplankton haplotypes—one with a toxicity trait and one without. The interaction between the two genetic phytoplankton groups is described by the Ricker competition model, taking into account the limitation of growth for each genotype due to resource availability, for example, nutrients, oxygen, light, etc. We used the Holling Type II functional response with predator saturation to model phytoplankton consumption by zooplankton. The birth rate and survival of zooplankton depend on their feeding success. The density-dependent regulation of the growth of the zooplankton population implicitly includes mortality due to the increase in the amount of toxic substances secreted by phytoplankton caused by the high density of the zooplankton. We assume that zooplankton (predators) predominantly consume the non-toxic haplotype of phytoplankton (prey) because phytoplankton can defend against predating.

The results obtained in this study show that within a significant part of the parameter space, the community primarily collapses due to prey overconsumption by the predator: only the toxic haplotype of phytoplankton survives, and its dynamics depend on the values of its population parameters. Note that the local extinction of a portion of the community at a specific spatial location is not a final state because, over time, non-toxic phytoplankton and zooplankton will again appear here due to their migration from neighboring areas.

The stability analysis has shown that the stability loss of a nontrivial fixed point, corresponding to the coexistence of polymorphic phytoplankton with zooplankton, can occur through a cascade of period-doubling bifurcations and the Neimark–Sacker bifurcation. The emergence of the Neimark–Sacker bifurcation leading to quasi-periodic oscillations is only feasible in the community model, whereas in the original Ricker equation describing the density of each species, destabilization of dynamics occurs via period-doubling bifurcation. The proposed eco-genetic model of planktonic community dynamics reveals long-period oscillations consistent with the results of field experiments. Variation in values of parameters characterizing phytoplankton and zooplankton can lead to significant changes in the dynamic mode of the community, for example, transitions from regular to quasi-periodic dynamics or from oscillations to stationary dynamics. Note that quasi-periodic dynamics can arise at not high phytoplankton growth rates, giving stable or periodic community dynamics. In this parametric area with multistability, the transition from regular dynamics to quasi-periodic and vice versa can occur due to variations in initial conditions caused by external factor impact changing the current species densities, which shifts the system into the attraction basin of another dynamics mode. In the multistability domain, a shift in both the dynamics mode and the evolution direction, changing community composition due to variations in initial conditions, is possible. Note that unpredictability in multispecies competition is discussed in the study by J. Huisman and F. Weissing [52]. The research demonstrates that the system may yield multiple alternative outcomes, with the dynamics leading to these alternatives possibly showcasing transient chaos. Moreover, the basins of attraction associated with these alternative outcomes exhibit an intermingled fractal geometry. As a result, the prediction of the winners in a multispecies competition community becomes a challenge [52].

From the perspective of evolutionary dynamics, the relationship between parameters that characterize competition among different phytoplankton haplotypes and self-limiting processes is more interesting. With φρ < a, the outcome of evolution can vary significantly depending on the competing genotypes’ fitnesses, leading to either monomorphism establishment or polymorphism maintenance. Here, a variation in the current densities of the community component can shift the observed dynamics, which changes the ratio for phytoplankton haplotype frequencies and consequently their total abundance, while the evolution direction remains unchanged. At φρ < a with phytoplankton monomorphism represented by a non-toxic haplotype, fluctuations in the current densities of algae and zooplankton can also shift the dynamics mode. In general, we can conclude that when φρ < α, the community composition and the evolution direction are determined by the fitness of phytoplankton haplotypes and the “predator pressure” alongside the carrying capacity. “Predator pressure” influences the non-toxic phytoplankton haplotype density and can be considered an indirect reduction in fitness. As a result, the haplotype with higher fitness fixes in the community. The variation in the current population sizes of community components can change the dynamics mode but not the evolution direction. In particular, oscillations can be replaced by stabilization, and vice versa. Such scenarios are possible due to saddle-node bifurcations, the occurrence of which primarily depends on the values of parameters x*, r, α, and A₂, characterizing the “predator–prey” interactions. As a result, population densities in the community can stabilize, reaching higher or lower values depending on the initial conditions with the same evolutionary direction.

With φρ ≥ a, a phenomenon known as the “bistability trap” can occur, where a more advantageous genotype cannot displace its genetically weaker competitor [49]. However, population decline resulting from external factors, leading to random processes referred to as “population bottlenecks”, can induce random changes in genetic composition, causing transition processes from one monomorphic state to another with greater fitness. It is important to note that in the case of a scenario of a bistability trap without a predator, the coexistence of both competing phytoplankton haplotypes is not possible [43]. However, the predator’s presence in the system maintains the genetic polymorphism of the prey, preventing the displacement of the haplotype with toxicity traits. These findings align well with the hypothesis proposed based on field observations [51], suggesting that blooming species are capable of evading microzooplankton predation during the early stages of blooms. In other words, a variation in the current community composition can lead to a shift in the evolution direction, where the haplotype with toxicity traits either becomes fixed or goes extinct.

In general, the eco-genetic discrete-time model of the planktonic community proposed in this study is relatively simple but demonstrates oscillations that are typical for predator–prey interactions. This system reveals dynamic modes that reflect the basic properties of experimental dynamics: one can observe the predator lagging behind the prey by approximately one-quarter of the period; such behavior is typical for a planktonic community without considering evolution. Taking into account the genetic heterogeneity of phytoplankton, even in the case of two genetic forms with and without the toxicity trait, allows us to observe both long-period antiphase oscillations of predator and prey and cryptic cycles. Cryptic cycles manifest as scenarios where the density of prey remains nearly constant while the density of predators fluctuates, and these phenomena are driven by rapid evolution that conceals trophic interactions among species.