Adaptive Dynamics of a Stoichiometric Phosphorus–Algae–Zooplankton Model with Environmental Fluctuations

Abstract

We present a stochastic evolutionary phosphorus–algae–zooplankton model with phosphorus recycling and originally investigate the patterns and outcomes of adaptive changes in algal cell size, under the influence of environmental fluctuations. The threshold that determines whether the stochastic model will ecologically persist or not is first obtained. We then introduce fitness functions with stochastic fluctuations and obtain the evolutionary conditions for continuously stable strategy (CSS) and evolutionary branching, confirmed by numerical simulation. Our results predict that environmental fluctuation could drive algal evolution toward smaller cell size. Algal cell size varies significantly with phosphorus input in the presence of zooplankton, but has no response to the changing phosphorus inflow without zooplankton, and evolutionary branching will never occur without zooplankton. With the existence of zooplankton that has a fixed trait, evolutionary branching occurs with small environmental fluctuation and moderate phosphorus inflow, and the existence of environmental fluctuation could narrow the cell size difference between the newly emerging algal species, while large fluctuation or extreme phosphorus inflow will result in CSS. Moreover, environmental fluctuation potentially benefits algal biodiversity in eutrophic environments, and oligotrophication inhibits algal diversity. For the coevolution of algae and zooplankton, evolutionary cycling could appear, i.e., algal cell size and zooplankton body size can coevolve to a stable limit cycle (the Red Queen dynamics) in an eutrophic environment. In oligotrophic or moderate phosphorus environments, the influence of environmental fluctuation on algal evolution in the coevolution process is similar to the scenario that algae evolves only.

Appendices

Appendix A: Proof of Theorem 2

Proof

Applying the It$\hat{\mathrm {o}}$’s formula to the second and third equations of A and Z in model (2), integrating both sides from 0 to t and then dividing t on both sides yield

$$\begin{aligned}&\frac{1}{t}\ln \frac{ A(t)}{A(0)}=\alpha \mu (x)\frac{1}{t}\int _{0}^{t}P(s)\mathrm {d}s-c(x,y)\frac{1}{t}\int _{0}^{t}Z(s)\mathrm {d}s-\eta (x,\rho _{1}) +\sigma _{2}(x,\rho _{1})\frac{M_{2}(t)}{t}, \end{aligned}$$

(A.47)

$$\begin{aligned}&\frac{1}{t}\ln \frac{ Z(t)}{Z(0)}=\delta c(x,y)\frac{1}{t}\int _{0}^{t}A(s)\mathrm {d}s-\kappa (y,\rho _{2})+ \sigma _{3}(y,\rho _{2})\frac{M_{3}(t)}{t}, \end{aligned}$$

(A.48)

where $M_{i}(t)=\int _{0}^{t}\mathrm {d}B_{i}(t)$, $i=1,2,3$. Let $V_{1}=P+Q(x)A+q(y)Z$, then

$$\begin{aligned} \begin{aligned} \mathrm {d}V_{1} =&\left[ eP_\mathrm{in}-eP-s_{1}(x)Q(x)A- s_{2}(y)q(y)Z\right] \mathrm {d}t+\sigma _{1}P\mathrm {d}B_{1}(t) \\&+\sigma _{2}(x,\rho _{1})Q(x)A\mathrm {d}B_{2}(t)+\sigma _{3}(y,\rho _{2})q(y)Z\mathrm {d}B_{3}(t)\\ \le&\big (eP_\mathrm{in}-h_{1}V_{1}\big )\mathrm {d}t+\sigma _{1}P\mathrm {d}B_{1}(t)+\sigma _{2}(x,\rho _{1})Q(x)A\mathrm {d}B_{2}(t) +\sigma _{3}(y,\rho _{2})q(y)Z\mathrm {d}B_{3}(t), \end{aligned} \end{aligned}$$

(A.49)

where $h_{1}=\min \{e,s_{1}(x),s_{2}(y)\}$. Integrating both sides of Eq. (A.49) from 0 to t and dividing t on both sides lead to

$$\begin{aligned} \frac{\psi _{1}}{t}=eP_\mathrm{in}-e\frac{1}{t}\int _{0}^{t}P(s)\mathrm {d}s -s_{1}(x)Q(x)\frac{1}{t}\int _{0}^{t}A(s)\mathrm {d}s -s_{2}(y)q(y)\frac{1}{t}\int _{0}^{t}Z(s)\mathrm {d}s,\nonumber \\ \end{aligned}$$

(A.50)

where $\psi (t)=V_{1}(t)-V_{1}(0)-\sigma _{1}N_{1}(t)-\sigma _{2}(x,\rho _{1})Q(x)N_{2}(t) -\sigma _{3}(y,\rho _{2})q(y)N_{3}(t),$ $N_{j}(t)=\int _{0}^{t}\varPhi _{j}(s)\mathrm {d}B_{j}(t)$, $j=1,2,3$. According to Theorem 1 and the strong law of large numbers for local martingales (Mao 2006), we have

$$\begin{aligned} \lim _{t\rightarrow \infty }\frac{M_{i}(t)}{t}=0,\ \ \lim _{t\rightarrow \infty }\frac{N_{j}(t)}{t}=0,\ \ \lim _{t\rightarrow \infty }\frac{\psi _{1}(t)}{t}=0,\ \ i,j=1,2,3,\ \ a.s.\nonumber \\ \end{aligned}$$

(A.51)

We next are going to prove Theorem 2 step by step based on Eqs. (A.47)–(A.51).

• (1) We first prove the first conclusion of Theorem 3. From (A.50),

  $$\begin{aligned} \frac{1}{t}\int _{0}^{t}P(s)\mathrm {d}s\le P_\mathrm{in}-\frac{\psi _{1}}{et}, \end{aligned}$$

  (A.52)

  Substituting (A.52) into (A.47) yields

  $$\begin{aligned} \frac{1}{t}\ln \frac{ A(t)}{A(0)}\le \alpha \mu (x) P_\mathrm{in}-\eta (x,\rho _{1}) -\frac{\alpha \mu (x)}{e}\frac{\psi _{1}}{t}+\sigma _{2}(x,\rho _{1})\frac{M_{2}(t)}{t}, \end{aligned}$$

  then when $P_\mathrm{in}<\frac{\eta (x,\rho _{1})}{\alpha \mu (x)}$, together with (A.51) we have $\limsup _{t\rightarrow \infty }\frac{\ln A(t)}{t}<0$, that is

  $$\begin{aligned} \lim _{t\rightarrow \infty }A(t)=0,\ \ a.s. \end{aligned}$$

  Consequently from (A.48), $\lim _{t\rightarrow \infty }Z(t)=0,\ a.s.$ Moreover, according to (A.50) and (A.51), $\langle P \rangle =P_\mathrm{in}$. The first conclusion of Theorem 2 is proved.

• (2) We now proceed to the proof of the second conclusion. From (A.50),

  $$\begin{aligned}&\frac{1}{t}\int _{0}^{t}P(s)\mathrm {d}s= P_\mathrm{in}-\frac{s_{1}(x)Q(x)}{e}\frac{1}{t}\int _{0}^{t}A(s)\mathrm {d}s\nonumber \\&\quad -\frac{s_{2}(y)q(y)}{e}\frac{1}{t}\int _{0}^{t}Z(s)\mathrm {d}s-\frac{\psi _{1}}{et}, \end{aligned}$$

  (A.53)

  Substituting (A.53) into (A.47) yields

  $$\begin{aligned} \frac{1}{t}\ln \frac{ A(t)}{A(0)}=&\alpha \mu (x) P_\mathrm{in}-\eta (x,\rho _{1})-\frac{\alpha \mu (x) s_{1}(x)Q(x)}{e}\frac{1}{t}\int _{0}^{t}A(s)\mathrm {d}s +\sigma _{2}(x,\rho _{1})\frac{M_{2}(t)}{t}\nonumber \\&-\frac{\alpha \mu (x)}{e}\frac{\psi _{1}}{t}-\left( \frac{\alpha \mu (x)s_{2}(y)q(y)}{e}+c(x,y)\right) \frac{1}{t}\int _{0}^{t}Z(s)\mathrm {d}s. \end{aligned}$$

  (A.54)

  Obviously,

  $$\begin{aligned}&\frac{1}{t}\ln \frac{ A(t)}{A(0)}\le \alpha \mu (x) P_\mathrm{in}-\eta (x,\rho _{1})-\frac{\alpha \mu (x) s_{1}(x)Q(x)}{e}\frac{1}{t}\int _{0}^{t}A(s)\mathrm {d}s\\&\quad -\frac{\alpha \mu (x)}{e}\frac{\psi _{1}}{t}+\sigma _{2}(x,\rho _{1})\frac{M_{2}(t)}{t}. \end{aligned}$$

  When $P_\mathrm{in}>\frac{\eta (x,\rho _{1})}{\alpha \mu (x)}$, it then follows from Lemma 4 in Liu and Bai (2016) that

  $$\begin{aligned} \overline{\langle A(x,y) \rangle }\le \frac{(\alpha \mu (x) P_\mathrm{in}-\eta (x,\rho _{1}))e}{\alpha \mu (x) s_{1}(x)Q(x)}=\frac{1}{s_{1}(x)Q(x)}\left( eP_\mathrm{in}-\frac{e\eta (x,\rho _{1})}{\alpha \mu (x)}\right) ,\ \ a.s.\nonumber \\ \end{aligned}$$

  (A.55)

  Substituting (A.55) into (A.48) leads to

  $$\begin{aligned} \limsup _{t\rightarrow \infty }\frac{\ln Z(t)}{t}\le&\delta c(x,y)\overline{\langle A(x,y) \rangle }-\kappa (y,\rho _{2})\\ =&\frac{\delta c(x,y)}{s_{1}(x)Q(x)}\lambda (x,y), \end{aligned}$$

  thus when $\lambda (x,y)<0$, $\limsup _{t\rightarrow \infty }\frac{\ln Z(t)}{t}<0$, i.e., $\lim _{t\rightarrow \infty }Z(t)=0$, a.s. Thus, for any $\varepsilon _{1}>0$, there exists $T_{1} >0$ and a set $\varOmega _{\varepsilon _{1}}$ with ${\mathcal {P}}(\varOmega _{\varepsilon _{1}})>1-\varepsilon _{1}$, such that for any $t>T_{1}$ and $\omega \in \varOmega _{\varepsilon _{1}}$, $\left( \frac{\alpha \mu (x)s_{2}(y)q(y)}{e}+c(x,y)\right) \frac{1}{t}\int _{0}^{t}Z(s)\mathrm {d}s< \varepsilon _{1}$ holds. It then follows from (A.54) that

  $$\begin{aligned} \frac{1}{t}\ln \frac{ A(t)}{A(0)}\ge & {} \alpha \mu (x) P_\mathrm{in}-\eta (x,\rho _{1})-\varepsilon _{1}-\frac{\alpha \mu (x) s_{1}(x)Q(x)}{e}\frac{1}{t}\int _{0}^{t}A(s)\mathrm {d}s\\&-\frac{\alpha \mu (x)}{e}\frac{\psi _{1}}{t}+\sigma _{2}(x,\rho _{1})\frac{M_{2}(t)}{t}. \end{aligned}$$

  Again, according to Lemma 4 in Liu and Bai (2016) and the arbitrariness of $\varepsilon _{1}$,

  $$\begin{aligned} \underline{\langle A(x,y) \rangle }\ge \frac{1}{s_{1}(x)Q(x)}\left( eP_\mathrm{in}-\frac{e\eta (x,\rho _{1})}{\alpha \mu (x)}\right) ,\ \ a.s. \end{aligned}$$

  (A.56)

  Combining (A.56) with (A.55), we have $\langle A(x) \rangle =\frac{1}{s_{1}(x)Q(x)}\left( eP_\mathrm{in}-\frac{e\eta (x,\rho _{1})}{\alpha \mu (x)}\right) ,$ a.s. Consequently, it is easy to obtain from (A.50) that $\langle P(x) \rangle =\frac{\eta (x,\rho _{1})}{\alpha \mu (x)}$, a.s. The second conclusion is thus proved.

• (3) We now are going to prove the last conclusion of Theorem 2. From (A.48) and together with Lemma 1,

  $$\begin{aligned} \overline{\langle A(x,y) \rangle }\le \frac{ \kappa (y,\rho _{2})}{\delta c(x,y)},\ \ a.s. \end{aligned}$$

  (A.57)

  Moreover, according to (A.54) we have

  $$\begin{aligned}&\left( \frac{\alpha \mu (x)s_{2}(y)q(y)}{e}+c(x,y)\right) \underline{\langle Z(x,y) \rangle }\nonumber \\&\ge \alpha \mu (x) P_\mathrm{in}-\eta (x,\rho _{1})-\frac{\alpha \mu (x) s_{1}(x)Q(x)}{e}\overline{\langle A(x,y) \rangle }. \end{aligned}$$

  (A.58)

  Together with (A.57), we can conclude that when $\lambda (x,y)>0$,

  $$\begin{aligned} \underline{\langle Z(x,y) \rangle }\ge \frac{\lambda (x,y)}{s_{2}(y)q(y)+ec(x,y)/\alpha \mu (x)}>0, \ \ a.s. \end{aligned}$$

  Then, we have $\limsup _{t\rightarrow \infty }\frac{\ln A(t)}{t}=0$, a.s., otherwise, $\lim _{t\rightarrow \infty }A(t)=0$, and consequently $\lim _{t\rightarrow \infty }Z(t)=0$, a.s. According to (A.47),

  $$\begin{aligned} \frac{1}{t}\ln \frac{ A(t)}{A(0)}\le \alpha \mu (x)\frac{1}{t}\int _{0}^{t}P(s)\mathrm {d}s-c(x,y)\frac{1}{t}\int _{0}^{t}Z(s)\mathrm {d}s +\sigma _{2}(x,\rho _{1})\frac{M_{2}(t)}{t}, \end{aligned}$$

  consequently,

  $$\begin{aligned} \underline{\langle P(x,y) \rangle }\ge \frac{c(x,y)}{\alpha \mu (x)}\overline{\langle Z(x,y)\rangle }>0,\ \ a.s. \end{aligned}$$

  Applying the same method to (A.48), we have $\underline{\langle A(x,y) \rangle }>0,\ \ a.s$. The proof is thus completed.

$\square$

Appendix B: Proof of Theorem 3

Before the proof of Theorem 3, we need the following lemma from Khasminskii (2012). Suppose that X(t) is a homogeneous Markov process in n-dimension Euclidean space $R^{n}$, satisfying the following stochastic differential equation:

$$\begin{aligned} \mathrm {d}X(t)=b(X)\mathrm {d}t+\sum _{r=1}^{k}\sigma _{r}(X)\mathrm {d}B_{r}(t), \end{aligned}$$

(B.59)

where $\sigma _{r}(X)=(\sigma _{r}^{1}(X),\sigma _{r}^{2}(X)\ldots , \sigma _{r}^{n}(X))^{T}$, $B(X)=(a_{ij}(X))_{n\times n}$ is the diffusion matrix of X(t) with $a_{ij}(X)=\sum _{r=1}^{k}\sigma _{r}^{i}(X)\sigma _{r}^{j}(X)$.

Lemma 2

(Khasminskii 2012). If there exists a bounded open domain $U\subset R^{n}$ with regular boundary, satisfying the following properties:

(H1):

The diffusion matrix B(x) is strictly positive definite for all $x\in U$;

(H2):

There exists a nonnegative $C^{2}$-function V(X) such that ${\mathscr {L}}V(X)$ is negative on $X\in R^{n}\setminus U$.

Then, the Markov process X(t) of the stochastic model (B.59) admits a unique stationary distribution $\pi (\cdot )$, and for any integrable function $f(\cdot )$ with regard to the measure $\pi$, the following equation holds,

$$\begin{aligned} {\mathcal {P}}\left( \lim _{t\rightarrow \infty }\frac{1}{t}\int _{0}^{t}f(X(t))\mathrm {d}t=\int _{R^{n}}f(x)\pi (\mathrm {d}x)\right) =1. \end{aligned}$$

In what follows, we are going to apply Lemma 2 to prove the existence of a unique ergodic stationary distribution for model (2).

Proof

Define a nonnegative ${\mathcal {C}}^{2}-$Lyapunov function

$$\begin{aligned} V_{2}=M_{2}V_{21}+V_{22}+V_{23}-V_{2}^{*}, \end{aligned}$$

$V_{2}^{*}$ is the minimum value of $V_{2}(P,A,Z)$ and

$$\begin{aligned} V_{21}=-m_{3}\ln Z-m_{4}\ln A+m_{5} Z-V_{1},~~~ V_{22}=\frac{1}{\theta _{1}+2}V_{1}^{\theta _{1}+2},~~~ V_{23}=-\ln P, \end{aligned}$$

where $m_{3}=\frac{s_{1}(x)Q(x)}{\delta c(x,y)},\ m_{4}=\frac{e}{ \alpha \mu (x) },\ m_{5}=(m_{4}c(x,y)+s_{2}(y)q(y))/(s_{2}(y)+d_{2})$, $0<\theta _{1}<1$ is a constant which is small enough such that

$$\begin{aligned} h_{1}(x,y):=\min \{e,s_{1}(x),s_{2}(y)\}-\frac{\theta _{1}+1}{2}\max \{\sigma _{1}^{2}, \sigma _{2}^{2}(x,\rho _{1}),\sigma _{3}^{2}(y,\rho _{2})\}>0, \end{aligned}$$

and the positive constant $M_{2}$ satisfying

$$\begin{aligned} -M_{2}\lambda (x,y)+f_{1}^{u}\le -2, \end{aligned}$$

where $f^{u}=\sup _{t\in (0,\infty )}f(t)$ and the function $f_{1}$ will be determined later.

By using the It$\hat{\mathrm {o}}$ formula, we have

$$\begin{aligned} \begin{aligned}&{\mathscr {L}}V_{21}\\&\quad =-s_{1}(x)Q(x)A+m_{3}\kappa (y,\rho _{2})-e P+m_{4}c(x,y)Z+m_{4}\eta (x,\rho _{1})\\&\qquad +m_{5}\delta c(x,y)AZ-m_{5}(s_{2}(y)+d_{2})Z -eP_\mathrm{in}+eP+s_{1}(x)Q(x)A+s_{2}(y)q(y)Z\\&\quad =m_{3}\kappa (y,\rho _{2})+m_{4}\eta (x,\rho _{1})\\&\qquad +m_{5}\delta c(x,y)AZ-eP_\mathrm{in}\\&\quad =-\lambda (x,y)+m_{5}\delta c(x,y)AZ.\\&{\mathscr {L}}V_{22}\\&\quad = eP_\mathrm{in}\big (P+Q(x)A+q(y)Z\big )^{\theta _{1}+1}\\&\qquad -\big (P+Q(x)A+q(y)Z\big )^{\theta _{1}+1}(eP+s_{1}(x)Q(x)A+s_{2}(y)q(y)Z)\\&\qquad +\frac{\theta _{1}+1}{2}\big (P+Q(x)A+q(y)Z\big )^{\theta _{1}}\big (\sigma _{1}^{2}P^{2}\\&\qquad +\sigma _{2}^{2}(x,\rho _{1})Q^{2}(x)A^{2}+\sigma _{3}^{2}(y,\rho _{2})q^{2}(y)Z^{2}\big )\\&\quad \le eP_\mathrm{in}\big (P+Q(x)A+q(y)Z\big )^{\theta _{1}+1}\\&\qquad -\min \{e,s_{1}(x),s_{2}(y)\}\big (P+Q(x)A+q(y)Z\big )^{\theta _{1}+2}\\&\qquad +\frac{(\theta +1)\max \{\sigma _{1}^{2},\sigma _{2}^{2}(x,\rho _{1}),\sigma _{3}^{2}(y,\rho _{2})\}}{2}\\&\qquad \big (P+Q(x)A+q(y)Z\big )^{\theta _{1}+2}\\&\quad =-h_{1}(x,y)(P+Q(x)A+q(y)Z\big )^{\theta _{1}+2}\\&\qquad +eP_\mathrm{in}\big (P+Q(x)A+q(y)Z\big )^{\theta _{1}+1}.\\&{\mathscr {L}}V_{23}\\&\quad =-\frac{eP_\mathrm{in}}{P}+e-\frac{d_{1}Q(x)A+ \big (Q(x)-\delta q(y)\big )c(x,y)AZ +d_{2}q(y) Z}{P}\\&\qquad +\alpha \mu (x)Q(x)A+\frac{\sigma _{1}^{2}}{2}\\&\quad \le -\frac{eP_\mathrm{in}}{P}+\alpha \mu (x)Q(x)A+e+\frac{\sigma _{1}^{2}}{2}. \end{aligned} \end{aligned}$$

Then,

$$\begin{aligned}&{\mathcal {L}}V_{2}\\&\quad \le -M_{2}\lambda (x,y)+M_{2}m_{5}\delta c(x,y)AZ-h_{1}(x,y)(P+Q(x)A+q(y)Z\big )^{\theta _{1}+2}\\&\qquad +eP_\mathrm{in}\big (P+Q(x)A+q(y)Z\big )^{\theta _{1}+1}\\&\qquad -\frac{eP_\mathrm{in}}{P}+\alpha \mu (x)Q(x)A+e+\frac{\sigma _{1}^{2}}{2}\\&\quad := -M_{2}\lambda (x,y)+M_{2}m_{5}\delta c(x,y)AZ-\frac{eP_\mathrm{in}}{P}+f_{1}(P,A,Z), \end{aligned}$$

where

$$\begin{aligned} f_{1}(P,A,Z)=&-h_{1}(x,y)(P+Q(x)A+q(y)Z\big )^{\theta _{1}+2}\\&+eP_\mathrm{in}\big (P+Q(x)A+q(y)Z\big )^{\theta _{1}+1} +\alpha \mu (x)Q(x)A+e+\frac{\sigma _{1}^{2}}{2}. \end{aligned}$$

Denote

$$\begin{aligned} \varUpsilon (N,{\tilde{P}},Z)=-M_{2}\lambda (x,y)+M_{2}m_{5}\delta c(x,y)AZ-\frac{eP_\mathrm{in}}{P}+f_{1}(P,A,Z), \end{aligned}$$

then

$$\begin{aligned} \varUpsilon (P,A,Z)\le {\left\{ \begin{array}{ll} \begin{aligned} \varUpsilon (0, A, Z)\rightarrow -\infty ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ &{}\text {as}\ \ P\rightarrow 0^{+},\\ \varUpsilon (\infty , {\tilde{P}}, Z)\rightarrow -\infty ,\ \ \ \ \ \ \ \ \ \ \ \ \ &{}\text {as}\ \ P\rightarrow \infty ,\\ -M_{2}\lambda (x,y)+f_{1}^{u}\le -1,\ \ \ \ \ \ \ \ &{}\text {as}\ \ A\rightarrow 0^{+},\\ \varUpsilon (P, \infty , Z)\rightarrow -\infty ,\ \ \ \ \ \ \ \ \ \ \ \ \ &{}\text {as}\ \ A\rightarrow \infty ,\\ -M_{2}\lambda (x,y)+f_{1}^{u}\le -1,\ \ \ \ \ \ \ \ &{}\text {as}\ \ Z\rightarrow 0^{+},\\ \varUpsilon (P, A, \infty )\rightarrow -\infty ,\ \ \ \ \ \ \ \ \ \ \ \ \ &{}\text {as}\ \ Z\rightarrow \infty . \end{aligned} \end{array}\right. } \end{aligned}$$

Thus, we can take $0<\varepsilon _{2} <1$ sufficiently small such that

$$\begin{aligned} {\mathscr {L}} V_{2}(t, P, A, Z)\le -1,\ \ \text {for}\ \text {any}\ \ (P,A,Z)\in R_{+}^{3}\setminus U_{1}, \end{aligned}$$

where $U_{1}=[\varepsilon _{2}, \frac{1}{\varepsilon _{2}}]\times [\varepsilon _{2}, \frac{1}{\varepsilon _{2}}]\times [\varepsilon _{2}, \frac{1}{\varepsilon _{2}}]$. Moreover, there exists a positive constant $M_{3}=\min \{\sigma _{1}^{2}P^{2},\sigma _{2}^{2}(x,\rho _{1})A^{2},\sigma _{3}^{2}(y,\rho _{2})Z^{2}\}$, such that

$$\begin{aligned} \sigma _{1}^{2}P^{2}\zeta _{1}^{2}+\sigma _{2}^{2}(x,\rho _{1})A^{2}\zeta _{2}^{2} +\sigma _{3}^{2}(y,\rho _{2})Z^{2}\zeta _{3}^{2}\ge M_{3}\Vert \zeta \Vert ^{2}, \end{aligned}$$

for all $(P, A, Z)\in U_{1}$, $\zeta =(\zeta _{1},\zeta _{2},\zeta _{3})\in R^{3}$. Then, based on Lemma 2, model (2) has a unique ergodic stationary distribution $\pi (.)$. Moreover, according to the ergodic property, the solution of model (2) $\varPhi (t)=(P(t),A(t),Z(t))$ satisfying

$$\begin{aligned} {\mathcal {P}}\left\{ \langle \varPhi _{i}(x,y)\rangle =\int _{R_{+}^{3}}\phi _{i}\pi (\mathrm {d}\phi _{1},\mathrm {d}\phi _{2},\mathrm {d}\phi _{3})\right\} =1,\ \ i=1,2,3. \end{aligned}$$

(B.60)

We then show that $\lim _{t\rightarrow \infty }\frac{\ln A(t)}{t}=0$, from (A.47), (B.60) and the strong law of large numbers for martingales, the limits of $\frac{1}{t}\ln A(t)$ exist. If $\lim _{t\rightarrow \infty }\frac{\ln A(t)}{t}\ne 0$, then together with Lemma 1, $\limsup _{t\rightarrow \infty }\frac{\ln A(t)}{t}<0$, consequently, $\lim _{t\rightarrow \infty }A(t)=0$, which contradicts with (B.60). Similarly, we have $\lim _{t\rightarrow \infty }\frac{\ln Z(t)}{t}=0$. Then, according to the It$\hat{\mathrm {o}}$ formula we have

$$\begin{aligned} {\left\{ \begin{array}{ll} \begin{aligned} 0=&{}\lim _{t\rightarrow \infty }\frac{\ln A(t)-\ln A(0)}{t}=\alpha \mu (x)\langle P(x,y)\rangle -c(x,y)\langle Z(x,y)\rangle -\eta (x,\rho _{1})\\ &{}+\sigma _{2}(x,\rho _{1})\lim _{t\rightarrow \infty }\frac{M_{2}(t)}{t},\\ 0=&{}\lim _{t\rightarrow \infty }\frac{\ln Z(t)-\ln Z(0)}{t}=\delta c(x,y) \langle A(x,y)\rangle -\kappa (y,\rho _{2})+\sigma _{3}(y,\rho _{2})\lim _{t\rightarrow \infty }\frac{M_{3}(t)}{t},\\ 0=&{}\lim _{t\rightarrow \infty }\frac{V_{11}(t)-V_{11}(0)}{t}= eP_\mathrm{in}-e\langle P(x,y)\rangle -s_{1}(x)Q(x)\langle A(x,y)\rangle - s_{2}(y)q(y)\langle Z(x,y)\\ &{}+\sigma _{1}\lim _{t\rightarrow \infty }\frac{N_{1}(t)}{t}+ \sigma _{2}(x,\rho _{1})Q(x)\lim _{t\rightarrow \infty }\frac{N_{2}(t)}{t} +\sigma _{3}(y,\rho _{2})q(y)\lim _{t\rightarrow \infty }\frac{N_{3}(t)}{t}. \end{aligned} \end{array}\right. } \end{aligned}$$

(B.61)

By using and the strong law of large numbers for martingales, $\lim _{t\rightarrow \infty }\frac{M_{i}(t)}{t}=\lim _{t\rightarrow \infty }\frac{N_{j}(t)}{t}=0$, $i=2,3$, $j=1,2,3$, then obviously (B.61) has a unique positive solution, provided that $\lambda (x,y)>0$, and

$$\begin{aligned} {\left\{ \begin{array}{ll} \begin{aligned} &{}\langle P(x,y)\rangle =\frac{c(x,y)\langle Z(x,y)\rangle +\eta (x,\rho _{1})}{\alpha \mu (x)},\\ &{}\langle A(x,y)\rangle =\frac{\kappa (y,\rho _{2})}{\delta c(x,y)},\\ &{}\langle Z(x,y)\rangle =\frac{\lambda (x,y)}{s_{2}(y)q(y)+\frac{e c(x,y)}{\alpha \mu (x)}}. \end{aligned} \end{array}\right. } \end{aligned}$$

(B.62)

This completes the proof of Theorem 3. $\square$

Appendix C: Proofs of Theorems 5 and 6

1.1 Proof of Theorem 5

Our main aim in this part is to derive the conditions under which $x^{*}$ is evolutionary stable and convergence stable. To achieve that, we first need the following results. From the first equation of (B.61), when $y=h$ is a constant we have

$$\begin{aligned} \alpha \mu (x)\langle P(x,h)\rangle -c(x)\langle Z(x,h)\rangle =\eta (x,\rho _{1}). \end{aligned}$$

(C.63)

Taking the partial derivative of both sides of (C.63) with respect to x leads to

$$\begin{aligned} \alpha \mu '(x)\langle P(x,h)\rangle +\alpha \mu (x)\langle P(x,h)\rangle '-c'(x,h)\langle Z(x,h)\rangle -c(x,h) \langle Z(x,h)\rangle '=\eta '(x,\rho _{1}).\nonumber \\ \end{aligned}$$

(C.64)

Combining with Eq. (25), we have

$$\begin{aligned} \begin{aligned} c(x^{*},h)\langle Z(x^*,h)\rangle '-\alpha \mu (x^{*})\langle P(x^{*},h)\rangle '=0. \end{aligned} \end{aligned}$$

(C.65)

We know that $x^{*}$ is an ESS, i.e., the fitness will achieve its maximum at $x^{*}$, provided that

$$\begin{aligned} \frac{\partial ^{2} F_{2}(x,x_\mathrm{mut})}{\partial x_\mathrm{mut}^{2}}\bigg |_{x_\mathrm{mut}=x=x^{*}}=\alpha \mu ''(x^{*})\langle P(x^{*},h)\rangle -c''(x^{*},h)\langle Z(x^{*},h)\rangle -\eta ''(x^{*},\rho _{1})<0.\nonumber \\ \end{aligned}$$

(C.66)

Taking the second partial derivative of both sides of (C.63) with respect to x leads to

$$\begin{aligned}&\alpha \mu ''(x^{*})\langle P(x^{*},h)\rangle -c''(x^{*},h)\langle Z(x^{*},h)\rangle -\eta ''(x^{*},\rho _{1})\nonumber \\&\quad = 2c'(x^{*},h)\langle Z(x^{*},h)\rangle '-2\alpha \mu '(x^{*})\langle P(x^{*},h)\rangle '+c(x,h)\langle Z(x^{*},h)\rangle ''\nonumber \\&\qquad -\alpha \mu (x)\langle P(x^{*},h)\rangle ''. \end{aligned}$$

(C.67)

Substituting (C.65) into (C.67), we have

$$\begin{aligned} \begin{aligned}&\frac{\partial ^{2} F_{2}(x,x_\mathrm{mut})}{\partial x_\mathrm{mut}^{2}}\bigg |_{x_\mathrm{mut}=x=x^{*}}\\&\quad =2c(x^{*},h)\langle Z(x^{*},h)\rangle '\left( \frac{c'(x^{*},h)}{c(x^{*},h)}- \frac{\mu '(x^{*})}{\mu (x^{*})}\right) -c(x^{*},h)\langle Z(x^{*},h)\rangle '\\&\qquad \left( \frac{\langle P(x^{*},h)\rangle ''}{\langle P(x^{*},h)\rangle '}-\frac{\langle Z(x^{*},h)\rangle ''}{\langle Z(x^{*},h)\rangle '}\right) \\&\quad =c(x^{*},h)\langle Z(x^{*},h)\rangle '(2h_{1}(x^{*})-h_{2}(x^{*})), \end{aligned} \end{aligned}$$

(C.68)

where

$$\begin{aligned} \begin{aligned} h_{1}(x)=\frac{c'(x,h)}{c(x,h)}- \frac{\mu '(x)}{\mu (x)}, \ \ \ \ \ h_{2}(x)=\frac{\langle P(x,h)\rangle ''}{\langle P(x,h)\rangle '}-\frac{\langle Z(x,h)\rangle ''}{\langle Z(x,h)\rangle '}. \end{aligned} \end{aligned}$$

(C.69)

Obviously, if $\langle Z(x^{*},h)\rangle '(2h_{1}(x^{*})-h_{2}(x^{*}))<0$ holds, then $x^{*}$ is an ESS. Once $x^{*}$ is an ESS, there is no possibility of further evolution changes by small mutation. Moreover, the singular strategy $x^{*}$ is convergence stable, provided that

$$\begin{aligned} \begin{aligned}&\frac{\partial D_{2}(x)}{\partial x}\bigg |_{x=x^{*}}\\&\quad =\alpha \mu ''(x^{*})\langle P(x^{*},h)\rangle -c''(x^{*},h)\langle Z(x^{*},h)\rangle \\&\qquad -\eta ''(x^{*},\rho _{1})+\alpha \mu '(x^{*})\langle P(x^{*},h)\rangle ' -c'(x^{*},h)\langle Z(x^{*},h)\rangle '\\&\quad =c'(x^{*},h)\langle Z(x^{*},h)\rangle '-\alpha \mu '(x^{*})\langle P(x^{*},h)\rangle '+c(x^{*},h)\langle Z(x^{*},h)\rangle ''-\alpha \mu (x^{*})\langle P(x^{*},h)\rangle ''\\&\quad = c(x^{*},h)\langle Z(x^{*},h)\rangle '(h_{1}(x^{*})-h_{2}(x^{*}))\\&\quad <0. \end{aligned} \end{aligned}$$

(C.70)

We then determine the sign of $\langle Z(x^{*},h)\rangle '$. From the third equation of (B.61) and the second equation of (B.62),

$$\begin{aligned} eP_\mathrm{in}-e\langle P(x,h)\rangle -s_{2}(h)q(h)\langle Z(x,h)\rangle =\frac{s_{1}(x)Q(x)\kappa (h,\rho _{2})}{\delta c(x,h)}, \end{aligned}$$

then taking the partial derivative of both sides with respect to x, we have

$$\begin{aligned} -e \langle P(x,h)\rangle '- s_{2}(h)q(h)\langle Z(x,h)\rangle '=\frac{\kappa (h,\rho _{2})}{\delta }g'(x), \end{aligned}$$

where $g(x)=\frac{s_{1}(x)}{c(x,h)}Q(x).$ Since $\left( \frac{s_{1}(x)}{c(x,h)}\right) '=\frac{2k_{1}x\exp (v(x-\theta h)^{2})}{c_\mathrm{m}}(vx^{2}-v\theta h x+1)> 0$ and $Q'(x)>0$, then $g'(x)>0$. This together with (C.65) we can conclude that $\langle Z(x^{*},h)\rangle '<0$ and $\langle P(x^{*},h)\rangle '<0$. Then, from (C.68) and (C.70) we know that when $h_{1}(x^{*})<h_{2}(x^{*})$, $x^{*}$ is a repellor. When $h_{1}(x^{*})>h_{2}(x^{*})$, $x^{*}$ is convergence stable; meanwhile, it is an ESS as long as $2h_{1}(x^{*})>h_{2}(x^{*})$ or an evolutionary branching point if $2h_{1}(x^{*})<h_{2}(x^{*})$. Theorem 5 is thus proved.

We then give some analysis about Table 2. Combining with Eqs. (C.63), (25) and (B.60) leads to

$$\begin{aligned} \left( \frac{\mu '(x^{*})}{\mu (x^{*})}-\frac{c'(x^{*},h)}{c(x^{*},h)}\right) \left( \frac{ \eta '(x^{*},\rho _{1})}{\eta (x^{*},\rho _{1})}-\frac{\mu '(x^{*})}{\mu (x^{*})}\right) >0, \end{aligned}$$

and

$$\begin{aligned} \left( \frac{\mu '(x^{*})}{\mu (x^{*})}-\frac{c'(x^{*},h)}{c(x^{*},h)}\right) \left( \frac{ \eta '(x^{*},\rho _{1})}{\eta (x^{*},\rho _{1})}-\frac{c'(x^{*},h)}{c(x^{*},h)}\right) >0, \end{aligned}$$

which implies that

$$\begin{aligned} \frac{\eta '(x^{*},\rho _{1})}{\eta (x^{*},\rho _{1})}> \frac{\mu '(x^{*})}{\mu (x^{*})}> \frac{c'(x^{*},h)}{c(x^{*},h)}\ \ \ \ \mathrm {or}\ \ \ \ \ \frac{c'(x^{*},h)}{c(x^{*},h)}> \frac{\mu '(x^{*})}{\mu (x^{*})}> \frac{ \eta '(x^{*},\rho _{1})}{\eta (x^{*},\rho _{1})}>0.\nonumber \\ \end{aligned}$$

(C.71)

It easy to see that $h_{1}(x^{*})\ne 0$, $x^{*}\in \left( 0,\ \min \left\{ \theta h, \sqrt{\frac{a_{3}}{a_{1}}}\right\} \right)$ or $x^{*}\in \left( \theta h,\ \infty \right)$.

From (C.69), the sign of $h_{1}(x^{*})$ is only determined by the sign of $\frac{c'(x^{*},h)}{c(x^{*},h)}-\frac{\mu '(x^{*})}{\mu (x^{*})}$. When $\frac{c'(x^{*},h)}{c(x^{*},h)}>\frac{\mu '(x^{*})}{\mu (x^{*})}>0$, $h_{1}(x^{*})>0$, satisfying $2h_{1}(x^{*})>h_{1}(x^{*})$. Then once $x^{*}$ is convergence stable, i.e., $h_{1}(x^{*})>h_{2}(x^{*})$, the inequality $2h_{1}(x^{*})>h_{2}(x^{*})$ is also holds, i.e., $x^{*}$ is also an ESS. At this situation, $x^{*}$ is either a CSS or a repellor, and $x^{*}\in (0, \min \{\theta h, \sqrt{a_{3}/a_{1}} \})$. When $\frac{\mu '(x^{*})}{\mu (x^{*})}>\frac{c'(x^{*},h)}{c(x^{*},h)}$, $h_{1}(x^{*})<0$, and when $h_{2}(x^{*})\ge 0$, $h_{1}(x^{*})< h_{2}(x^{*})$ always holds, thus $x^{*}$ is a repellor. When $h_{2}(x^{*})<0$, it can be divided into three scenarios: (1) $h_{1}(x^{*})< h_{2}(x^{*})<0$, at this situation, $x^{*}$ is a repellor; (2) $h_{2}(x^{*})< 2h_{1}(x^{*})<0$, at this point, $x^{*}$ is a CSS; (3) $h_{1}(x^{*})> h_{2}(x^{*})>2h_{1}(x^{*})$, then $x^{*}$ is an evolutionary branching point, which implies the algal population will split up into two species with different cell size. This confirms Table 2.

1.2 Proof of Theorem 6

In this part, we mainly investigate how environmental fluctuation and phosphorus content affect the trend of algal evolution in the presence of zooplankton. According to (25),

$$\begin{aligned} \begin{aligned}&\frac{\partial D_{2}(x^{*})}{\partial \rho _{1}}\\&\quad =\alpha \mu '(x^{*})\frac{\partial \langle P(x^{*},h)\rangle }{\partial \rho _{1}}-c'(x^{*},h)\frac{\partial \langle Z(x^{*},h)\rangle }{\partial \rho _{1}}-\frac{\partial ^{2}\eta (x,\rho _{1})}{\partial x\partial \rho _{1}}\bigg |_{x=x^{*}}\\&\quad =\alpha \frac{\mu '(x^{*})}{\mu (x^{*})}\mu (x^{*})\frac{\partial \langle P(x^{*},h)\rangle }{\partial \rho _{1}}-\frac{c'(x^{*},h)}{c(x^{*},h)}c(x^{*},h)\frac{\partial \langle Z(x^{*},h)\rangle }{\partial \rho _{1}}-\frac{\partial ^{2}\eta (x,\rho _{1})}{\partial x\partial \rho _{1}}\bigg |_{x=x^{*}}. \end{aligned} \end{aligned}$$

(C.72)

From (21),

$$\begin{aligned} \frac{\mathrm {\partial }\langle Z(x,h)\rangle }{\mathrm {\partial }\rho _{1}}=-\frac{\frac{e}{\alpha \mu (x)}\sigma _{2}(x,\rho _{1})\frac{\partial \sigma _{2}(x,\rho _{1})}{\partial \rho _{1}}}{s_{2}(h)q(h)+\frac{ec(x,h)}{\alpha \mu (x)}}<0, \end{aligned}$$

and

$$\begin{aligned} \frac{\mathrm {\partial }\langle P(x,h)\rangle }{\mathrm {\partial }\rho _{1}}=\frac{\sigma _{2}(x,\rho _{1})\frac{\partial \sigma _{2}(x,\rho _{1})}{\partial \rho _{1}}}{\alpha \mu (x)}\left( -\frac{\frac{ec(x,h)}{\alpha \mu (x)}}{s_{2}(h)q(h)+\frac{ec(x,h)}{\alpha \mu (x)}}+1\right) >0. \end{aligned}$$

Moreover, from the first equation of (B.61),

$$\begin{aligned} \alpha \mu (x)\langle P(x,h)\rangle -c(x,h)\langle Z(x,h)\rangle =\eta (x,\rho _{1}), \end{aligned}$$

then

$$\begin{aligned} \alpha \mu (x)\frac{\partial \langle P(x,h)\rangle }{\partial \rho _{1}}-c(x,h)\frac{\partial \langle Z(x,h)\rangle }{\partial \rho _{1}}=\frac{\partial \eta (x,\rho _{1})}{\partial \rho _{1}}. \end{aligned}$$

It then follows from (C.72), when $\frac{\mu '(x^{*})}{\mu (x^{*})} >\frac{c'(x^{*},h)}{c(x^{*},h)}$,

$$\begin{aligned} \begin{aligned} \frac{\partial D_{2}(x^{*})}{\partial \rho _{1}}<&\frac{\mu '(x^{*})}{\mu (x^{*})}\left[ \alpha \mu (x^{*})\frac{\partial \langle P(x^{*},h)\rangle }{\partial \rho _{1}}-c(x^{*},h)\frac{\partial \langle Z(x^{*},h)\rangle }{\partial \rho _{1}}\right] -\frac{\partial ^{2}\eta (x,\rho _{1})}{\partial x\partial \rho _{1}}\bigg |_{x=x^{*}}\\ =&\frac{\mu '(x^{*})}{\mu (x^{*})}\frac{\partial \eta (x^{*},\rho _{1})}{\partial \rho _{1}}-\frac{\partial ^{2}\eta (x,\rho _{1})}{\partial x\partial \rho _{1}}\bigg |_{x=x^{*}}\\ =&\frac{a_{3}-a_{1}x^{*2}}{(a_{1}x^{*2}+a_{2}x^{*}+a_{3})x^{*}} \frac{x^{*4}k_{3}^{2}k_{4}\rho _{1}}{(k_{4}+\rho _{1})^{3}} -\frac{4x^{*3}k_{3}^{2}k_{4}\rho _{1}}{(k_{4}+\rho _{1})^{3}}\\ =&\left( \frac{-5a_{1}x^{*2}-4a_{2}x^{*}-3a_{3}}{a_{1}x^{*2}+a_{2}x^{*}+a_{3}}\right) \frac{x^{*3}k_{3}^{2}k_{4}\rho _{1}}{(k_{4}+\rho _{1})^{3}}\\ <&0. \end{aligned} \end{aligned}$$

Otherwise, $\frac{c'(x^{*},h)}{c(x^{*},h)}>\frac{\mu '(x^{*})}{\mu (x^{*})}>0$, since $v< \frac{4}{\theta ^{2}h^{2}}$, then

$$\begin{aligned} \begin{aligned} \frac{\partial D_{2}(x^{*})}{\partial \rho _{1}}<&\frac{c'(x^{*},h)}{c(x^{*},h)}\frac{\partial \eta (x^{*},\rho _{1})}{\partial \rho _{1}}-\frac{\partial ^{2}\eta (x,\rho _{1})}{\partial x\partial \rho _{1}}\bigg |_{x=x^{*}}\\ =&2v(\theta h-x^{*})\frac{x^{*4}k_{3}^{2}k_{4}\rho _{1}}{(k_{4}+\rho _{1})^{3}} -\frac{4x^{*3}k_{3}^{2}k_{4}\rho _{1}}{(k_{4}+\rho _{1})^{3}}\\ =&\left( -2vx^{*2}+2v\theta h x^{*}-4\right) \frac{x^{*3}k_{3}^{2}k_{4}\rho _{1}}{(k_{4}+\rho _{1})^{3}}\\ <&0. \end{aligned} \end{aligned}$$

Since our concern is the CSS or evolutionary branching point, then $-\frac{\partial D_{1}(x^{*})}{\partial x^{*}}>0.$ By the implicit function theorem, we have

$$\begin{aligned} \frac{\mathrm {d}x^{*}}{\mathrm {d}\rho _{1}}=\frac{\frac{\partial D_{2}(x^{*})}{\partial \rho _{2}}}{-\frac{\partial D_{2}(x^{*})}{\partial x^{*}}}<0. \end{aligned}$$

Moreover, from Eq. (21), we have $\frac{\partial \langle P(x,h)\rangle }{\partial P_\mathrm{in}}=\frac{c(x,h)}{\alpha \mu (x)}\frac{\partial \langle Z(x,h)\rangle }{\partial P_\mathrm{in}}$ and $\frac{\partial \langle Z(x,h)\rangle }{\partial P_\mathrm{in}}=\frac{e}{s_{2}(h)q(h)+e c(x,h)/\alpha \mu (x)}>0$. Then,

$$\begin{aligned} \frac{\partial D_{2}(x^{*})}{\partial P_\mathrm{in}}&=\alpha \mu '(x)\frac{\partial \langle P(x^{*},h)\rangle }{\partial P_\mathrm{in}} -c'(x^{*},h)\frac{\partial \langle Z(x^{*},h)\rangle }{\partial P_\mathrm{in}}\\&=c(x^{*},h)\frac{\partial \langle Z(x^{*},h)\rangle }{\partial P_\mathrm{in}}\left( \frac{\mu '(x^{*})}{\mu (x^{*})}-\frac{c'(x^{*},h)}{c(x^{*},h)}\right) \\&=-c(x^{*},h)h_{1}(x^{*})\frac{\partial \langle Z(x^{*},h)\rangle }{\partial P_\mathrm{in}}. \end{aligned}$$

Then if $h_{1}(x^{*})<0$,

$$\begin{aligned} \frac{\mathrm {d}x^{*}}{\mathrm {d}P_\mathrm{in}}=\frac{\frac{\partial D_{2}(x^{*})}{\partial P_\mathrm{in}}}{-\frac{\partial D_{2}(x^{*})}{\partial x^{*}}}>0, \end{aligned}$$

otherwise $h_{1}(x^{*})>0,$ and $\frac{\mathrm {d}x^{*}}{\mathrm {d}P_\mathrm{in}}<0$. This completes the proof.

Appendix D: Proofs of Theorems 7 and 8

1.1 Proof of Theorem 7

Proof

Since the coefficients of model (26) are locally Lipschitz in $R_{+}^{4}$, then there exists a unique local solution $(P(t),A_{1}(t),A_{2}(t),Z(t))$ of model (26) on the interval $(0,\tau _{e})$, where $\tau _{e}$ denotes the explosion time. In order to prove $\tau _{e}=\infty$, we only need to construct a nonnegative ${\mathcal {C}}^{2}$-function $V_{3}$ satisfying ${\mathcal {L}}V_{3}\le M_{3}$, where $M_{3}$ is a positive constant (Mao et al. 2002). Define the nonnegativity function $V_{3}: R_{+}^{4}\rightarrow R_{+}$ by

$$\begin{aligned} V_{3}=m_{6}V_{31}-V_{32}-V_{3}^{*}, \end{aligned}$$

where $V_{31}=P+\sum _{i=1}^{2}Q(x_{i})A_{i}+qZ$, $V_{32}=m_{7}\ln P+\sum _{i=1}^{2}\ln A_{i}+\ln Z$, $V_{3}^{*}$ is the minimum value point of $V_{3}(P,A_{1},A_{2},Z)$ and $m_{6}=\frac{2\check{c}}{s_{2}q}$, $m_{7}=\frac{\delta {\hat{c}} }{\alpha \check{u} \check{Q}}.$ Then by using the It$\hat{\mathrm {o}}$ formula, we have

$$\begin{aligned} {\mathcal {L}}V_{3}\le m_{6} eP_\mathrm{in}+\sum _{i=1}^{2}\eta (x_{i},\rho _{1})+m_{7}e+\kappa (h,\rho _{2}):=M_{3}. \end{aligned}$$

Arguing similarly as in Mao et al. (2002), we can obtain the first part of Theorem 7. Moreover, by using the It$\hat{\mathrm {o}}$ formula to model (26) again, we have

$$\begin{aligned} \begin{aligned}&\mathrm {d}V_{31}\\&\quad = \left[ eP_\mathrm{in}-eP-\sum _{i=1}^{2}s_{1}(x_{i})Q(x_{i})A_{i}- s_{2}(h)q(h)Z\right] \mathrm {d}t\\&\qquad +\sigma _{1}P\mathrm {d}B_{1}(t)+\sigma _{3}(h,\rho _{2})q(h)Z\mathrm {d}B_{4}(t)\\&\qquad +\sum _{i=1}^{2}\sigma _{2}(x_{i},\rho _{1})Q(x_{i})A_{i}\mathrm {d}B_{i+1}(t)\\&\quad \le \big (eP_\mathrm{in}-h_{2}V_{31}\big )\mathrm {d}t+\sigma _{1}P\mathrm {d}B_{1}(t)\\&+\qquad \sum _{i=1}^{2}\sigma _{2}(x_{i},\rho _{1})Q(x_{i})A_{i}\mathrm {d}B_{i+1}(t) +\sigma _{3}(h,\rho _{2})q(h)Z\mathrm {d}B_{4}(t), \end{aligned} \end{aligned}$$

where $h_{2}=\min \{e,{\hat{s}}_{1},s_{2}(h)\}$. Considering the following model

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \mathrm {d}Y = \big (eP_\mathrm{in}-h_{2}V_{31}\big )\mathrm {d}t+\sigma _{1}P\mathrm {d}B_{1}(t) +\sum _{i=1}^{2}\sigma _{2}(x_{i},\rho _{1})Q(x_{i})A_{i}\mathrm {d}B_{i+1}(t)\\ +\sigma _{3}(h,\rho _{2})q(h)Z\mathrm {d}B_{4}(t),\\ \displaystyle { Y(0)=V_{31}(0).} \end{array} \right. \end{aligned}$$

(D.73)

Then, the solution of Eq. (D.73) has the following form

$$\begin{aligned} Y(t)=\frac{eP_\mathrm{in}}{h_{2}}+Y(0)-\left( \frac{eP_\mathrm{in}}{h_{2}}\right) \exp (-h_{2}t)+M(t), \end{aligned}$$

where

$$\begin{aligned} M(t)=&\sigma _{1}\int _{0}^{t}\exp [-h_{2}(t-s)]P(s)\mathrm {d}B_{1}(s)+\sum _{i=1}^{2}\sigma _{2}(x_{i},\rho _{1})Q(x_{i})\\&\int _{0}^{t}\exp [-h_{2}(t-s)]A_{i}(s)\mathrm {d}B_{i+1}(s)\\&+\sigma _{3}(h,\rho _{2})q(h)\int _{0}^{t}\exp [-h_{2}(t-s)] Z(s)\mathrm {d}B_{4}(s) \end{aligned}$$

is a local martingale satisfying $M(0)=0$, a.s. The stochastic comparison theorem implies $V_{31}(t)\le Y(t)$, a.s.

Define $Y(t)=Y(0)+ B(t)-U(t) + M(t),$ where $B(t) = \frac{eP_\mathrm{in}}{h_{2}} [1-\exp (-h_{2}t)],\ U(t) = Y(0)[1- \exp (-h_{2}t)],\ B(0)=U(0)=0$. Obviously, B(t) and U(t) are continuous adapted increasing processes. With the aid of the nonnegative semimartingale convergence theorem (Mao 2006), we have $\lim _{t\rightarrow \infty } Y(t) < \infty$, a.s. Thus, $\lim _{t\rightarrow \infty } V_{31}(t) < \infty$, a.s. This completes the proof. $\square$

1.2 Proof of Theorem 8

Let

$$\begin{aligned} {\tilde{A}}=\sum _{i=1}^{2}Q(x_{i})A_{i}, \end{aligned}$$

then model (26) becomes a three-dimensional model. We first prove that when ${\hat{\lambda }}(x_{i},h)>0$, $i=1,2$, the new three-dimensional model admits a unique stationary distribution. The proof is similar to the proof of Theorem 3 in “Appendix B”. Define a nonnegative ${\mathcal {C}}^{2}-$Lyapunov function

$$\begin{aligned} V_{4}=M_{3}V_{41}+V_{42}+V_{43}-V_{4}^{*}, \end{aligned}$$

$V_{4}^{*}$ is the minimum value of $V_{4}(P,{\tilde{A}},Z)$ and

$$\begin{aligned}&V_{41}=m_{8} Z-m_{9}\ln Z-m_{10}\ln {\tilde{A}}-(P+{\tilde{A}}+q(h)Z), \\&V_{32}=\frac{1}{\theta _{2}+2}(P+{\tilde{A}}+q(h)Z)^{\theta _{2}+2},\ \ \ \ V_{43}=-\ln P, \end{aligned}$$

where $m_{8}=(m_{10}\check{c}(x_{i},h)+s_{2}(h)q(h))/(s_{2}(h)+d_{2})$, $m_{9}=\frac{\check{s}_{1}\check{Q}}{\delta {\hat{c}}(x_{i},h)},$ $m_{10}=\frac{e}{ \alpha {\hat{\mu }} },$ $0<\theta _{2}<1$ is a constant which is small enough such that

$$\begin{aligned} h_{3}(x_{1}, x_{2}):=\min \{e,{\hat{s}}_{1},s_{2}\}-\frac{\theta _{2}+1}{2}\max \{\sigma _{1}^{2},{\check{\sigma }}_{2}^{2},\sigma _{3}^{2}\}>0, \end{aligned}$$

and the positive constant $M_{3}$ satisfying

$$\begin{aligned} -M_{3}{\hat{\lambda }}(x_{i},h)+f_{2}^{u}\le -2, \end{aligned}$$

where $f_{2}=-h_{3}(x_{1}, x_{2})(P+{\tilde{A}}+q(h)Z\big )^{\theta _{1}+2}+eP_\mathrm{in}\big (P+{\tilde{A}}+q(h)Z\big )^{\theta _{1}+1}+\alpha {\check{\mu }} {\tilde{A}}+e+\frac{\sigma _{1}^{2}}{2}$. By using the It$\hat{\mathrm {o}}$ formula, we can compute that

$$\begin{aligned} {\mathcal {L}}V_{4}\le&-M_{3}\left[ {\hat{\lambda }}(x_{i},h)-\frac{m_{8}\delta \check{c}(x_{i},h) }{ {\hat{Q}} }{\tilde{A}}Z\right] \\&\qquad -h_{3}(x_{1}, x_{2})(P+{\tilde{A}}+q(h)Z\big )^{\theta _{1}+2}+eP_\mathrm{in}\big (P+{\tilde{A}}+q(h)Z\big )^{\theta _{1}+1}\\&\qquad -\frac{eP_\mathrm{in}}{P}+\alpha {\check{\mu }} {\tilde{P}}+e+\frac{\sigma _{1}^{2}}{2}\\&\quad = -M_{3}\left[ {\hat{\lambda }}(x_{i},h)-\frac{m_{8}\delta \check{c}(x_{i},h) }{ {\hat{Q}} }{\tilde{A}}Z\right] -\frac{eP_\mathrm{in}}{P}+f_{2}(P,{\tilde{A}},Z). \end{aligned}$$

Then following the same logic as the proof in “Appendix B”, we can obtain that, when ${\hat{\lambda }}(x_{i},h)>0$, $i=1,2$, the new three-dimensional model has a unique stationary distribution $\nu (.)$ and it is ergodic. According to the ergodic property, the solution of the new three-dimensional model $\varPsi (t)=(P(t),{\tilde{A}}(t),Z(t))$ satisfying

$$\begin{aligned} {\mathcal {P}}\left\{ \langle \varPsi _{j}(x_{1},x_{2})\rangle =\int _{R_{+}^{3}}\psi _{j}\nu (\mathrm {d}\psi _{1},\mathrm {d}\psi _{2},\mathrm {d}\psi _{3}) \right\} =1,\ \ j=1,2,3. \end{aligned}$$

(D.74)

This suggests that the new three-dimensional model could achieve a relative stable state, i.e., the long-time mean persistent level of the concentration of dissolved inorganic phosphorus, the amount of organophosphorus in the entire algal species and the zooplankton carbon density remain unchanged a.s., no matter how many algal species exists in the model. If one of the algal species goes to extinction, then the dimorphic model (26) will ultimately degenerate to the original monomorphic model (20), and the following equations hold,

$$\begin{aligned} \langle P(x_{1},x_{2})\rangle =\langle P(x,h)\rangle ,\ \ \langle {\tilde{A}}(x_{1},x_{2})\rangle =Q(x)\langle A(x,h)\rangle , \ \ \langle Z(x_{1},x_{2})\rangle =\langle Z(x,h)\rangle ,\ \ \ a.s.\nonumber \\ \end{aligned}$$

(D.75)

We next prove that when $\min \{x_1,x_2\}<x<\max \{x_1,x_2\}$ holds, model (26) is stable in mean. To achieve that, we only need to prove $\langle A_{1}(x_1,x_2)\rangle >0$ and $\langle A_{2}(x_1,x_2)\rangle >0$. It is easy to compute that

$$\begin{aligned}&\mathrm {d}(P+{\tilde{A}}+q(h)Z)\\&\quad =\big [eP_\mathrm{in}-e P-s_{2}(h)q(h) Z-s_{1}(x_{1})Q(x_{1})A_{1}-s_{1}(x_{2})Q(x_{2})A_{2}\big ]\mathrm {d}t+\sigma _{1}N\mathrm {d}B_{1}(t)\\&\qquad +\sigma _{2}(x_{1},\rho _{1})Q(x_{1})P_{1}\mathrm {d}B_{2}(t) +\sigma _{2}(x_{2},\rho _{1})Q(x_{2})P_{3}\mathrm {d}B_{2}(t)+\sigma _{3}(h,\rho _{2})q(h)Z\mathrm {d}B_{4}(t). \end{aligned}$$

Integrating the above equation from 0 to t and dividing by t on both sides lead to

$$\begin{aligned} eP_\mathrm{in}-\frac{e}{t}\int _{0}^{t} P(s)\mathrm {d}s-\frac{s_{2}(h)q(h)}{t} \int _{0}^{t}Z(s)\mathrm {d}s-\sum _{i=1}^{2}\frac{s_{1}(x_{i})Q(x_{i})}{t}\int _{0}^{t}A_{i}(s)\mathrm {d}s=\frac{\psi _{2}(t)}{t},\nonumber \\ \end{aligned}$$

(D.76)

where $\psi _{2}(t)=R(t)-R(0)-H(t)$, $H(t)=\sigma _{1}\int _{0}^{t}N(s)\mathrm {d}B_{1}(s)+\sum _{i=1}^{2}\sigma _{2}(x_{i},\rho )Q(x_{i})\int _{0}^{t}P_{i}(s) \mathrm {d}B_{1+i}(s)+\sigma _{3}(h,\rho _{2})q(h)\int _{0}^{t}Z(s)\mathrm {d}B_{4}(s)$, $R(t)=P(t)+\sum _{i=1}^{2}Q(x_{i})A_{i}(t)+q(h)Z(t)$. According to Theorem 1 and the strong law of large numbers for local martingales (Mao 2006):

$$\begin{aligned} \lim _{t\rightarrow \infty }\frac{\psi _{2}(t)}{t}=0, \end{aligned}$$

then together with (D.75), (B.62) and the third equation of (B.61), we have

$$\begin{aligned}&\frac{s_{1}(x)Q(x)\kappa (h,\rho _{2})}{\delta c(x,h)}\nonumber \\&\quad =eP_\mathrm{in}-e\langle P(x,h)\rangle -s_{2}q(h)\langle Z(x,h)\rangle =eP_\mathrm{in}-e\langle P(x_{1},x_{2})\rangle -s_{2}(h)q(h)\langle Z(x_{1},x_{2})\rangle \nonumber \\&\quad =\lim _{t\rightarrow \infty }\sum _{i=1}^{2}\frac{s_{1}(x_{i})Q(x_{i})}{t}\int _{0}^{t}A_{i}(s)\mathrm {d}s. \end{aligned}$$

(D.77)

Moreover, from (D.75) and (B.62),

$$\begin{aligned} \begin{aligned} \frac{\kappa (h,\rho _{2}) Q(x)}{\delta c(x,h)}&=Q(x)\langle A(x,h)\rangle =\lim _{t\rightarrow \infty }\sum _{i=1}^{2}\frac{Q(x_{i})}{t}\int _{0}^{t}A_{i}(s)\mathrm {d}s. \end{aligned} \end{aligned}$$

(D.78)

Noticing the fact that, for any $f(t),g(t)<\infty$,

$$\begin{aligned} \liminf _{t\rightarrow \infty }[f(t)+g(t)]\le \liminf _{t\rightarrow \infty }f(t)+\limsup _{t\rightarrow \infty }g(t)\le \limsup _{t\rightarrow \infty }[f(t)+g(t)], \end{aligned}$$

and if $\lim _{t\rightarrow \infty }[f(t)+g(t)]$ exists, $\liminf _{t\rightarrow \infty }f(t)+\limsup _{t\rightarrow \infty }g(t)=\lim _{t\rightarrow \infty }[f(t)+g(t)]$. Then, together with (D.78) and (D.77),

$$\begin{aligned} {\left\{ \begin{array}{ll} Q(x_1)\overline{\langle A_{1}(x_1,x_2)\rangle }+Q(x_2)\underline{\langle A_{2}(x_1,x_2)\rangle }=\frac{\kappa (h,\rho _{2}) Q(x)}{\delta c(x,h)},\\ s_{1}(x_{1})Q(x_{1})\overline{\langle A_{1}(x_1,x_2)\rangle }+s_{1}(x_{2})Q(x_{2})\underline{\langle A_{2}(x_1,x_2)\rangle }=\frac{s_{1}(x)Q(x)\kappa (h,\rho _{2})}{\delta c(x,h)}, \end{array}\right. } \end{aligned}$$

(D.79)

and

$$\begin{aligned} {\left\{ \begin{array}{ll} Q(x_1)\underline{\langle A_{1}(x_1,x_2)\rangle }+Q(x_2)\overline{\langle A_{2}(x_1,x_2)\rangle }=\frac{\kappa (h,\rho _{2}) Q(x)}{\delta c(x,h)},\\ s_{1}(x_{1})Q(x_{1})\underline{\langle A_{1}(x_1,x_2)\rangle }+s_{1}(x_{2})Q(x_{2})\overline{\langle A_{2}(x_1,x_2)\rangle }=\frac{s_{1}(x)Q(x)\kappa (h,\rho _{2})}{\delta c(x,h)}. \end{array}\right. } \end{aligned}$$

(D.80)

When $\min \{x_1,x_2\}<x<\max \{x_1,x_2\}$, by solving (D.79) and (D.80) we can obtain that:

$$\begin{aligned}&\overline{\langle A_{1}(x_1,x_2)\rangle }=\underline{\langle A_{1}(x_1,x_2)\rangle }=\frac{\kappa (h,\rho _{2}) Q(x)}{\delta c(x,h)Q(x_1)}\frac{s(x_2)-s(x)}{s(x_2)-s(x_1)},\\&\overline{\langle A_{2}(x_1,x_2)\rangle }=\underline{\langle A_{2}(x_1,x_2)\rangle }=\frac{\kappa (h,\rho _{2}) Q(x)}{\delta c(x,h)Q(x_2)}\frac{s(x_1)-s(x)}{s(x_1)-s(x_2)}. \end{aligned}$$

According to the description of x-functions in Table 1, $s'(x)>0$, then $\min \{s(x_1),s(x_2)\}<s(x)<\max \{s(x_1),s(x_2)\}$. Consequently, we can conclude that $\langle A_{1}(x_1,x_2)\rangle >0$ and $\langle A_{2}(x_1,x_2)\rangle >0$.

Our next step is to explicitly express the persistence level of model (26) in terms of $x_1$ and $x_2$ only. Following similar methods as in the derivation of (B.62), the solution of model (26) satisfying $\lim _{t\rightarrow \infty }\frac{\ln A_{i}(t)}{t}=0$, $\lim _{t\rightarrow \infty }\frac{\ln Z(t)}{t}=0$, $i=1,2$. Together with (D.77) leads to:

$$\begin{aligned} {\left\{ \begin{array}{ll} \begin{aligned} &{}\alpha \mu (x_1)\langle P(x_1,x_2)\rangle -c(x_1,h)\langle Z(x_1,x_2)\rangle =\eta (x_1,\rho _{1}),\\ &{}\alpha \mu (x_2)\langle P(x_1,x_2)\rangle -c(x_2,h)\langle Z(x_1,x_2)\rangle =\eta (x_2,\rho _{1}),\\ &{}\delta c(x_1,h) \langle A_1(x_1,x_2)\rangle +\delta c(x_2,h) \langle A_2(x_1,x_2)\rangle =\kappa (h,\rho _{2}),\\ &{}e\langle P(x_1,x_2)\rangle +s_{1}(x_1)Q(x_1)\langle A_1(x_1,x_2)\rangle +s_{1}(x_2)Q(x_2)\langle A_2(x_1,x_2)\rangle +\\ {} &{} s_{2}(y)q(y)\langle Z(x_1,x_2)\rangle =eP_\mathrm{in}. \end{aligned} \end{array}\right. } \end{aligned}$$

(D.81)

Through simple calculation, (D.81) has a unique positive solution as follows:

$$\begin{aligned} \begin{aligned}&\langle P(x_{1},x_{2})\rangle =\frac{c(x_{2},h)\eta (x_{1},\rho _{1})-c(x_{1},h)\eta (x_{2},\rho _{1})}{\alpha c(x_{2},h)\mu (x_{1})-\alpha c(x_{1},h)\mu (x_{2})},\\&\langle A_1(x_{1},x_{2})\rangle =\frac{s_{1}(x_{2})Q(x_{2})\kappa (h,\rho _{2})-\delta c(x_{2},h)f_{3}}{\delta c(x_{1},h)s_{1}(x_{2})Q(x_{2})-\delta c(x_{2},h)s_{1}(x_{1})Q(x_{1})},\\&\langle A_2(x_{1},x_{2})\rangle =\frac{\delta c(x_{1},h)f_{3}-s_{1}(x_{1})Q(x_{1})\kappa (h,\rho _{2})}{\delta c(x_{1},h)s_{1}(x_{2})Q(x_{2})-\delta c(x_{2},h)s_{1}(x_{1})Q(x_{1})},\\&\langle Z(x_{1},x_{2})\rangle =\frac{\mu (x_{2}) \eta (x_{1},\rho _{1})-\mu (x_{1})\eta (x_{2},\rho _{1})}{c(x_{2})\mu (x_{1})-c(x_{1})\mu (x_{2})}, \end{aligned} \end{aligned}$$

(D.82)

where $f_{3}=eP_\mathrm{in}-e\langle P(x_{1},x_{2})\rangle -s_{2}(h)q(h)\langle Z(x_{1},x_{2})\rangle$. The proof is completed.

Remark 10

Note that when zooplankton species also evolves, Theorem 7 and 8 also hold, and the long-time mean persistent level becomes:

$$\begin{aligned} \begin{aligned}&\langle P(x_{1},x_{2},y)\rangle =\frac{c(x_{2},y)\eta (x_{1},\rho _{1})-c(x_{1},y)\eta (x_{2},\rho _{1})}{\alpha c(x_{2},y)\mu (x_{1})-\alpha c(x_{1},y)\mu (x_{2})},\\&\langle A_1(x_{1},x_{2},y)\rangle =\frac{s_{1}(x_{2})Q(x_{2})\kappa (y,\rho _{2})-\delta c(x_{2},y)f_{4}}{\delta c(x_{1},y)s_{1}(x_{2})Q(x_{2})-\delta c(x_{2},y)s_{1}(x_{1})Q(x_{1})},\\&\langle A_2(x_{1},x_{2},y)\rangle =\frac{\delta c(x_{1},y)f_{4}-s_{1}(x_{1})Q(x_{1})\kappa (y,\rho _{2})}{\delta c(x_{1},y)s_{1}(x_{2})Q(x_{2})-\delta c(x_{2},y)s_{1}(x_{1})Q(x_{1})},\\&\langle Z(x_{1},x_{2},y)\rangle =\frac{\mu (x_{2}) \eta (x_{1},\rho _{1})-\mu (x_{1})\eta (x_{2},\rho _{1})}{c(x_{2},y)\mu (x_{1})-c(x_{1},y)\mu (x_{2})}, \end{aligned} \end{aligned}$$

(D.83)

where $f_{4}=eP_\mathrm{in}-e\langle P(x_{1},x_{2},y)\rangle -s_{2}(y)q(y)\langle Z(x_{1},x_{2},y)\rangle$.

Appendix E: Proof of Theorem 9

In this part, we mainly investigate the convergence stability and the evolutionary stability of $(x_{1}^{*}, x_{2}^{*})$. Taking the partial derivatives of the first and second equations of (D.81), with respect to $x_{1}$, $x_{2}$, respectively, leads to

$$\begin{aligned}&\alpha \mu '(x_1)\langle P(x_1,x_2)\rangle +\alpha \mu (x_1)\frac{\partial \langle P(x_1,x_2)\rangle }{\partial x_1}\\&\quad -c'(x_1,h)\langle Z(x_1,x_2)\rangle -c(x_1,h)\frac{\partial \langle Z(x_1,x_2)\rangle }{\partial x_1}= \eta '(x_1,\rho _{1}), \\&\alpha \mu '(x_2)\langle P(x_1,x_2)\rangle +\alpha \mu (x_2)\frac{\partial \langle P(x_1,x_2)\rangle }{\partial x_2}\\&\quad -c'(x_2,h)\langle Z(x_1,x_2)\rangle -c(x_2,h)\frac{\partial \langle Z(x_1,x_2)\rangle }{\partial x_2}=\eta '(x_2,\rho _{1}). \end{aligned}$$

Together with (30) leads to

$$\begin{aligned} \begin{aligned} c(x_{i}^{*},h)\frac{\partial \langle Z(x_{1},x_{2})\rangle }{\partial x_{i}}\bigg |_{x_{i}=x_{i}^{*}}=\alpha \mu (x_{i}^{*})\frac{\partial \langle P(x_{1},x_{2})\rangle }{\partial x_{i}}\bigg |_{x_{i}=x_{i}^{*}},\ \ i=1,2. \end{aligned} \end{aligned}$$

(E.84)

Also from the first and second equations of (D.81), we can compute that

$$\begin{aligned} \begin{aligned} c(x_{i},h)\frac{\partial \langle Z(x_{1},x_{2})\rangle }{\partial x_{j}}&=\alpha \mu (x_{i})\frac{\partial \langle P(x_{1},x_{2})\rangle }{\partial x_{j}}, \ i,j=1,2,\ i\ne j. \end{aligned} \end{aligned}$$

(E.85)

Combining Eq. (E.85) with (E.84) leads to,

$$\begin{aligned} c(x_{i}^{*},h)\left( \frac{\mu (x_{j}^{*})}{\mu (x_{i}^{*})}-\frac{c(x_{j}^{*},h)}{c(x_{i}^{*},h)}\right) \frac{\partial \langle Z(x_{1},x_{2})\rangle }{\partial x_{i}}\bigg |_{x_{i}=x_{i}^{*}}=0, \end{aligned}$$

(E.86)

where $i,j=1,2$ and $i\ne j$. Moreover, from the expression of $\langle Z(x_{1},x_{2})\rangle$ in (D.82), we know that

$$\begin{aligned}{}[c(x_{i},h)\mu (x_{j})-c(x_{j},h)\mu (x_{i})][\eta (x_{j},\rho _{1})\mu (x_{i})-\eta (x_{i},\rho _{1})\mu (x_{j})]>0, \end{aligned}$$

otherwise, it will contradict with (D.74). Then, from (E.86), we have

$$\begin{aligned} \begin{aligned} \frac{\partial \langle Z(x_{1},x_{2})\rangle }{\partial x_{i}}\bigg |_{x_{i}=x_{i}^{*}}=\frac{\partial \langle P(x_{1},x_{2})\rangle }{\partial x_{i}}\bigg |_{x_{i}=x_{i}^{*}}=0,\ \ i=1,2. \end{aligned} \end{aligned}$$

(E.87)

Let

$$\begin{aligned} h_{1i}(x_{1}, x_{2})=&c'(x_{i},h)\frac{\partial \langle Z(x_{1},x_{2})\rangle }{\partial x_{i}}-\alpha \mu '(x_{i})\frac{\partial \langle P(x_{1},x_{2})\rangle }{\partial x_{i}},\\ h_{2i}(x_{1}, x_{2})=&\alpha \mu (x_{i})\frac{\partial ^{2}\langle P(x_{1},x_{2})\rangle }{\partial x_{i}^{2}}-c(x_{i},h)\frac{\partial ^{2}\langle Z(x_{1},x_{2})\rangle }{\partial x_{i}^{2}}. \end{aligned}$$

Obviously, $h_{1i}(x_{1}, x_{2})=0$, $i=1,2$. The evolutionary singular dimorphism $(x_{1}^{*}, x_{2}^{*})$ is evolutionary stable if for both $i\in \{1,2\}$, the following inequality holds,

$$\begin{aligned} \begin{aligned} \frac{\partial F_{3}^{2}(x_{1},x_{2},x_\mathrm{mut})}{\partial x_\mathrm{mut}^{2}}\bigg | _{y=x=x_{i}^{*}}=&\alpha \mu ''(x_{i}^{*})\langle P(x_{1}^{*},x_{2}^{*})\rangle -c''(x_{i}^{*},h)\langle Z(x_{1}^{*},x_{2}^{*})\rangle -\frac{\partial ^{2}\eta (x_{i},\rho _{1})}{\partial x_{i}^{2}}\bigg |_{x_{i}=x_{i}^{*}}\\ =&2h_{1i}(x_{1}^{*}, x_{2}^{*})- h_{2i}(x_{1}^{*}, x_{2}^{*})\\ =&- h_{2i}(x_{1}^{*}, x_{2}^{*})\\ <&0. \end{aligned} \end{aligned}$$

(E.88)

Furthermore, the local convergence stability at the evolutionary singular dimorphism $(x_{1}^{*}, x_{2}^{*})$ is determined by the Jacobian matrix of evolutionary dynamics (29) at this point. The Jacobian matrix ${\mathcal {J}}$ of (29) at $(x_{1}^{*}, x_{2}^{*})$ is given by:

$$\begin{aligned} \begin{aligned} {\mathcal {J}}&=\left[ \begin{aligned} \displaystyle { m_{1}(x_{1}^{*},x_{2}^{*})\frac{\partial D_{31}(x_{1},x_{2})}{\partial x_{1}}~~~~ m_{1}(x_{1}^{*},x_{2}^{*})\frac{\partial D_{31}(x_{1},x_{2})}{\partial x_{2}}}\\ \displaystyle { m_{2}(x_{1}^{*},x_{2}^{*})\frac{\partial D_{32}(x_{1}^{*},x_{2}^{*})}{\partial x_{1}}~~~~ m_{2}(x_{1}^{*},x_{2}^{*})\frac{\partial D_{32}(x_{1},x_{2})}{\partial x_{2}}} \end{aligned} \right] _{x_{1}=x_{1}^{*},\ x_{2}=x_{2}^{*}}\\&=\left[ \begin{aligned} \displaystyle {m_{1}(x_{1}^{*},x_{2}^{*})\frac{\partial D_{31}(x_{1},x_{2})}{\partial x_{1}}\bigg |_{\begin{array}{c} x_{1}=x_{1}^{*}\\ x_{2}=x_{2}^{*} \end{array}}~~~~~~~~~~~~~~~~~~~~~0~~~~~~~~~~~~~~~~~~~~}\\ \displaystyle {~~~~~~~~~~~~~~~0~~~~~~~~~~~~~~~~~~~~~~~~~ m_{2}(x_{1}^{*},x_{2}^{*})\frac{\partial D_{32}(x_{1},x_{2})}{\partial x_{2}}\bigg |_{\begin{array}{c} x_{1}=x_{1}^{*}\\ x_{2}=x_{2}^{*} \end{array}}} \end{aligned} \right] \\&=\left[ \begin{aligned} \displaystyle {-m_{1}(x_{1}^{*},x_{2}^{*})h_{21}(x_{1}^{*}, x_{2}^{*})~~~~~~~~~~~~~~~~~~~~0~~~~~~~~~~~~}\\ \displaystyle {~~~~~~~~~~~~~~0~~~~~~~~~ -m_{2}(x_{1}^{*},x_{2}^{*})h_{22}(x_{1}^{*}, x_{2}^{*})} \end{aligned} \right] , \end{aligned} \end{aligned}$$

(E.89)

then if the Jacobian matrix satisfies $\mathrm {det}({\mathcal {J}}) > 0$ and $\mathrm {tr}({\mathcal {J}}) < 0$, i.e., $h_{2i}(x_{1}^{*}, x_{2}^{*})>0$ hold for both $i\in \{1,2\}$, the evolutionary singular dimorphism is locally convergence stable. From the above analysis, it is obvious that once the evolutionary singular dimorphism $(x_{1}^{*}, x_{2}^{*})$ is convergence stable, it must also be evolutionary stable, i.e., it is a CSS. Otherwise, i.e., there at least exists one $i\in \{1,2\}$, such that $h_{2i}(x_{1}^{*}, x_{2}^{*})<0$, the evolutionary singular dimorphism is a repellor. This completes the proof.

Appendix F: Proof of Remark 6

Proof

When $s_{2}(y)=\frac{1}{\beta _{3}}\exp (k_{2}y)$, in order to prove $y^{*}$ is always an ESS, we only need to prove that inequality (38) is always holds. Combing the second equation of (35) and the expression of $\langle A(x,y)\rangle$ in (9) leads to

$$\begin{aligned} \langle A(x^{*},y^{*})\rangle '_{y}=0, \end{aligned}$$

(F.90)

and

$$\begin{aligned} \frac{ \kappa '(y^*,\rho _2)}{\delta c'_{y}(x^{*},y^{*})}=\frac{ \kappa (y^*,\rho _2)}{\delta c(x^*,y^*)}=\langle A(x^{*},y^{*})\rangle , \end{aligned}$$

(F.91)

where

$$\begin{aligned}&c(x,y)=c_\mathrm{m}\exp (-v(x-\theta y)^{2}),\ \ \ \kappa (y)=\frac{1}{\beta _{3}}\exp (k_{2}y)+d_{2}+\frac{{\tilde{m}}}{2}\exp (-2k_{7}y),\\&{\tilde{m}}=\left( \frac{k_{5}\rho _{2}}{k_{6}+\rho _{2}}\right) ^{2}. \end{aligned}$$

By calculation,

$$\begin{aligned}&c'_{y}(x^{*},y^{*}) =2c_\mathrm{m}\theta v(x^{*}-\theta y^{*})\exp (-v(x^{*}-\theta y^{*})^{2}), \end{aligned}$$

(F.92)

$$\begin{aligned}&c''_{y}(x^{*},y^{*}) =2v\theta ^{2}c_\mathrm{m}\big (2v(x^{*}-\theta y^{*})^{2}-1\big )\exp (-v(x^{*}-\theta y^{*})^{2}), \end{aligned}$$

(F.93)

$$\begin{aligned}&\kappa '(y^{*}) =\frac{k_{2}}{\beta _3}\exp (k_{2}y^{*})-{\tilde{m}}k_{7}\exp (-2k_{7}y^{*}), \end{aligned}$$

(F.94)

$$\begin{aligned}&\kappa ''(y^{*}) =\frac{k_{2}^{2}}{\beta _3}\exp (k_{2}y^{*})+2{\tilde{m}}k_{7}^{2}\exp (-2k_{7}y^{*})>0. \end{aligned}$$

(F.95)

If $0<2v(x^{*}-\theta y^{*})^{2}\le 1$, $c''_{y}(x^{*},y^{*})\le 0$, then according to (38) and (F.95), we have,

$$\begin{aligned} \frac{\partial ^2 F_{5}(x,y,y_\mathrm{mut})}{\partial y_\mathrm{mut}^2}\bigg |_{\begin{array}{c} x=x^{*} \end{array}}<0. \end{aligned}$$

We then prove if $2v(x^{*}-\theta y^{*})^{2}> 1$, $\frac{\partial ^2 F_{5}(x,y,y_\mathrm{mut})}{\partial y_\mathrm{mut}^2}\bigg |_{\begin{array}{c} x=x^{*} \end{array}}<0$ also holds.

When $2v(x^{*}-\theta y^{*})^{2}> 1$, substituting (F.91)-(F.95) into (38) leads to

$$\begin{aligned} \begin{aligned}&\frac{\partial ^2 F_{5}(x,y,y_\mathrm{mut})}{\partial y_\mathrm{mut}^2}\bigg |_{\begin{array}{c} x=x^{*}\\ y_\mathrm{mut}=y=y^{*} \end{array}}\\&\quad =\delta c''_{y}(x^{*},y^{*}) \langle A(x^{*},y^{*})\rangle -\kappa ''(y^*,\rho _2)\\&\quad = c''_{y}(x^{*},y^{*})\frac{\kappa '(y^*,\rho _2) }{c'_{y}(x^{*},y^{*})} -\kappa ''(y^*,\rho _2)\\&\quad =\theta \big (2v(x^{*}-\theta y^{*})^{2}-1\big ) \frac{\frac{k_{2}}{\beta _3}\exp (k_{2}y^{*})-{\tilde{m}}k_{7}\exp (-2k_{7}y^{*}) }{x^{*}-\theta y^{*}} -\kappa ''(y^*,\rho _2)\\&\quad <2\theta v(x^{*}-\theta y^{*})^{2} \frac{\frac{k_{2}}{\beta _3}\exp (k_{2}y^{*})-{\tilde{m}}k_{7}\exp (-2k_{7}y^{*}) }{x^{*}-\theta y^{*}} -\kappa ''(y^*,\rho _2). \end{aligned} \end{aligned}$$

(F.96)

We then divide our proof into two cases.

(1). $x^*>\theta y^*$, then $\frac{k_{2}}{\beta _3}\exp (k_{2}y^{*})-{\tilde{m}}k_{7}\exp (-2k_{7}y^{*})>0$ and from (F.91),

$$\begin{aligned} k_{2}>\frac{\frac{k_{2}}{\beta _3}\exp (k_{2}y^{*})-{\tilde{m}}k_{7}\exp (-2k_{7}y^{*})}{\frac{1}{\beta _3}\exp (k_{2}y^{*})+d_{2}+\frac{{\tilde{m}}}{2}\exp (-2k_{7}y^{*})} =2\theta v(x^{*}-\theta y^{*})>0. \end{aligned}$$

(F.97)

Substitute (F.97) into (F.96), we have

$$\begin{aligned} \begin{aligned} \frac{\partial ^2 F_{5}(x,y,y_\mathrm{mut})}{\partial y_\mathrm{mut}^2}\bigg |_{\begin{array}{c} x=x^{*}\\ y_\mathrm{mut}=y=y^{*} \end{array}}&<2\theta v(x^{*}-\theta y^{*})\frac{k_{2}}{\beta _3}\exp (k_{2}y^{*})-\frac{k_{2}^{2}}{\beta _3}\exp (k_{2}y^{*})\\&<\frac{k_{2}^{2}}{\beta _3}\exp (k_{2}y^{*})-\frac{k_{2}^{2}}{\beta _3}\exp (k_{2}y^{*})\\&=0. \end{aligned} \end{aligned}$$

(F.98)

(2). $x^*<\theta y^*$, then $\frac{k_{2}}{\beta _3}\exp (k_{2}y^{*})-{\tilde{m}}k_{7}\exp (-k_{7}y^{*})<0$ and from (F.91),

$$\begin{aligned} 2k_{7}>\frac{{\tilde{m}}k_{7}\exp (-2k_{7}y^{*})-\frac{k_{2}}{\beta _3}\exp (k_{2}y^{*})}{\frac{1}{\beta _3}\exp (k_{2}y^{*})+d_{2}+\frac{{\tilde{m}}}{2}\exp (-2k_{7}y^{*})} =2\theta v(\theta y^{*}-x^{*})>0. \end{aligned}$$

(F.99)

then

$$\begin{aligned} \begin{aligned} \frac{\partial ^2 F_{5}(x,y,y_\mathrm{mut})}{\partial y_\mathrm{mut}^2}\bigg |_{\begin{array}{c} x=x^{*}\\ y_\mathrm{mut}=y=y^{*} \end{array}}&<2 v\theta (\theta y^{*}-x^{*}) {\tilde{m}}k_{7}\exp (-2k_{7}y^{*})-2{\tilde{m}}k_{7}^{2}\exp (-2k_{7}y^{*})\\&<2 {\tilde{m}}k_{7}^{2}\exp (-2k_{7}y^{*})-2{\tilde{m}}k_{7}^{2}\exp (-2k_{7}y^{*})\\&=0. \end{aligned} \end{aligned}$$

(F.100)

Then, $\frac{\partial ^2 F_{5}(x,y,y_\mathrm{mut})}{\partial y_\mathrm{mut}^2}\bigg |_{\begin{array}{c} x=x^{*}\\ y_\mathrm{mut}=y=y^{*} \end{array}}<0$ always holds. The proof is thus completed. $\square$