A discrete-time nutrients-phytoplankton-oysters mathematical model of a bay ecosystem*

Abstract

Populations are generally censused daily, weekly, monthly or annually. In this paper, we introduce a discrete-time nutrients-phytoplankton-oysters (NPO) model that describes the interactions of nutrients, phytoplankton and oysters in a bay ecosystem. We compute the threshold parameter ${\mathcal{R}}_N$ for persistence of phytoplankton with or without oysters. When ${\mathcal{R}}_N<1$, then both phytoplankton and oysters populations go extinct. However, when ${\mathcal{R}}_N>1$, we show that the model may exhibit two scenarios: (1) a locally asymptotically stable equilibrium with positive values of nutrients and phytoplankton with oysters missing, and (2) a locally asymptotically stable interior equilibrium with positive values of nutrients, phytoplankton and oysters. We use sensitivity analysis to study the impact of human and environmental factors on the model. We use examples to illustrate that some human activities and environmental factors can force the interior equilibrium to undergo a Neimark–Sacker bifurcation which generates phytoplankton blooms with oscillations in oysters population and nutrients level.

Keywords:

• Discrete-time
• nutrients
• phytoplankton
• oysters

Mathematics Subject Classification:

• 92-10

1. Introduction

Estuaries, for example the Chesapeake Bay (the largest estuary in the US), provide diverse habitats for wildlife and aquatic life while shielding adjacent communities against flooding, decrease waterways and pollution. Furthermore, estuaries provide opportunities for commercial and recreational activities that support local economies [2]. The observed decrease in the Chesapeake Bay's health is associated with the increase of the human population in its watershed [3]. The large population causes waste, utilization of the resources and the consequent changes of the character of the land, water and air of the bay [1]. For example, loss of habitat, overfishing and disease have depleted oyster population in the Chesapeake Bay. Excess nutrients cause algae blooms that block sunlight which endangers aquatic species. The excess nutrients and sedimentation threaten bay water quality [1].

Interactions of nutrients and phytoplankton have been investigated in several research papers using continuous-time systems of ordinary differential equations models. For example, Fan and Glibert in [12] and Hood et al. in [14] used continuous-time models to study the impact of nutrients on the formation of phytoplankton blooms. A continuous-time model of nutrients-phytoplankton (NP) with phytoplankton bloom was analysed by Huppert et al. in [15, 16]. In [8–10, 13, 19, 21], nutrients-phytoplankton-zooplankton (NPZ) continuous-time models were studied. The authors demonstrated that a wide range of dynamical behaviours from equilibrium to oscillatory dynamics can be generated in NPZ models. Moreover, the stability analysis of a more general type of continuous-time nutrients-phytoplankton-oysters models was studied by Saijnders and Bazin [20]. Ziyadi et al. in [25], used sensitivity analysis in a continuous-time nutrients-phytoplankton-oysters (NPO) model to investigate the effect of human and environmental factors on the NPO model dynamics.

Populations are typically censused at discrete-time intervals, for example, daily, weekly, monthly or annually. To capture the discrete nature of population surveillance data, in this paper, unlike in [25], we introduce a discrete-time nutrients-phytoplankton-oysters model that describes nutrients, phytoplankton and oysters interaction dynamics in a bay ecosystem. The model is a system of three difference equations. In our NPO model, we use the threshold parameter, ${\mathcal{R}}_N$, to study persistence of phytoplankton with or without oysters in the system. In particular, we obtain that when ${\mathcal{R}}_N<1$, then both phytoplankton and oysters go extinct. However, when ${\mathcal{R}}_N>1$, we show that the model is capable of exhibiting two scenarios: (1) oyster population goes extinct or the occurrence of a locally asymptotically stable equilibrium point with positive values of nutrients level and phytoplankton biomas with oysters missing, and (2) the nutrients, phytoplankton and oyster persist or the occurrence of a locally asymptotically stable interior equilibrium with positive values of nutrients, phytoplankton and oysters. As in [25], we use sensitivity analysis to study the impact of human and environmental factors on the discrete-time NPO model.

This paper is organized as follows: In Section 2, we introduce the discrete-time mathematical model that describes interaction dynamics of nutrients, phytoplankton and oysters in a bay ecosystem. We establish the boundedness of nutrients level and all populations in Section 3. In Section 4, we give the condition for persistence of nutrients in the model. In addition, we give an illustrative example that shows that increasing model parameters, such as nutrients flow, can force persistence of the phytoplankton population. Furthermore, we perform sensitivity analysis. In Section 5, we give conditions for persistence of both nutrients and phytoplankton in the model. Also, we give an illustrative example that shows that increasing model parameters, such as nutrients flow, can force persistence of oysters population. In Section 6, we present illustrative examples showing that changing model parameters such as increasing the flow rate of nutrients or increasing the survival of phytoplankton can lead to the occurrence of Neimark–Sacker bifurcation. Finally, we summarize our results in the Discussion section.

2. Discrete-time nutrients-phytoplankton-oysters model

In this section, we construct the discrete-time NPO model. However, we first list the model variables and parameters in Tables 1 and 2, respectively.

Table 1. List of model variables.

Variable	Description	Units
t	Time ($t=0,1,2,\dots$)	$\mathrm{d}$ (Day)
$N_t$	Nutrients density at time t	$\mathrm{mg}\mathrm{N}{\mathrm{l}}^{-1}$
$P_t$	Phytoplankton density at time t	$\mathrm{mg}\mathrm{N}{\mathrm{l}}^{-1}$
$O_t$	Oysters density at time t	$\mathrm{mg}\mathrm{N}{\mathrm{l}}^{-1}$

Table 2. Table of model parameters.

Parameter	Description	Unit	Range	Baseline	Reference
γ	Flow rate of N	$\mathrm{mg}\mathrm{N}{\mathrm{l}}^{-1}{\mathrm{d}}^{-1}$	0.01 – 5	2.505	[25]
${\mu }_N$	Escapement rate of N	${\mathrm{d}}^{-1}$	0.1 – 1	0.7	Estimated
${\mu }_P$	Survival rate of P	${\mathrm{d}}^{-1}$	0.1 – 1	0.6	Estimated
${\mu }_O$	Survival rate of O	${\mathrm{d}}^{-1}$	0.01 – 0.1	0.03	Estimated
α	N uptake by P	${\mathrm{d}}^{-1}$	0.1 – 1.5	0.8	[12]
δ	P filtration by O	${\mathrm{d}}^{-1}$	0.1 – 2	1.05	[25]

In [25], Figure 1 was used to construct a continuous-time ordinary differential equations model for the interactions of nutrients, phytoplankton and oysters in a bay ecosystem.

Next, we use Figure 1 together with the variables and parameters of Tables 1 and 2 to introduce a discrete-time NPO model in three steps. First, we introduce a discrete-time equation that describes the nutrients level of our discrete-time NPO model.

2.1. (i) Nutrients level equation

The nutrients level at time $(t+1)$ in the bay ecosystem is modelled by Equation (1(1) $N_{t+1}=\underset{\begin{array}{c}N\mathrm{flow}\mathrm{not}\\ \mathrm{consumed}\mathrm{by}P\end{array}}{\underbrace{\gamma {\mathrm{e}}^{-\alpha P_t}}}+\underset{\begin{array}{c}N\mathrm{not}\mathrm{lost}\mathrm{to}\mathrm{sinking}\\ \mathrm{and}\mathrm{denitrification}\end{array}}{\underbrace{{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t}}},$(1) ), (1) $N_{t+1}=\underset{\begin{array}{c}N\mathrm{flow}\mathrm{not}\\ \mathrm{consumed}\mathrm{by}P\end{array}}{\underbrace{\gamma {\mathrm{e}}^{-\alpha P_t}}}+\underset{\begin{array}{c}N\mathrm{not}\mathrm{lost}\mathrm{to}\mathrm{sinking}\\ \mathrm{and}\mathrm{denitrification}\end{array}}{\underbrace{{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t}}},$(1) for $t=0,1,2,3,\dots$, where $0<{\mu }_N<1$, $\gamma >0$ and $\alpha >0$. In the unit time interval, first term of Equation (1(1) $N_{t+1}=\underset{\begin{array}{c}N\mathrm{flow}\mathrm{not}\\ \mathrm{consumed}\mathrm{by}P\end{array}}{\underbrace{\gamma {\mathrm{e}}^{-\alpha P_t}}}+\underset{\begin{array}{c}N\mathrm{not}\mathrm{lost}\mathrm{to}\mathrm{sinking}\\ \mathrm{and}\mathrm{denitrification}\end{array}}{\underbrace{{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t}}},$(1) ) represents the proportion of nutrients flow that was not consumed by phytoplankton, and second term models nutrients not consumed by phytoplankton and not lost to sinking and denitrification. As in the continuous-time NPO model of [25], when there are no phytoplankton and $P_t=0$, then the nutrients level is asymptotically constant in model (1(1) $N_{t+1}=\underset{\begin{array}{c}N\mathrm{flow}\mathrm{not}\\ \mathrm{consumed}\mathrm{by}P\end{array}}{\underbrace{\gamma {\mathrm{e}}^{-\alpha P_t}}}+\underset{\begin{array}{c}N\mathrm{not}\mathrm{lost}\mathrm{to}\mathrm{sinking}\\ \mathrm{and}\mathrm{denitrification}\end{array}}{\underbrace{{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t}}},$(1) ).

Figure 1. Interactions of nutrients, phytoplankton and oysters in a bay ecosystem [25].

In step 2, we introduce the discrete-time equation that describes the phytoplankton biomass at time $(t+1)$.

2.2. (ii) Phytoplankton biomass equation

The phytoplankton biomass at time $(t+1)$ in the bay ecosystem is modelled by Equation (2(2) $P_{t+1}=\underset{\begin{array}{c}P\mathrm{gain}\mathrm{by}N\mathrm{uptake}\\ \mathrm{and}\mathrm{escaped}O\mathrm{filtration}\end{array}}{\underbrace{N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}}}+\underset{\begin{array}{c}P\mathrm{not}\mathrm{lost}\mathrm{to}\mathrm{sinking}\\ \mathrm{and}\mathrm{predation}\end{array}}{\underbrace{{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t}}},$(2) ), (2) $P_{t+1}=\underset{\begin{array}{c}P\mathrm{gain}\mathrm{by}N\mathrm{uptake}\\ \mathrm{and}\mathrm{escaped}O\mathrm{filtration}\end{array}}{\underbrace{N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}}}+\underset{\begin{array}{c}P\mathrm{not}\mathrm{lost}\mathrm{to}\mathrm{sinking}\\ \mathrm{and}\mathrm{predation}\end{array}}{\underbrace{{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t}}},$(2) for $t=0,1,2,3,\dots$, where $0<{\mu }_P<1$, $\beta >0$ and $\delta >0$. In the unit time interval, first term of Equation (2(2) $P_{t+1}=\underset{\begin{array}{c}P\mathrm{gain}\mathrm{by}N\mathrm{uptake}\\ \mathrm{and}\mathrm{escaped}O\mathrm{filtration}\end{array}}{\underbrace{N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}}}+\underset{\begin{array}{c}P\mathrm{not}\mathrm{lost}\mathrm{to}\mathrm{sinking}\\ \mathrm{and}\mathrm{predation}\end{array}}{\underbrace{{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t}}},$(2) ) represents gain of phytoplankton from nutrients uptake while escaping from oysters filtration, and second term represents phytoplankton not consumed by oysters and not lost to sinking and predation. As in the continuous-time NPO model of [25], the phytoplankton population goes extinct when $N_t=0$ and the nutrients level is very low.

In step 3, we introduce the discrete-time equation that describes the oysters biomass at time $(t+1)$.

2.3. (iii) Oysters biomass equation

The oysters biomass at time $(t+1)$ in the bay ecosystem is modelled by Equation (3(3) $O_{t+1}=\underset{O\mathrm{gain}\mathrm{from}\mathrm{filtration}\mathrm{of}P}{\underbrace{P_t(1-{\mathrm{e}}^{-\delta O_t})}}+\underset{\begin{array}{c}O\mathrm{not}\mathrm{lost}\mathrm{tomortality},\\ \mathrm{harvest},\mathrm{disease}\mathrm{and}\mathrm{respiration}\end{array}}{\underbrace{{\mu }_O{O}_t}},$(3) ), (3) $O_{t+1}=\underset{O\mathrm{gain}\mathrm{from}\mathrm{filtration}\mathrm{of}P}{\underbrace{P_t(1-{\mathrm{e}}^{-\delta O_t})}}+\underset{\begin{array}{c}O\mathrm{not}\mathrm{lost}\mathrm{tomortality},\\ \mathrm{harvest},\mathrm{disease}\mathrm{and}\mathrm{respiration}\end{array}}{\underbrace{{\mu }_O{O}_t}},$(3) for $t=0,1,2,3,\dots$ where $0<{\mu }_0<1$, and $\eta >0$. In the unit time interval, first term of Equation (3(3) $O_{t+1}=\underset{O\mathrm{gain}\mathrm{from}\mathrm{filtration}\mathrm{of}P}{\underbrace{P_t(1-{\mathrm{e}}^{-\delta O_t})}}+\underset{\begin{array}{c}O\mathrm{not}\mathrm{lost}\mathrm{tomortality},\\ \mathrm{harvest},\mathrm{disease}\mathrm{and}\mathrm{respiration}\end{array}}{\underbrace{{\mu }_O{O}_t}},$(3) ) represents the proportion of oysters gained by filtration of phytoplankton, and second term models the oysters not lost to mortality, harvest, disease and respiration at time $(t+1)$. Also, as in the continuous-time NPO model of [25], the oyster population goes extinct when $P_t=0$ and there are no phytoplankton.

From Equations (1(1) $N_{t+1}=\underset{\begin{array}{c}N\mathrm{flow}\mathrm{not}\\ \mathrm{consumed}\mathrm{by}P\end{array}}{\underbrace{\gamma {\mathrm{e}}^{-\alpha P_t}}}+\underset{\begin{array}{c}N\mathrm{not}\mathrm{lost}\mathrm{to}\mathrm{sinking}\\ \mathrm{and}\mathrm{denitrification}\end{array}}{\underbrace{{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t}}},$(1) ), (2(2) $P_{t+1}=\underset{\begin{array}{c}P\mathrm{gain}\mathrm{by}N\mathrm{uptake}\\ \mathrm{and}\mathrm{escaped}O\mathrm{filtration}\end{array}}{\underbrace{N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}}}+\underset{\begin{array}{c}P\mathrm{not}\mathrm{lost}\mathrm{to}\mathrm{sinking}\\ \mathrm{and}\mathrm{predation}\end{array}}{\underbrace{{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t}}},$(2) ) and (3(3) $O_{t+1}=\underset{O\mathrm{gain}\mathrm{from}\mathrm{filtration}\mathrm{of}P}{\underbrace{P_t(1-{\mathrm{e}}^{-\delta O_t})}}+\underset{\begin{array}{c}O\mathrm{not}\mathrm{lost}\mathrm{tomortality},\\ \mathrm{harvest},\mathrm{disease}\mathrm{and}\mathrm{respiration}\end{array}}{\underbrace{{\mu }_O{O}_t}},$(3) ), our discrete-time NPO model is the following system of three difference equations. (NPO) $\{\begin{array}{rl}{N}_{t+1}& =\gamma {\mathrm{e}}^{-\alpha P_t}+{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t},\\ P_{t+1}& =N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}+{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t},\\ O_{t+1}& =P_t(1-{\mathrm{e}}^{-\delta O_t})+{\mu }_O{O}_t,\end{array}$(NPO) where $t=0,1,2,3,\dots$, with non-negative initial conditions $N_0\ge 0$, $P_0\ge 0$ and $O_0\ge 0$. NPO model predicts the vectors of population levels of nutrients, phytoplankton and oysters at time $(t+1)$, $(N_{t+1},P_{t+1},O_{t+1})$, from knowledge of their respective levels at time t, $(N_t,P_t,O_t)$, where $t=0,1,2,\dots$. The unit of time depends on the specific application. Baseline values for the discrete-time NPO model parameters are given in Table 2.

3. Boundedness of orbits

In this section, we establish the well-posedness of the discrete-time NPO model. That is, we show that the discrete-time NPO model has no unbounded growth for all non-negative initial conditions.

Let (4) $\{\begin{array}{rl}{N}_{t+1}& =F_1(N_t,P_t,O_t)=\gamma {\mathrm{e}}^{-\alpha P_t}+{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t},\\ P_{t+1}& =F_2(N_t,P_t,O_t)=N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}+{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t},\\ O_{t+1}& =F_3(N_t,P_t,O_t)=P_t(1-{\mathrm{e}}^{-\delta O_t})+{\mu }_O{O}_t,\end{array}$(4) for $t=0,1,2,3,\dots$. The set of density sequences generated by Model (4(4) $\{\begin{array}{rl}{N}_{t+1}& =F_1(N_t,P_t,O_t)=\gamma {\mathrm{e}}^{-\alpha P_t}+{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t},\\ P_{t+1}& =F_2(N_t,P_t,O_t)=N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}+{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t},\\ O_{t+1}& =F_3(N_t,P_t,O_t)=P_t(1-{\mathrm{e}}^{-\delta O_t})+{\mu }_O{O}_t,\end{array}$(4) ) is equivalent to the set of iterates of the map $F:[0,\mathrm{\infty })\times [0,\mathrm{\infty })\times [0,\mathrm{\infty })⟶[0,\mathrm{\infty })\times [0,\mathrm{\infty })\times [0,\mathrm{\infty })$defined by $F(N,P,O)=(F_1(N,P,O),F_2(N,P,O),F_3(N,P,O)).$

Lemma 3.1

In Model (4(4) $\{\begin{array}{rl}{N}_{t+1}& =F_1(N_t,P_t,O_t)=\gamma {\mathrm{e}}^{-\alpha P_t}+{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t},\\ P_{t+1}& =F_2(N_t,P_t,O_t)=N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}+{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t},\\ O_{t+1}& =F_3(N_t,P_t,O_t)=P_t(1-{\mathrm{e}}^{-\delta O_t})+{\mu }_O{O}_t,\end{array}$(4) ), if $N_0\ge 0$, $P_0\ge 0$ and $O_0\ge 0$, then $N_t\ge 0$, $P_t\ge 0$ and $O_t\ge 0$ for all $t\ge 0$. Furthermore, every non-negative initial condition has a bounded orbit in Model (4(4) $\{\begin{array}{rl}{N}_{t+1}& =F_1(N_t,P_t,O_t)=\gamma {\mathrm{e}}^{-\alpha P_t}+{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t},\\ P_{t+1}& =F_2(N_t,P_t,O_t)=N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}+{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t},\\ O_{t+1}& =F_3(N_t,P_t,O_t)=P_t(1-{\mathrm{e}}^{-\delta O_t})+{\mu }_O{O}_t,\end{array}$(4) ).

Proof.

Clearly, $N_t\ge 0$, $P_t\ge 0$ and $O_t\ge 0$ whenever $N_0\ge 0$, $P_0\ge 0$ and $O_0\ge 0$.

The inequalities $0\le N_{t+1}\le \gamma +{\mu }_N{N}_t$ for $t\ge 0$ imply, by induction argument, that $0\le N_t\le y_t$ where $y_t$ is the solution of the initial value problem $y_{t+1}=\gamma +{\mu }_N{y}_t,y_0=N_0.$Because $0\le {\mu }_N<1$, this linear difference equation has a globally attracting equilibrium $y^{\ast }=\frac{\gamma }{1-{\mu }_N}$, i.e. $\underset{t\to \mathrm{\infty }}{lim}{y}_t=y^{\ast }$. Thus, $0\le \underset{t\to \mathrm{\infty }}{limsup}{N}_t\le \underset{t\to \mathrm{\infty }}{lim}{y}_t=y^{\ast }$which implies $N_t$ is bounded. Let $N^{\ast }$ be a bound for $N_t$, i.e. $0\le N_t\le N^{\ast }$ for all $t\le 0$.

From $0\le P_{t+1}\le N_t+{\mu }_P{P}_t\le N^{\ast }+{\mu }_P{P}_t$ for $t\ge 0$, we conclude, by an analogous argument, that $P_t$ is bounded.

Finally, by an analogous argument used for $P_t$, one can argue that $O_t$ is bounded.

Hence, the orbit of every non-negative initial condition is bounded in Model (4(4) $\{\begin{array}{rl}{N}_{t+1}& =F_1(N_t,P_t,O_t)=\gamma {\mathrm{e}}^{-\alpha P_t}+{\mu }_N{N}_t{\mathrm{e}}^{-\alpha P_t},\\ P_{t+1}& =F_2(N_t,P_t,O_t)=N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}+{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t},\\ O_{t+1}& =F_3(N_t,P_t,O_t)=P_t(1-{\mathrm{e}}^{-\delta O_t})+{\mu }_O{O}_t,\end{array}$(4) ).

4. Nutrients-only equilibrium point: stability and sensitivity analysis

Next, we compute the nutrients-only equilibrium point of the NPO model. Then, we establish conditions for its stability and perform local and global sensitivity analysis.

Model NPO has a nutrients-only equilibrium point $E_N=(\frac{\gamma }{1-{\mu }_N},0,0).$

4.1. $E_N$: stability

To study the stability of $E_N$, we introduce the threshold parameter ${\mathcal{R}}_N=\frac{\alpha \gamma }{(1-{\mu }_P)(1-{\mu }_N)}.$Note that ${\mathcal{R}}_N>0$. ${\mathcal{R}}_N<1$ is equivalent to $\gamma <\frac{(1-{\mu }_N)(1-{\mu }_P)}{\alpha }$ and ${\mathcal{R}}_N>1$ is equivalent to $\gamma >\frac{(1-{\mu }_N)(1-{\mu }_P)}{\alpha }$.

When γ is sufficiently ‘small’, then ${\mathcal{R}}_N<1$, and $E_N$ is locally asymptotically stable. However, when γ is ‘large’, then ${\mathcal{R}}_N>1$, and $E_N$ is unstable. Using baseline parameter values of Table 2, ${\mathcal{R}}_N=16.7$. When ${\mathcal{R}}_N>1$, the Model NPO has an only oysters-free equilibrium point. We summarize this stability result in the following theorem.

Theorem 4.1

In Model NPO, the phytoplankton-, oysters-free boundary fixed point, $E_N$, is locally asymptotically stable and small initial nutrients, phytoplankton and oysters levels lead to extinction of phytoplankton and oysters when ${\mathcal{R}}_N<1$. However, $E_N$ is unstable when ${\mathcal{R}}_N>1$.

Proof.

The Jacobian matrix of the NPO model evaluated at $E_N$ is ${\mathcal{J}}_{|E_N}=\left[\begin{array}{ccc}{\mu }_N& -{\displaystyle \frac{\alpha \gamma }{1-{\mu }_N}}& 0\\ 0& {\displaystyle \frac{\alpha \gamma }{1-{\mu }_N}}+{\mu }_P& 0\\ 0& 0& {\mu }_0\end{array}\right].$The eigenvalues of ${\mathcal{J}}_{|E_N}$ are: ${\mu }_N$, $\frac{\alpha \gamma }{1-{\mu }_N}+{\mu }_P$ and ${\mu }_0$. Since $0<{\mu }_N<1$ and $0<{\mu }_0<1$, the stability of $E_N$ is determined by $\frac{\alpha \gamma }{1-{\mu }_N}+{\mu }_P$. ${\mathcal{R}}_N<1⟺\frac{\alpha \gamma }{(1-{\mu }_N)}+{\mu }_P<1\mathrm{and}{\mathcal{R}}_N>1⟺\frac{\alpha \gamma }{(1-{\mu }_N)}+{\mu }_P>1.$Hence, ${\mathcal{R}}_N<1$ implies $E_N$ is locally asymptotically stable and ${\mathcal{R}}_N>1$ implies $E_N$ is unstable.

Corollary 4.2

In Model NPO, the phytoplankton-, oysters-free boundary equilibrium point, $E_N$, is globally asymptotically stable and independently of positive initial population numbers, both phytoplankton and oysters go extinct when ${\mathcal{R}}_N<1$.

Proof.

Since $1-{\mathrm{e}}^{-\alpha P_t}\le \alpha P_t$ and ${\mathrm{e}}^{-\delta O_t}\le 1$, $\begin{array}{rl}{P}_{t+1}& =N_t(1-{\mathrm{e}}^{-\alpha P_t}){\mathrm{e}}^{-\delta O_t}+{\mu }_P{P}_t{\mathrm{e}}^{-\delta O_t}\\ & \le \alpha N_t{P}_t+{\mu }_P{P}_t\\ & \le (\alpha N_t+{\mu }_P)P_t.\end{array}$Since $\underset{t\to \mathrm{\infty }}{limsup}{N}_t\le y^{\ast }=\frac{\gamma }{1-{\mu }_N}$this implies, for arbitrary $\epsilon >0$, there exists a time $T(\epsilon )>0$ such that $0\le N_t\le y^{\ast }+\epsilon$ for all $t\ge T(\epsilon )$. Since ${\mathcal{R}}_N<1$ it follows that $\alpha \frac{\gamma }{1-{\mu }_N}+{\mu }_P<1$and we can choose ε so small that (5) $\alpha \frac{\gamma }{1-{\mu }_N}+{\mu }_P+\alpha \epsilon <1.$(5) Since $0\le P_{t+1}\le (\alpha N_t+{\mu }_P)P_t$which for $t\ge T(\epsilon )$ implies $\begin{array}{rl}0& \le P_{t+1}\le (\alpha y^{\ast }+\alpha \epsilon +{\mu }_P)P_t\\ & =(\alpha \frac{\gamma }{1-{\mu }_N}+{\mu }_P+\alpha \epsilon )P_t\end{array}$and, by (5(5) $\alpha \frac{\gamma }{1-{\mu }_N}+{\mu }_P+\alpha \epsilon <1.$(5) ), that $\underset{t\to \mathrm{\infty }}{lim}{P}_t=0$.

Since ${\displaystyle \underset{t\to \mathrm{\infty }}{lim}{P}_t=0}$, then $\mathrm{\exists }{k}_1>0$ such that $P_t<\frac{1-{\mu }_0}{\delta }$, $\mathrm{\forall }t\ge k_1$. Hence, ${\displaystyle \underset{t\to \mathrm{\infty }}{lim}{O}_t=0}$.

By Theorem 4.1, $E_N$ is a locally asymptotically stable boundary equilibrium point of Model NPO. Hence, $E_N$ is globally stable when ${\mathcal{R}}_N<1$.

4.2. $E_N$: sensitivity analysis

Next, we use the local sensitivity analysis index to estimate the relative change in a state variable caused when the value of a specific model parameter changes. Consequently, as in [7, 23–26], we use the following definition of normalized forward sensitivity index to carry out the local sensitivity analysis and compute normalized sensitivity indices.

Definition 4.3

[7]

The normalized forward sensitivity index of a variable, u, that depends differentiably on a parameter, q, is defined as: ${\mathrm{{\rm Y}}}_q^u:=\frac{\mathrm{\partial }u}{\mathrm{\partial }q}\times \frac{q}{u}.$

We summarize the computations of the normalized sensitivity indices of equilibrium nutrients level at $E_N$ in Table 3.

Table 3. Normalized sensitivity indices of equilibrium nutrients level at $E_N$ to parameters evaluated at baseline parameter values of Table 2 and their order of importance.

Parameter	Sensitivity index	Order of importance
γ	1	2
${\mu }_N$	2.3333	1
${\mu }_P$, α, ${\mu }_O$, δ	0	3

Sensitivity indices of Table 3 indicate that increasing (respectively, decreasing) nutrients flow, γ, by $1\mathrm{\%}$ will increase (respectively, decrease) equilibrium nutrients level at $E_N$ by $1\mathrm{\%}$. Increasing (respectively, decreasing) escapement rate of nutrients, ${\mu }_N$, will increase (respectively, decrease) nutrients equilibrium level at $E_N$, by $2.3333\mathrm{\%}$.

Sensitivity indices of Table 3 suggest that nutrients level in the nutrients-only equilibrium has two important parameters: (1) flow rate of nutrients, γ and (2) escapement rate of nutrients, ${\mu }_N$. This is consistent with the sensitivity analysis results in [25] for the continuous-time NPO model that nutrients level in the phytoplankton- and oyster-free only equilibrium also has two important parameters: (1) nutrients flow and (2) nutrients loss [25].

Now, we use the global sensitivity analysis at $E_N$ to measure the effect generated by one specific model parameter on a model variable, while all the other model parameters vary [5, 6, 18]. Similar to the global sensitivity approach used in [25], we adopt the ANOVA (ANalysis Of VAriance) decomposition-based method called Sobol' sensitivity method [22]. As in [25], we perform the global sensitivity analysis using Sobol' method with the free software tool, Global Sensitivity Analysis Toolbox (GSAT) [4]. GSAT is performed using the Matlab software. We compute Sobol' sensitivity indices of type first-order and total-order. The samplings of our model parameter values are based on quasi-random sequences and are used to run 200,000 simulations. The ranges of model parameter are given in Table 2.

Sobol' sensitivity results of Table 4 indicate that the escapement rate of nutrients, ${\mu }_N$, is the most important parameter and contributes about $72\mathrm{\%}$ of the nutrients equilibrium level at $E_N$ total variance. There is no large difference between first and total Sobol' index for each parameter. Therefore, there is no significant interaction between the parameters. However, phytoplankton biomass and oysters biomass at the fixed point $E_N$ are not sensitive to changes in parameter values, and therefore remain fixed at zero. This is in contrast to the Sobol' sensitivity results for the continuous-time NPO model [25] where nutrients flow and nutrients loss are the most important parameters and contribute respectively about $24\mathrm{\%}$ and $57\mathrm{\%}$ to the total variance of the nutrients level in the nutrients-only equilibrium.

4.3. ${\mathcal{R}}_N$: sensitivity analysis

Now, we use sensitivity analysis to study the impact of parameter changes on ${\mathcal{R}}_N$.

Table 4. First-order and total-order Sobol sensitivity indices of equilibrium nutrients level at $E_N$, related to changes in parameters with ranges listed in Table 2 and their order of importance.

Parameter	Sobol' first-order sensitivity	Order of importance	Sobol' total-order sensitivity	Order of importance
γ	0	2	0.2800	2
${\mu }_N$	0.7200	1	0.9999	1
${\mu }_P$, α, ${\mu }_O$, δ	0	3	0	3

Sensitivity indices of Table 5 show that increasing (respectively, decreasing) the escapement rate of nutrients, ${\mu }_N$, by $1\mathrm{\%}$ will increase (respectively, decrease) ${\mathcal{R}}_N$ by $2.3333\mathrm{\%}$. Increasing (respectively, decreasing) survival rate of phytoplankton, ${\mu }_P$, will increase (respectively, decrease) ${\mathcal{R}}_N$, by $1.5\mathrm{\%}$.

Note that, because the list of parameters for the continuous- and discrete-time NPO models are not identical, the threshold parameter ${\mathcal{R}}_N$ is different for the two models [25]. However, in both models when ${\mathcal{R}}_N<1$, nutrients persist while phytoplankton and oysters population go extinct. However, when ${\mathcal{R}}_N>1$ both nutrients and phytoplankton persist with or without oysters.

Table 5. Normalized sensitivity indices of ${\mathcal{R}}_N$, to parameters evaluated at baseline of Table 2 and their order of importance.

Parameter	Sensitivity index	Order of importance
γ, α	1	3
${\mu }_N$	2.3333	1
${\mu }_P$	1.5	2
${\mu }_O$, δ	0	4

Sensitivity indices of Table 5 show that the threshold parameter ${\mathcal{R}}_N$ has two important parameters: (1) escapement rate of nutrients, ${\mu }_N$, (2) survival rate of phytoplankton, ${\mu }_P$. While sensitivity analysis of in [25] for the continuous-time NPO model show that ${\mathcal{R}}_N$ has two important parameters: (1) removal of phytoplankton by sinking and predation and (2) the maximum nutrients uptake by phytoplankton. In Table 6, we summarize our results on the global sensitivity analysis of ${\mathcal{R}}_N$.

Table 6. First-order and total-order Sobol sensitivity indices ${\mathcal{R}}_N$ related to changes in parameters with ranges listed in Table 2 and their order of importance.

Parameter	Sobol' first-order sensitivity	Order of importance	Sobol' total-order sensitivity	Order of importance
γ	0	3	0.4282	3
${\mu }_N$	0.2177	1	0.5256	2
${\mu }_P$	0.1759	2	0.8972	1
α	0	4	0.0291	4
${\mu }_O$	0	5	0	5
δ	0	6	0	6

Sobol' sensitivity results of Table 6 show that escapement rate of nutrients, ${\mu }_N$, contributes $21.77\mathrm{\%}$ and the survival rate of phytoplankton, ${\mu }_p$, contributes $17.59\mathrm{\%}$ of variance of ${\mathcal{R}}_N$. Also, there is a large difference between first and total Sobol' indices for nutrients flow, γ, nutrients escapement rate, ${\mu }_N$, phytoplankton survival rate, ${\mu }_P$ and nutrients uptake by phytoplankton, α. Hence, each of these parameters interact strongly with the others.

Sobol' sensitivity results of Table 6 illustrate that ${\mathcal{R}}_N$ has two important parameters: (1) escapement rate of nutrients and (2) survival rate of phytoplankton. While Sobol' sensitivity results in [25] for the continuous-time NPO model showed that ${\mathcal{R}}_N$ has two important parameters: (1) removal of phytoplankton by sinking and predation and (2) the maximum nutrients uptake by phytoplankton.

4.4. Illustrative example: $E_N$ and ${\mathcal{R}}_N$

Now, we use the following specific example to illustrate the stability of the nutrients-only equilibrium point, $E_N$, in Model NPO.

Example 4.1

In Model NPO, let $\alpha =0.2$, $\gamma =2$, ${\mu }_N=0.2$, ${\mu }_P=0.2$, ${\mu }_O=0.09$ and $\delta =0.5$.

In Example 4.1, ${\mathcal{R}}_N=0.625<1$. As predicted, by Corollary 4.2, Figure 2 shows that $E_N=(2.5,0,0)$ is globally asymptotically stable.

To study the impact of high nutrients inflow on Figure 2, we increase γ from $\gamma =2$ to $\gamma =4$ and keep all the other parameters fixed at their current values in Example 4.1. Then ${\mathcal{R}}_N$ increases from ${\mathcal{R}}_N=0.625$ to ${\mathcal{R}}_N=1.25>1$. Figure 3 shows that $E_N=(5,0,0)$ is unstable and the initial condition $(5,0.001,0)$ limits on an only oysters-free equilibrium point $(4.2642,0.6466,0)$. Thus, high nutrients flow promotes the growth of phytoplankton population.

Figure 2. ${\mathcal{R}}_N=0.625<1$ and Model NPO has a globally asymptotically stable nutrients-only equilibrium point, $E_N=(2.5,0,0)$ where $\alpha =0.2$, $\gamma =2$, ${\mu }_N=0.2$, ${\mu }_P=0.2$, ${\mu }_O=0.09$ and $\delta =0.5$.

5. Oysters-free equilibrium point: stability

By Corollary 4.2, ${\mathcal{R}}_N<1$ implies $E_N$ is globally asymptotically stable and Model NPO has no equilibrium points with positive levels of phytoplankton and oysters. To ensure persistence of phytoplankton, we only consider ${\mathcal{R}}_N>1$. In this case, Model NPO has an only oysters-free equilibrium point with positive level of phytoplankton, $E_{NP}=(\overline{N},\overline{P},0)$; where $\overline{N}$ and $\overline{P}$ are positive the solutions of (6) $\{\begin{array}{l}\overline{N}(1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}})=\gamma {\mathrm{e}}^{-\alpha \overline{P}},\\ \overline{N}(1-{\mathrm{e}}^{-\alpha \overline{P}})=(1-{\mu }_P)\overline{P}.\end{array}$(6) That is, at $(\overline{N},\overline{P},0)$, (7) $\{\begin{array}{rl}\overline{N}& =\frac{\gamma {\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}},\\ \frac{{\mathrm{e}}^{\alpha \overline{P}}-1}{{\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N}& =\frac{(1-{\mu }_P)}{\gamma }\overline{P}{\mathrm{e}}^{\alpha \overline{P}},\end{array}$(7) where $\overline{N}$ and $\overline{P}$ are positive.

Figure 3. ${\mathcal{R}}_N=1.25>1$ and Model NPO has an unstable nutrients-only equilibrium point, $E_N=(5,0,0)$ and an asymptotically stable only oysters-free equilibrium $(4.2642,0.6466,0)$, where $\alpha =0.2$, $\gamma =4$, ${\mu }_N=0.2$, ${\mu }_P=0.2$, ${\mu }_O=0.09$ and $\delta =0.5$.

Let $f_1(P)=\frac{{\mathrm{e}}^{\alpha \overline{P}}-1}{{\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N}\text{and}{f}_2(P)=\frac{(1-{\mu }_P)}{\gamma }\overline{P}{\mathrm{e}}^{\alpha \overline{P}}.$Notice that $f_1(0)=f_2(0)=0$. Clearly, $f_1^{\prime }(P)>0$ and $f_1^{″}(P)<0$ for all $P\ge 0$. Hence, the graph of $f_1$ is increasing and concave down in $[0,\mathrm{\infty })$. Also, $f_2^{\prime }(P)>0$ and $f_2^{″}(P)>0$ for all $P\ge 0$. Hence, the graph of $f_2$ is increasing and concave up in $[0,\mathrm{\infty })$. Since, $f_1^{\prime }(0)=\frac{\alpha }{1-{\mu }_N}$ and $f_2^{\prime }(0)=\frac{1-{\mu }_P}{\gamma }$, ${\mathcal{R}}_N>1$ implies $f_2^{\prime }(0)<f_1^{\prime }(0)$. Since ${\displaystyle \underset{P\to \mathrm{\infty }}{lim}{f}_1(P)=1}$ and ${\displaystyle \underset{P\to \mathrm{\infty }}{lim}{f}_2(P)=\mathrm{\infty }}$, using the monotonicity and concavity of the graphs of $f_1$ and $f_2$, we obtain that the graphs of $f_1$ and $f_2$ intersect at exactly one point, $\overline{P}$, in $(0,\mathrm{\infty })$, whenever ${\mathcal{R}}_N>1$ (see Figure 4). That is ${\mathcal{R}}_N>1$ implies Model NPO has a unique oysters-free equilibrium point, $E_{NP}$, with positive level of nutrients and phytoplankton. We summarize this in the following result.

Figure 4. The graphs of $f_1(P)$ and $f_2(P)$ intersect at P = 0 and at $\overline{P}=1.4279$ where $\alpha =0.8$, $\gamma =4$, ${\mu }_N=0.75$, ${\mu }_P=0.2$, ${\mu }_O=0.09$, $\delta =0.5$.

Theorem 5.1

When ${\mathcal{R}}_N>1$, Model NPO has a unique oysters-free equilibrium point, $E_{NP}=(\overline{N},\overline{P},0)$, where $\overline{N}$ and $\overline{P}$ are the unique positive solutions of (6(6) $\{\begin{array}{l}\overline{N}(1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}})=\gamma {\mathrm{e}}^{-\alpha \overline{P}},\\ \overline{N}(1-{\mathrm{e}}^{-\alpha \overline{P}})=(1-{\mu }_P)\overline{P}.\end{array}$(6) ) and (7(7) $\{\begin{array}{rl}\overline{N}& =\frac{\gamma {\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}},\\ \frac{{\mathrm{e}}^{\alpha \overline{P}}-1}{{\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N}& =\frac{(1-{\mu }_P)}{\gamma }\overline{P}{\mathrm{e}}^{\alpha \overline{P}},\end{array}$(7) ).

5.1. $E_{NP}$: stability

In this section, the asymptotic stability of the oysters-free equilibrium $E_{NP}$ will be analysed. We claim the following result:

Theorem 5.2

Let ${\mathcal{R}}_N>1$ and $E_{NP}=(\overline{N},\overline{P},0)$, where $\overline{N}$ and $\overline{P}$ are solutions of (6(6) $\{\begin{array}{l}\overline{N}(1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}})=\gamma {\mathrm{e}}^{-\alpha \overline{P}},\\ \overline{N}(1-{\mathrm{e}}^{-\alpha \overline{P}})=(1-{\mu }_P)\overline{P}.\end{array}$(6) ) and (7(7) $\{\begin{array}{rl}\overline{N}& =\frac{\gamma {\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}},\\ \frac{{\mathrm{e}}^{\alpha \overline{P}}-1}{{\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N}& =\frac{(1-{\mu }_P)}{\gamma }\overline{P}{\mathrm{e}}^{\alpha \overline{P}},\end{array}$(7) ), and the following two inequalities are satisfied.

1. $0<\overline{P}<\frac{1-{\mu }_0}{\delta }$,

2. $\begin{array}{rl}& {\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}+\frac{\alpha \gamma {\mathrm{e}}^{-2\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}+{\mu }_P\\ & <1+\frac{\alpha \gamma {\mu }_N{\mathrm{e}}^{-3\alpha \overline{P}}-({\mu }_P{\mu }_N^2+\alpha \gamma ){\mathrm{e}}^{-2\alpha \overline{P}}+({\mu }_P{\mu }_N+\alpha \gamma ){\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}<2.\end{array}$

Then, the oysters-free equilibrium point, $E_{NP}$, is asymptotically stable and small initial nutrients, phytoplankton and oysters levels lead to extinction of oysters. However, $E_{NP}$ is unstable and persistence of nutrients, phytoplankton and oysters occur, when either inequality (i) or (ii) is not satisfied.

Proof.

To study stability of the oysters-free equilibrium point, $E_{NP}$, we compute the eigenvalues of the Jacobian matrix of the model evaluated $E_{NP}=(\overline{N},\overline{P},0)$, where $\overline{N}$ and $\overline{P}$ are solutions of (6(6) $\{\begin{array}{l}\overline{N}(1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}})=\gamma {\mathrm{e}}^{-\alpha \overline{P}},\\ \overline{N}(1-{\mathrm{e}}^{-\alpha \overline{P}})=(1-{\mu }_P)\overline{P}.\end{array}$(6) ) and (7(7) $\{\begin{array}{rl}\overline{N}& =\frac{\gamma {\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}},\\ \frac{{\mathrm{e}}^{\alpha \overline{P}}-1}{{\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N}& =\frac{(1-{\mu }_P)}{\gamma }\overline{P}{\mathrm{e}}^{\alpha \overline{P}},\end{array}$(7) ). The Jacobian matrix of NPO model evaluated at $E_{NP}$ is ${\mathcal{J}}_{|E_{NP}}=\left[\begin{array}{ccc}{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}& {\displaystyle \frac{-\alpha \gamma {\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}}& 0\\ 1-{\mathrm{e}}^{-\alpha \overline{P}}& {\displaystyle \frac{\alpha \gamma {\mathrm{e}}^{-2\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}}+{\mu }_P& -{\displaystyle \frac{\delta \gamma {\mathrm{e}}^{-\alpha \overline{P}}(1-{\mathrm{e}}^{-\alpha \overline{P}})}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}}+{\mu }_P\overline{P}\\ 0& 0& \delta \overline{P}+{\mu }_0\end{array}\right].$Hence, ${\mathcal{J}}_{|E_{NP}}$ has three eigenvalues ${\lambda }_1$, ${\lambda }_2$ and ${\lambda }_3$, where

• ${\lambda }_1$ and ${\lambda }_2$ are the eigenvalues of the matrix $A=\left[\begin{array}{cc}{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}& {\displaystyle \frac{-\alpha \gamma {\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}}\\ 1-{\mathrm{e}}^{-\alpha \overline{P}}& {\displaystyle \frac{\alpha \gamma {\mathrm{e}}^{-2\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}}+{\mu }_P\end{array}\right].$

• ${\lambda }_3=\delta \overline{P}+{\mu }_0$.

Since the model parameters are considered to be positive, then $|{\lambda }_3|<1⟺0<\overline{P}<\frac{1-{\mu }_0}{\delta }$which means that $|{\lambda }_3|<1$ is equivalent to the condition (i). By the Jury stability criteria [11, 17], $\begin{array}{rl}& \left|{\lambda }_1\right|<1\mathrm{and}\left|{\lambda }_2\right|<1⟺|\mathrm{tr}(A)|<1+\det (A)<2\\ & ⟺|{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}+\frac{\alpha \gamma {\mathrm{e}}^{-2\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}+{\mu }_P|\\ & <1+\frac{\alpha \gamma {\mu }_N{\mathrm{e}}^{-3\alpha \overline{P}}-({\mu }_P{\mu }_N^2+\alpha \gamma ){\mathrm{e}}^{-2\alpha \overline{P}}+({\mu }_P{\mu }_N+\alpha \gamma ){\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}<2,\end{array}$where $tr(A)$ is the trace of A and $det(A)$ is determinant of A. $0<{\mu }_N<1$, then $\begin{array}{rl}\left|{\lambda }_1\right|& <1\mathrm{and}\left|{\lambda }_2\right|<1⟺{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}+\frac{\alpha \gamma {\mathrm{e}}^{-2\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}+{\mu }_P\\ & <1+\frac{\alpha \gamma {\mu }_N{\mathrm{e}}^{-3\alpha \overline{P}}-({\mu }_P{\mu }_N^2+\alpha \gamma ){\mathrm{e}}^{-2\alpha \overline{P}}+({\mu }_P{\mu }_N+\alpha \gamma ){\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}<2,\end{array}$Which means that $|{\lambda }_1|<1\mathrm{and}|{\lambda }_2|<1$ are equivalent to condition (ii). When $|{\lambda }_1|<1|{\lambda }_2|<1\mathrm{and}|{\lambda }_3|<1$ then $E_{NP}$ is asymptotically stable. $E_{NP}$ is unstable when either of the conditions is not satisfied. Hence the result of the Theorem.

5.2. Illustrative example: $E_{NP}$ and ${\mathcal{R}}_N$

Now, we use the following specific example to illustrate the stability of the oysters-free equilibrium point, $E_{NP}$, in Model NPO.

Example 5.1

In Model NPO, let $\alpha =0.2$, $\gamma =2$, ${\mu }_N=0.9$, ${\mu }_P=0.7$, ${\mu }_O=0.09$ and $\delta =0.2$.

In Example 5.1, ${\mathcal{R}}_N=13.33>1$ and the Jacobian matrix of NPO model evaluated at $E_{NP}$, ${\mathcal{J}}_{|E_{NP}}$ has three eigenvalues, ${\lambda }_1=0.4811+0.4364i$, ${\lambda }_2=0.4811-0.4364i$ and ${\lambda }_3=0.7250$. Note that, $|{\lambda }_1|=|{\lambda }_2|=0.6495<1$ and $|{\lambda }_3|<1$. So conditions (i) and (ii) of Theorem 5.2 are satisfied. As predicted by Theorem 5.2, Figure 5 shows that the only oysters-free equilibrium point, $E_{NP}=(2.0263,3.1750,0)$, is locally asymptotically stable.

To study the impact of high nutrients flow on Figure 5, we increase γ from $\gamma =2$ to $\gamma =4$ and keep all the other parameters fixed at their current values in Example 5.1. Then ${\mathcal{R}}_N$ increases form ${\mathcal{R}}_N=13.33$ to ${\mathcal{R}}_N=26.67>1$. The eigenvalues of the Jacobian matrix evaluated at $E_{NP}$, ${\mathcal{J}}_{|E_{NP}}$, are ${\lambda }_1=0.4329+0.5295i$, ${\lambda }_2=0.4329-0.5295i$ and ${\lambda }_3=1.0507$. Note that $|{\lambda }_3|>1$. So, condition (i) of Theorem 5.2 is not satisfied. As predicted by Theorem 5.2, Figure 6 shows that $E_{NP}=(2.3342,4.8037,0)$ is unstable and the initial condition $(2.3342,4.8037,0.0001)$ limits at the interior equilibrium point $(2.4817,4.6050,0.1204)$. Thus, high nutrients flow promotes phytoplankton growth which serves as food for the survival of the oysters.

Figure 5. Model NPO has the locally asymptotically stable and the oysters-free equilibrium point, $E_{NP}=(2.0263,3.1750,0)$, where $\alpha =0.2$, $\gamma =2$, ${\mu }_N=0.9$, ${\mu }_P=0.7$, ${\mu }_O=0.09$ and $\delta =0.2$.

6. Interior equilibrium point: stability

To find the interior equilibrium point of Model NPO, $E_{NPO}$, we solve (8) $\{\begin{array}{l}\gamma {\mathrm{e}}^{-\alpha P}=(1-{\mu }_N{\mathrm{e}}^{-\alpha P})N\\ N(1-{\mathrm{e}}^{-\alpha P}){\mathrm{e}}^{-\delta O}=(1-{\mu }_P{\mathrm{e}}^{-\delta O})P\\ P(1-{\mathrm{e}}^{-\delta O})=(1-{\mu }_O)O\end{array}$(8) for $(N,P,O)\in (0,\mathrm{\infty })\times (0,\mathrm{\infty })\times (0,\mathrm{\infty })$.

Figure 6. ${\mathcal{R}}_N=26,67>1$ and Model NPO has an unstable oysters-free equilibrium point, $E_{NP}=(2.3342,4.8037,0)$ and a stable interior equilibrium point at $E_{NPO}=(2.4817,4.6050,0.1204)$, where $\alpha =0.2$, $\gamma =4$, ${\mu }_N=0.9$, ${\mu }_P=0.7$, ${\mu }_O=0.09$ and $\delta =0.2$.

Then $N=\frac{\gamma {\mathrm{e}}^{-\alpha P}}{1-{\mu }_N{\mathrm{e}}^{-\alpha P}}\mathrm{and}\frac{\gamma {\mathrm{e}}^{-\alpha P}(1-{\mathrm{e}}^{-\alpha P})}{(1-{\mu }_N{\mathrm{e}}^{-\alpha P})}{\mathrm{e}}^{-\delta O}=(1-{\mu }_P{\mathrm{e}}^{-\delta O})P.$Thus, ${\mathrm{e}}^{-\delta O}=\frac{P({\mathrm{e}}^{\alpha P}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha P})+{\mu }_{P}P({\mathrm{e}}^{\alpha P}-{\mu }_N)}.$That is, $O=-\frac{1}{\delta }\ln \left(\frac{P({\mathrm{e}}^{\alpha P}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha P})+{\mu }_{P}P({\mathrm{e}}^{\alpha P}-{\mu }_N)}\right).$By substituting for O in the third equation of (8(8) $\{\begin{array}{l}\gamma {\mathrm{e}}^{-\alpha P}=(1-{\mu }_N{\mathrm{e}}^{-\alpha P})N\\ N(1-{\mathrm{e}}^{-\alpha P}){\mathrm{e}}^{-\delta O}=(1-{\mu }_P{\mathrm{e}}^{-\delta O})P\\ P(1-{\mathrm{e}}^{-\delta O})=(1-{\mu }_O)O\end{array}$(8) ), we obtain $\frac{\delta }{(1-{\mu }_O)}P=-\frac{\ln \left(\frac{P({\mathrm{e}}^{\alpha P}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha P})+{\mu }_{P}P({\mathrm{e}}^{\alpha P}-{\mu }_N)}\right)}{1-\frac{P({\mathrm{e}}^{\alpha P}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha P})+{\mu }_{P}P({\mathrm{e}}^{\alpha P}-{\mu }_N)}}.$Let $f_1(P)=\frac{\delta }{(1-{\mu }_O)}P\mathrm{and}{f}_2(P)=-\frac{\ln \left(\frac{P({\mathrm{e}}^{\alpha P}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha P})+{\mu }_{P}P({\mathrm{e}}^{\alpha P}-{\mu }_N)}\right)}{1-\frac{P({\mathrm{e}}^{\alpha P}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha P})+{\mu }_{P}P({\mathrm{e}}^{\alpha P}-{\mu }_N)}}.$Notice that $f_1(0)=0$ and $f_1^{\prime }(P)>0$ for all P>0 and ${\displaystyle \underset{P\to \mathrm{\infty }}{lim}{f}_1(P)=\mathrm{\infty }}$. When ${\mathcal{R}}_N>1$, then ${\displaystyle \underset{P\to 0^{+}}{lim}{f}_2(P)=R_N(1-{\mu }_p)+{\mu }_p>0,}$${\displaystyle \underset{P\to \mathrm{\infty }}{lim}{f}_2(P)=0}$ and $f_2^{\prime }(P)<0$ for all $P\ge 0$. Consequently the graphs of $f_1$ and $f_2$ intersect at exactly one point, $\overline{P}$, in $(0,\mathrm{\infty })$ whenever ${\mathcal{R}}_N>1$ (see Figure 7). That is ${\mathcal{R}}_N>1$ implies Model NPO has a unique interior equilibrium point $E_{NPO}=(\overline{N},\overline{P},\overline{O})$ where (9) $\{\begin{array}{rl}\overline{N}& =\frac{\gamma {\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}\\ \overline{P}& =\frac{(1-{\mu }_O)}{\delta }\frac{\ln \left(\frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}\right)}{\frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}-1}\\ \overline{O}& =-\frac{1}{\delta }\ln \left(\frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}\right).\end{array}$(9) We summarize this in the following results.

Figure 7. When $\alpha =0.2$, $\gamma =4$, ${\mu }_N=0.9$, ${\mu }_P=0.7$, ${\mu }_O=0.09$ and $\delta =0.2$ then, the functions $f_1(P)=\frac{\delta }{(1-{\mu }_O)}P$ and $f_2(P)=-\frac{\ln (\frac{P({\mathrm{e}}^{\alpha P}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha P})+{\mu }_{P}P({\mathrm{e}}^{\alpha P}-{\mu }_N)})}{1-\frac{P({\mathrm{e}}^{\alpha P}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha P})+{\mu }_{P}P({\mathrm{e}}^{\alpha P}-{\mu }_N)}}$ intersect at P = 4.6189. Hence, Model NPO has an interior equilibrium point $E_{NPO}=(2.4817,4.6050,0.1204)$.

Theorem 6.1

When ${\mathcal{R}}_N>1$, Model NPO has a unique interior equilibrium point, $E_{NPO}=(\overline{N},\overline{P},\overline{O})$, with positive levels of nutrients, phytoplankton and oysters, where $(\overline{N},\overline{P},\overline{O})$ is the solution of (9(9) $\{\begin{array}{rl}\overline{N}& =\frac{\gamma {\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}\\ \overline{P}& =\frac{(1-{\mu }_O)}{\delta }\frac{\ln \left(\frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}\right)}{\frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}-1}\\ \overline{O}& =-\frac{1}{\delta }\ln \left(\frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}\right).\end{array}$(9) ) in $(0,\mathrm{\infty })\times (0,\mathrm{\infty })\times (0,\mathrm{\infty })$.

6.1. Illustrative example: $E_{NPO}$ and Neimark–Sacker bifurcation

An analysis of the local stability of the interior equilibrium point, $E_{NPO}=(\overline{N},\overline{P},\overline{O})$ is not tractable, so we show, by numerical examples, that $E_{NPO}$ can be either stable or unstable.

We evaluate the Jacobian matrix of the model evaluated at $E_{NPO}$, ${\mathcal{J}}_{|E_{NPO}}$. $\begin{array}{rl}{\mathcal{J}}_{|E_{NPO}}& =[\begin{array}{cc}{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}& -\alpha \gamma {\mathrm{e}}^{-\alpha \overline{P}}-{\displaystyle \frac{\alpha {\mu }_N\gamma {\mathrm{e}}^{-2\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}}\\ {\displaystyle \frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)(1-{\mathrm{e}}^{-\alpha \overline{P}})}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}}& {\displaystyle \frac{(\frac{\alpha \gamma {\mathrm{e}}^{-2\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}+{\mu }_P)\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}}\\ 0& 1-{\displaystyle \frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}}\end{array}\\ & \begin{array}{c}0\\ {\displaystyle \frac{(\frac{-\delta \gamma {\mathrm{e}}^{-\alpha \overline{P}}(1-{\mathrm{e}}^{-\alpha \overline{P}})}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}-\delta {\mu }_P\overline{P})\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}}\\ {\displaystyle \frac{\delta {\overline{P}}^2({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}}+{\mu }_O\end{array}]\end{array}$where $\overline{N}$, $\overline{P}$ and $\overline{O}$ are solutions of (9(9) $\{\begin{array}{rl}\overline{N}& =\frac{\gamma {\mathrm{e}}^{-\alpha \overline{P}}}{1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}}}\\ \overline{P}& =\frac{(1-{\mu }_O)}{\delta }\frac{\ln \left(\frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}\right)}{\frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}-1}\\ \overline{O}& =-\frac{1}{\delta }\ln \left(\frac{\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}{\gamma (1-{\mathrm{e}}^{-\alpha \overline{P}})+{\mu }_P\overline{P}({\mathrm{e}}^{\alpha \overline{P}}-{\mu }_N)}\right).\end{array}$(9) ).

$E_{NPO}$ is locally asymptotically stable if $|\lambda |<1$ for all the eigenvalues, λ, of ${\mathcal{J}}_{|E_{NPO}}$. Otherwise, $E_{NPO}$ is unstable.

Now, we use the following two specific examples to illustrate the stability of the interior equilibrium point, $E_{NPO}$, of Model NPO.

Example 6.1

In Model NPO, let $\alpha =1.2$, $\gamma =3$, ${\mu }_N=0.2$, ${\mu }_P=0.6$, ${\mu }_O=0.09$ and $\delta =0.9$.

In Example 6.1, ${\mathcal{R}}_N=11.25>1$ and the Jacobian matrix of NPO model evaluated at $E_{NPO}$, ${\mathcal{J}}_{|E_{NPO}}$, has three eigenvalues, ${\lambda }_1=0.4796+0.7733i$, ${\lambda }_2=0.4796-0.7733i$ and ${\lambda }_3=0.7730$. Note that, $|{\lambda }_1|=|{\lambda }_2|=0.91<1$ and $|{\lambda }_3|<1$. Figure 8 shows that the interior equilibrium point, $E_{NPO}=(0.8542,1.0930,0.1754)$, is locally asymptotically stable. That is, small initial nutrients, phytoplankton and oysters values close to $E_{NPO}=(0.8542,1.0930,0.1754)$ levels lead to the persistence of nutrients, phytoplankton and oysters in Model NPO at $E_{NPO}$.

To study the impact of high nutrients flow on 8, we increase γ from $\gamma =3$ to $\gamma =4.5$ and keep all the other parameters fixed at their current values in Example 6.1. Then ${\mathcal{R}}_N$ increases from ${\mathcal{R}}_N=11.25$ to ${\mathcal{R}}_N=16.8750$, $|{\lambda }_1|=|{\lambda }_2|=1.0214>1$. Figure 9 illustrates the emergence of a Neimark–Sacker closed curve. As γ increase to $\gamma =4.5$, the interior equilibrium point undergoes a Neimark–Sacker bifurcation.

High flow of nutrients into a bay ecosystem is capable of forcing phytoplankton blooms with oscillations in oysters via Neimark–Sacker bifurcation.

Figure 8. Model NPO has a locally asymptotically stable interior equilibrium point, $E_{NPO}=(0.8542,1.0930,0.1754)$, where $\alpha =1.2$, $\gamma =3$, ${\mu }_N=0.2$, ${\mu }_P=0.6$, ${\mu }_O=0.09$ and $\delta =0.9$.

Example 6.2

In Model NPO, let $\alpha =1.2$, $\gamma =3.5$, ${\mu }_N=0.2$, ${\mu }_P=0.4$, ${\mu }_O=0.09$ and $\delta =0.9$.

Figure 9. Neimark–Sacker bifurcation closed curve in N-P, N-O and P-O planes at $\gamma =4.5$, where the other parameters are fixed at their current values in Example 6.1.

In Example 6.2, ${\mathcal{R}}_N=8.75>1$ and the Jacobian matrix of NPO model, ${\mathcal{J}}_{|E_{NP}}$ has three eigenvalues, ${\lambda }_1=0.4046+0.8746i$, ${\lambda }_2=0.4046-0.8746i$ and ${\lambda }_3=0.8805$. Note that, $|{\lambda }_1|=|{\lambda }_2|=0.9637<1$ and $|{\lambda }_3|<1$. Figure 10 shows that the interior equilibrium point, $E_{NPO}=(1.0369,1.0617,0.1094)$, is locally asymptotically stable. That is, small initial nutrients, phytoplankton and oysters close to $E_{NPO}=(1.0369,1.0617,0.1094)$ levels lead to the persistence of nutrients, phytoplankton and oysters in Model NPO at $E_{NPO}$.

To illustrate the instability of $E_{NPO}$ in Example 6.2. and Figure 10, we increase ${\mu }_p$ from ${\mu }_p=0.4$ to ${\mu }_p=0.9$ and keep all the other parameters fixed at their current values in Example 6.2. Then ${\mathcal{R}}_N$ increases from ${\mathcal{R}}_N=8.75$ to ${\mathcal{R}}_N=52.5$, $|{\lambda }_3|=2.3576>1$. Figure 11 illustrates the emergence of a Neimark–Sacker bifurcation.

Figure 10. Model NPO has the locally asymptotically stable and the interior equilibrium point, $E_{NPO}=(1.0369,1.0617,0.1094)$, where $\alpha =1.2$, $\gamma =3.5$, ${\mu }_N=0.2$, ${\mu }_P=0.4$, ${\mu }_O=0.09$ and $\delta =0.9$.

7. Discussion

In this paper, we introduced a discrete time mathematical model and used it to study the dynamics of interactions of nutrients, phytoplankton and oysters occurring in a bay ecosystem. We established that when the threshold level ${\mathcal{R}}_N<1$, then extinction of both populations of phytoplankton and oysters is attained and the discrete time model stabilizes with a positive value of nutrients. Also, we showed that when ${\mathcal{R}}_N>1$, the discrete time model is capable of exhibiting: (a) an equilibrium point where only nutrients and phytoplankton persist while oysters go extinct and (b) an equilibrium point where nutrients, phytoplankton and oysters are all persistent (interior equilibrium point). However, when ${\mathcal{R}}_N>1$, both equilibrium points with oysters only missing and interior equilibrium point are given implicitly as solutions of systems (6(6) $\{\begin{array}{l}\overline{N}(1-{\mu }_N{\mathrm{e}}^{-\alpha \overline{P}})=\gamma {\mathrm{e}}^{-\alpha \overline{P}},\\ \overline{N}(1-{\mathrm{e}}^{-\alpha \overline{P}})=(1-{\mu }_P)\overline{P}.\end{array}$(6) ) and (8(8) $\{\begin{array}{l}\gamma {\mathrm{e}}^{-\alpha P}=(1-{\mu }_N{\mathrm{e}}^{-\alpha P})N\\ N(1-{\mathrm{e}}^{-\alpha P}){\mathrm{e}}^{-\delta O}=(1-{\mu }_P{\mathrm{e}}^{-\delta O})P\\ P(1-{\mathrm{e}}^{-\delta O})=(1-{\mu }_O)O\end{array}$(8) ), respectively. Furthermore, we presented a set of conditions in Theorem 5.2 that guarantee the local stability of the oysters-free only equilibrium point.

Figure 11. A Neimark–Sacker bifurcation closed curve emerges in N-P, N-O and P-O planes at ${\mu }_p=0.9$, where the other parameters are fixed at their current values in Example 6.2.

We performed sensitivity analysis on the nutrients-only equilibrium point. Results of the local sensitivity analysis are the following:

• Increasing nutrients flow, γ, by $1\mathrm{\%}$ will increase nutrients level at equilibrium by $1\mathrm{\%}$.

• Increasing escapement rate of nutrients, ${\mu }_N$, will increase nutrients level at equilibrium by $2.3333\mathrm{\%}$.

Results of the global sensitivity analysis using Sobol' method at the phytoplankton- and oysters- free equilibrium point are the following:

• Escapement rate of nutrients, ${\mu }_N$, is the most important parameter and contributes about $100\mathrm{\%}$ of the total variance at the the nutrients equilibrium level.

Sensitivity analysis of both the oysters-free only equilibrium and the interior equilibrium was not performed due to the lack of explicit formulas of the oysters-free only equilibrium. We used examples to illustrate that human activities that lead to increase nutrients flow rate and increase phytoplankton survival in the bay ecosystem can force the interior equilibrium point to undergo a Neimark–Sacker bifurcation which generates phytoplankton blooms with corresponding oscillations in oysters population and nutrients level. Similar oscillatory bifurcations were also observed in the continuous time model for nutrients-phytoplankton-oysters interactions [25], where increasing either nutrients flow or phytoplankton removal rates forced the interior equilibrium point to undergo Hopf bifurcation.

The discrete-time and continuous-time NPO model have similarities that we summarize below:

• The nutrients level is asymptotically constant when there are no phytoplankton and $P_t=0$.

• The phytoplankton population goes extinct when $N_t=0$ and the nutrients level is very low.

• The oyster population goes extinct when $P_t=0$ and there are no phytoplankton.

• The threshold parameter ${\mathcal{R}}_N$ formulas are different in the two models. However, in both models when ${\mathcal{R}}_N<1$, then nutrients persist while phytoplankton and oysters population go extinct, and when ${\mathcal{R}}_N>1$, then both nutrients and phytoplankton persist with or without oysters.

• Results of local sensitivity indices for the nutrients-only equilibrium is consistent in both models.

• Dynamics behaviours of discrete-time and continuous-time NPO model can range from equilibrium to oscillatory dynamics via bifurcation.

However, results of Sobols' sensitivity indices of the nutrients-only equilibrium showed that the flow of nutrients contribute to the total variance in the continuous-time NPO model while it does not contribute in the discrete-time NPO model.

This discrete-time NPO model is a first step in using difference equations to model interactions of nutrients, phytoplankton and oysters in a bay ecosystem. For simplicity, the model does not take into account seasonality, temperature, light which play a role in the feeding patterns of phytoplankton and oyster populations. Including some of these factors in the discrete-time NPO model is an open question.

Acknowledgement

I am very thankful to Dr. Abdul-Aziz Yakubu for helpful discussions during the early stages of the project, particularly in the formulation of the model.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Najat Ziyadi was partially supported by a Simons Foundation Grant.