Plankton Model with Time Delayed Nutrient Recycling Sue Ann - - PowerPoint PPT Presentation

Plankton Model with Time Delayed Nutrient Recycling

Sue Ann Campbell, Matthew Kloosterman and Francis Poulin

Department of Applied Mathematics University of Waterloo

Southern Ontario Dynamics Day April 12, 2013

Outline

1

Introduction/Background

2

Existence of Equilibrium Points

3

Stability of Equilibrium Points No Delay With Delay

4

Conclusions and Implications

Introduction

Plankton are free floating organisms found in oceans and lakes which form the bottom of the food chain.

Introduction

Phytoplankton are plankton which carry out photosynthesis examples: diatoms, golden algae, green algae and cyanobacteria

Introduction

Zooplankton are plankton that feed on phytoplankton examples: jelly fish, small crustaceans and insect larvae

Motivation

Why study plankton? Plankton form the bottom of the ocean food chain.

Motivation

Why study plankton? Plankton form the bottom of the ocean food chain. Phytoplankton can exhibit blooms which can be harmful to ecosystem and humans.

Motivation

Why study plankton? Plankton form the bottom of the ocean food chain. Phytoplankton can exhibit blooms which can be harmful to ecosystem and humans. Phytoplankton are very important in the transfer of carbon dioxide from the atmosphere to the ocean.

Model

Closed model with three compartments: dissolved nutrient - N(t) phytoplankton - P(t) zooplankton - Z(t) (measured by amount of limiting nutrient/nitrogen)

P .J.S. Franks (2002) J Oceanogr. 58:379-387.

Model

Closed model with three compartments: dissolved nutrient - N(t) phytoplankton - P(t) zooplankton - Z(t) (measured by amount of limiting nutrient/nitrogen) Nutrient Zooplankton Phytoplankton Recycling Recycling Uptake Grazing

P .J.S. Franks (2002) J Oceanogr. 58:379-387.

Model

Model with three compartments: dissolved nutrient - N(t) phytoplankton - P(t) zooplankton - Z(t) (measured by amount of limiting nutrient/nitrogen) N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t)

Model

Model with three compartments: dissolved nutrient - N(t) phytoplankton - P(t) zooplankton - Z(t) (measured by amount of limiting nutrient/nitrogen) N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) phytoplankton nutrient uptake

Model

Model with three compartments: dissolved nutrient - N(t) phytoplankton - P(t) zooplankton - Z(t) (measured by amount of limiting nutrient/nitrogen) N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) zooplankton grazing on phytoplankton

Model

Model with three compartments: dissolved nutrient - N(t) phytoplankton - P(t) zooplankton - Z(t) (measured by amount of limiting nutrient/nitrogen) N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) zooplankton grazing on phytoplankton nutrient recycling

Model

Model with three compartments: dissolved nutrient - N(t) phytoplankton - P(t) zooplankton - Z(t) (measured by amount of limiting nutrient/nitrogen) N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) zooplankton and phytoplankton death

Model

Model with three compartments: dissolved nutrient - N(t) phytoplankton - P(t) zooplankton - Z(t) (measured by amount of limiting nutrient/nitrogen) N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) zooplankton and phytoplankton death nutrient recycling

Model Parameters

Parameter Meaning Units µ phytoplankton maximum growth rate day−1 λ phytoplankton death rate day−1 g zooplankton maximum grazing rate day−1 γ zooplankton assimilation efficiency δ zooplankton death rate day−1 N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t)

Functional Response

Nutrient uptake by phytoplankton: µP(t)f(N(t)) f(0) = 0, f ′(N) ≥ 0, f ′′(N) ≤ 0, lim

N→∞ f(N) = 1 (Michaelis-Menten/Type II)

W.C. Gentleman & A.B. Neuheimer (2008)

• J. Plankton Research 30(11) 1215-1231.

Functional Response

Nutrient uptake by phytoplankton: µP(t)f(N(t)) f(0) = 0, f ′(N) ≥ 0, f ′′(N) ≤ 0, lim

N→∞ f(N) = 1 (Michaelis-Menten/Type II)

Zooplankton grazing on phytoplankton: gZ(t)h(P(t)) h(0) = 0, h′(P) ≥ 0, lim

P→∞ h(P) = 1 (Type II or III)

W.C. Gentleman & A.B. Neuheimer (2008)

• J. Plankton Research 30(11) 1215-1231.

Functional Response

Nutrient uptake by phytoplankton: µP(t)f(N(t)) f(0) = 0, f ′(N) ≥ 0, f ′′(N) ≤ 0, lim

N→∞ f(N) = 1 (Michaelis-Menten/Type II)

Zooplankton grazing on phytoplankton: gZ(t)h(P(t)) h(0) = 0, h′(P) ≥ 0, lim

P→∞ h(P) = 1 (Type II or III)

W.C. Gentleman & A.B. Neuheimer (2008)

• J. Plankton Research 30(11) 1215-1231.

Model

N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) Nutrient Zooplankton Phytoplankton Recycling Recycling Uptake Grazing

Model

N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) Include distributed time delay in recycling

Model

N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) Include distributed time delay in recycling N′(t) = ∞ [λP(t − u) + δZ(t − u) + (1−γ)gZ(t − u)h(P(t − u))]η(u) du −µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) where ∞ η(u) du = 1, τ = ∞ u η(u) du (mean delay)

Model

N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) Include distributed time delay in recycling N′(t) = ∞ [λP(t − u) + δZ(t − u) + (1−γ)gZ(t − u)h(P(t − u))]η(u) du −µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) where ∞ η(u) du = 1, τ = ∞ u η(u) du (mean delay) Recycling time is u ∈ [0, ∞) with probability η(u).

Distributions

Gamma distribution: η(u) = up−1 p

τ

p e−pu/τ Γ(p) Uniform distribution: η(u) =

• 1

2W ,

τ − W ≤ u ≤ τ + W 0, elsewhere , Tent distribution: η(u) =   

u+W−τ W 2

, τ − W ≤ u ≤ τ

−u+W+τ W 2

, τ ≤ u ≤ τ + W 0, elsewhere . Discrete delay: η(u) = δ(u − τ)

Distributions (τ = 2)

g 0.3 0.1 u 0.6 0.5 0.4 0.2 0.0 5 4 3 2 1

Gamma (p = 1, 2, 4, 8) Uniform (W = 0.5, 1, 2) Tent (W = 0.5, 1, 2)

Conservation Laws

Model with no delay: N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t)

Conservation Laws

Model with no delay: N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) Total nutrient in system is conserved. N(t) + P(t) + Z(t) = NT (constant)

Conservation Laws

Model with delay: N′(t) = ∞ [λP(t − u) + δZ(t − u) + (1−γ)gZ(t − u)h(P(t − u))]η(u) du −µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t)

Conservation Laws

Model with delay: N′(t) = ∞ [λP(t − u) + δZ(t − u) + (1−γ)gZ(t − u)h(P(t − u))]η(u) du −µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) Total nutrient in system is conserved. NT = N(t) + P(t) + Z(t) + ∞ t

t−u

[λP(v) + δZ(v) + (1 − γ)gZ(v)h(P(v))]η(u) dv du

• nutrient being recycled

Equilibrium Points

Model with delay: N′(t) = ∞ [λP(t − u) + δZ(t − u) + (1−γ)gZ(t − u)h(P(t − u))]η(u) du −µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) Equilibrium points: (N(t), P(t), Z(t)) = (N∗, P∗, Z ∗) constant.

Equilibrium Points

Model with delay: N′(t) = ∞ [λP(t − u) + δZ(t − u) + (1−γ)gZ(t − u)h(P(t − u))]η(u) du −µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) Equilibrium points: (N(t), P(t), Z(t)) = (N∗, P∗, Z ∗) constant. Must satisfy µP∗f(N∗) − gZ ∗h(P∗) − λP∗ = γgZ ∗h(P∗) − δZ ∗ =

Equilibrium Points

Model with delay: N′(t) = ∞ [λP(t − u) + δZ(t − u) + (1−γ)gZ(t − u)h(P(t − u))]η(u) du −µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) Equilibrium points: (N(t), P(t), Z(t)) = (N∗, P∗, Z ∗) constant. Must satisfy µP∗f(N∗) − gZ ∗h(P∗) − λP∗ = γgZ ∗h(P∗) − δZ ∗ = and conservation law: NT = N∗ + P∗ + Z ∗ + [λP∗ + δZ ∗ + (1 − γ)gZ ∗h(P∗)]τ

Equilibrium Points - Existence and Uniqueness

For each value of NT there exists a unique equilibrium point

• f each of the following types:

Trivial: (NT, 0, 0) - lies in positive orthant if NT > 0

Equilibrium Points - Existence and Uniqueness

For each value of NT there exists a unique equilibrium point

• f each of the following types:

Trivial: (NT, 0, 0) - lies in positive orthant if NT > 0 Phytoplankton: ( ˆ N, ˆ P, 0) where ˆ N = f −1(λ/µ), ˆ P = NT − f −1(λ/µ) 1 + λτ lies in positive orthant if NT > NT1 = f −1

λ µ

Equilibrium Points - Existence and Uniqueness

For each value of NT there exists a unique equilibrium point

• f each of the following types:

Trivial: (NT, 0, 0) - lies in positive orthant if NT > 0 Phytoplankton: ( ˆ N, ˆ P, 0) where ˆ N = f −1(λ/µ), ˆ P = NT − f −1(λ/µ) 1 + λτ lies in positive orthant if NT > NT1 = f −1

λ µ

• Coexistence: (N∗, P∗, Z ∗) where

P∗ = h−1 δ γg

• Z ∗ = γP∗

δ [µf(N∗) − λ] NT = N∗ + h−1 δ γg 1 − γλ δ + γ δ + τ

• µf(N∗)
• lies in positive orthant if NT > NT2 = f −1

λ µ

• + (1 + λτ)h−1

δ γg

Equilibrium Points - Stability with No Delay

N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) N(t) + P(t) + Z(t) = NT Using linearization and invariance of axes can show, for fixed NT

Equilibrium Points - Stability with No Delay

N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) N(t) + P(t) + Z(t) = NT Using linearization and invariance of axes can show, for fixed NT If 0 < NT < NT1 then (NT, 0, 0) is globally (asymptotically) stable. If NT1 < NT < NT2 then ( ˆ N, ˆ P, 0) is globally (asymptotically) stable, (NT, 0, 0) is unstable. If NT2 < NT, then (NT, 0, 0) and ( ˆ N, ˆ P, 0) are unstable. Stability of (N∗, P∗, Z ∗) depends on form of h(P).

Equilibrium Points - Stability with No Delay

N′(t) = λP(t) + δZ(t) + (1 − γ)gZ(t)h(P(t)) − µP(t)f(N(t)) P′(t) = µP(t)f(N(t)) − gZ(t)h(P(t)) − λP(t) Z ′(t) = γgZ(t)h(P(t)) − δZ(t) N(t) + P(t) + Z(t) = NT Using linearization and invariance of axes can show, for fixed NT h(P) type II: Exists NT3 such that.

If NT2 < NT < NT3 then (N∗, P∗, Z ∗) is asymptotically stable. When NT3 = NT then characteristic equation has a pair of pure imaginary roots. If NT3 < NT then (N∗, P∗, Z ∗) is unstable.

h(P) type III: Stability depends P∗ = h−1( δ

γg )

If h(P∗) ≤ P∗h′(P∗) and NT2 < NT then (N∗, P∗, Z ∗) asymptotically stable. If h(P∗) > P∗h′(P∗) then stability of (N∗, P∗, Z ∗) is as for Type II.

Model Parameters

Parameter Meaning Value µ phytoplankton maximum growth rate 5.9 day−1 λ phytoplankton death rate 0.017 day−1 KN half saturation constant for N uptake 1.0 µM N g zooplankton maximum grazing rate 7 day−1 γ zooplankton assimilation efficiency 0.7 δ zooplankton death rate 0.17 day−1 KP half saturation constant for Z grazing on P 1.0 µM N Functional response for phytoplankton nutrient uptake: f(N) =

N N+KN

Functional response for zooplankton grazing: h(P) =

P P+KP (Type II)

• r

h(P) =

P2 P2+K 2

P (Type III)

References: F.J. Poulin & P .J.S. Franks J. Plankton Research 32(8) (2010) 1121-1130. A.E. Edwards J. Plankton Research 23(4) (2001) 389-413.

Model without Delay

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 −0.01 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

NT1 NT2

NT P

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 −5 5 10 15 20 x 10

−3

NT2

NT

Z

Transcritical bifurcations at NT = NT1 and NT = NT2

Model without Delay (Type II Functional Response)

0.2 0.4 0.6 0.8 1 1.2 1.4 0.2 0.4 0.6 0.8 1 1.2 1.4 N T N P Z

Hopf bifurcation at NT = NT3

Model without Delay (Type II Functional Response)

Numerical Simulations

1 2 3 4 5 6 7 8 9 x 10

−4

1 2 3 4 5 6 7 8 9 x 10

−4

P Z

NT < NT1

0.005 0.01 0.015 0.02 0.025 0.005 0.01 0.015 0.02 0.025

P Z

NT1 < NT < NT2

0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5

P Z

NT2 < NT < NT3

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 1.2

P Z

NT > NT3

Model with Delay - Stability

Theorem: Equilibrium (NT, 0, 0) is stable/unstable if NT

< > NT1 = f −1 λ µ

• Equilibrium ( ˆ

N, ˆ P, 0) is stable for any NT and distribution satisfying f −1

λ µ

• + 2λ(1+τλ)

µa

< NT < f −1

λ µ

• + (1 + λτ)h−1

δ γg

• = NT2(τ)

Equilibrium ( ˆ N, ˆ P, 0) is unstable for any NT and distribution satisfying NT > NT2(τ)

Model with Delay - Stability

Theorem: Equilibrium (NT, 0, 0) is stable/unstable if NT

< > NT1 = f −1 λ µ

• Equilibrium ( ˆ

N, ˆ P, 0) is stable for any NT and distribution satisfying f −1

λ µ

• + 2λ(1+τλ)

µa

< NT < f −1

λ µ

• + (1 + λτ)h−1

δ γg

• = NT2(τ)

Equilibrium ( ˆ N, ˆ P, 0) is unstable for any NT and distribution satisfying NT > NT2(τ)

• Proof. Characteristic equation of linearization about (NT, 0, 0):

s (s + δ) (s − µf(NT) + λ) = 0. Characteristic equation of linearization about ( ˆ N, ˆ P, 0): (s − γgd + δ)[s2 + µˆ Pas + µˆ Paλ(1 − ˆ η(s))] = 0. where a = f ′( ˆ N), d = h(ˆ P). Apply Rouché’s Theorem.

Model with Delay - Stability

Theorem: Equilibrium (NT, 0, 0) is stable/unstable if NT

< > NT1 = f −1 λ µ

• Equilibrium ( ˆ

N, ˆ P, 0) is stable for any NT and distribution satisfying f −1

λ µ

• + 2λ(1+τλ)

µa

< NT < f −1

λ µ

• + (1 + λτ)h−1

δ γg

• = NT2(τ)

Equilibrium ( ˆ N, ˆ P, 0) is unstable for any NT and distribution satisfying NT > NT2(τ)

100 200 300 400 500 0.01 0.02 0.03 0.04 0.05

τ NT

1a 1b 2 3 4

Model with Delay - Stability of ( ˆ N, ˆ P, 0)

Characteristic equation for ( ˆ N, ˆ P, 0): (s − γgd + δ)[s2 + µˆ Pas + µˆ Paλ(1 − ˆ η(s))] = 0. Boundary of stability region corresponds to points in parameter space where characteristic equation has a pair of pure imaginary roots.

Model with Delay - Stability of ( ˆ N, ˆ P, 0)

Characteristic equation for ( ˆ N, ˆ P, 0): (s − γgd + δ)[s2 + µˆ Pas + µˆ Paλ(1 − ˆ η(s))] = 0. Boundary of stability region corresponds to points in parameter space where characteristic equation has a pair of pure imaginary roots. Set s = ±iω in second factor: −ω2 + µˆ Paiω + µˆ Paλ(1 − ˆ η(iω))] = 0.

Model with Delay - Stability of ( ˆ N, ˆ P, 0)

Characteristic equation for ( ˆ N, ˆ P, 0): (s − γgd + δ)[s2 + µˆ Pas + µˆ Paλ(1 − ˆ η(s))] = 0. Boundary of stability region corresponds to points in parameter space where characteristic equation has a pair of pure imaginary roots. Set s = ±iω in second factor: −ω2 + µˆ Paiω + µˆ Paλ(1 − ˆ η(iω))] = 0. Parameterizing distribution, η(u), in terms of mean delay, τ, define C(ω, τ) = Re[ˆ η(iω)], S(ω, τ) = −Im[ˆ η(iω)] then boundary is determined by −ω2 + µˆ Paλ[1 − C(ω, τ)] = ω + λS(ω, τ) = ⇒ τ = τc(ˆ P) ω = ωc(ˆ P) and NT = f −1(λ/µ) + [1 + λτc(ˆ P)]ˆ P

Model with Delay - Stability of ( ˆ N, ˆ P, 0)

Gamma distribution with p = 1, 2: no solution for τc, ωc ( ˆ N, ˆ P, 0) stable for any NT and τ satisfying NT < NT2(τ)

Model with Delay - Stability of ( ˆ N, ˆ P, 0)

Gamma distribution with p = 1, 2: no solution for τc, ωc ( ˆ N, ˆ P, 0) stable for any NT and τ satisfying NT < NT2(τ) Discrete delay: boundary of stability region is given by τc(ˆ P) =   

1 ωc

• π − sin−1

− ωc

λ

• if 0 < µˆ

Pa ≤ λ

1 ωc

• 2π + sin−1

− ωc

λ

• if λ < µˆ

Pa < 2λ. NTc(ˆ P) = f −1

λ µ

• + [1 + λτc(ˆ

P)]ˆ P where ωc =

• 2µˆ

Paλ − (µˆ Pa)2.

Model with Delay - Stability of ( ˆ N, ˆ P, 0)

Exact region of stability

100 200 300 400 500 0.01 0.02 0.03 0.04 0.05

τ NT

1a 1b 2 3 4

Model with Delay - Stability of ( ˆ N, ˆ P, 0)

Exact region of stability

100 200 300 400 500 0.01 0.02 0.03 0.04 0.05

τ NT

1a 1b 2 3 4

Parameter values: as before (Poulin & Franks (2010) Other parameters: τ ∼ 5 − 250 days; NT ∼ 1 − 15 mmol N m−3 (A.E. Edwards J. Plankton Research 23(4) (2001) 389-413).

Model with Delay - Stability of (N∗, P∗, Z ∗)

Recall: P∗ = h−1

δ γg

• , Z ∗ = γP∗

δ [µf(N∗) − λ]

NT = N∗ + h−1 δ γg 1 − γλ δ + γ δ + τ

• µf(N∗)
• Characteristic equation:

s3 + a2(N∗)s2 + a1(N∗)s + a0(N∗) + [b1(N∗)s + b0(N∗)]ˆ η(s) = 0

Model with Delay - Stability of (N∗, P∗, Z ∗)

Recall: P∗ = h−1

δ γg

• , Z ∗ = γP∗

δ [µf(N∗) − λ]

NT = N∗ + h−1 δ γg 1 − γλ δ + γ δ + τ

• µf(N∗)
• Characteristic equation:

s3 + a2(N∗)s2 + a1(N∗)s + a0(N∗) + [b1(N∗)s + b0(N∗)]ˆ η(s) = 0 Characteristic equation with s = ±iω is equivalent to B(ω, N∗) C(ω, τ) S(ω, τ)

• = y(ω, N∗)

⇒ ω = ωc(N∗), τ = τc(N∗) Determines boundary of region of stability in τ, NT parameter space

Model with Discrete Delay - Stability of (N∗, P∗, Z ∗)

C(ω, τ) = cos(ωτ), S(ω, τ) = sin(ωτ)

5 10 15 20 0.5 1 1.5 1 3 5 1 3 5

τ NT

100 200 300 400 500 0.5 1 1.5 2 2.5 3 3.5 4 1 3 5

τ NT

Type II - h(P) = P P + KP Type III - h(P) = P2 P2 + K 2

P

Physical values: τ ∼ 5 − 250 days; NT ∼ 1 − 15 mmol N m−3

Model with Distributed Delay - Stability of (N∗, P∗, Z ∗)

5 10 15 20 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

p = 1 Gamma Distribution

τ NT

p = 2 p = 5 p = 10 p = 20

1 3 5 5 10 15 20 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 1 3 5

τ NT

Uniform Distribution

5 10 15 20 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 1 3

Tent Distribution

τ NT

5

W = 0.001, 1, 3, 4, 5, 6 W = .001, 1, 3, 5, 6, 7

Model with Distributed Delay - Stability of (N∗, P∗, Z ∗)

5 10 15 20 0.5 1 1.5

τ NT

Variance: 1 day2

1 3 5 Gamma Uniform Tent Discrete Existence Boundary 5 10 15 20 0.5 1 1.5 2 1 3 5

τ NT

Variance: 8 day2

Gamma Uniform Tent Discrete Existence Boundary

Model with Gamma Distributed Delay - Simulations

Simulations p = 20, NT = 0.5

50 100 150 200 250 300 0.1762 0.1764 0.1766 0.1768

N t (Days)

50 100 150 200 250 300 0.0359 0.0359 0.036

P t (Days)

50 100 150 200 250 300 0.1284 0.1284 0.1285 0.1286 0.1286

Z t (Days)

50 100 150 200 250 300 0.142 0.143 0.144 0.145

N t (Days)

50 100 150 200 250 300 0.035 0.0355 0.036 0.0365

P t (Days)

50 100 150 200 250 300 0.1065 0.107 0.1075 0.108

Z t (Days)

50 100 150 200 250 300 0.1144 0.1146 0.1148 0.115

N t (Days)

50 100 150 200 250 300 0.0359 0.0359 0.0359 0.036 0.036

P t (Days)

50 100 150 200 250 300 0.0872 0.0873 0.0873 0.0874 0.0875 0.0875

Z t (Days)

τ = 5 τ = 8 τ = 12

5 10 15 20 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

p = 1 Gamma Distribution

τ NT

p = 2 p = 5 p = 10 p = 20

1 3 5

Conclusions

Characterized total nutrient needed to sustain phytoplankton in isolation and in coexistence with zooplankton (NT1 < NT2). NT2 increases as time delay in recycling increases, as less biomass is available to sustain organisms. Type II functional response for phytoplankton grazing is “less stable” than type III in the following sense:

With type II coexistence equilibrium can be destabilized for sufficiently large total nutrient (NT > NT3), leading to oscillations If NT < NT3 type III needs larger delay to destabilize coexistence equilibrium.

Small delay can be stabilizing - coexistence equilibrium is stable for larger values of NT. If variance in distribution of delays is small, then actual distribution not important in determining stability. Kloosterman, Campbell & Poulin J. Mathematical Biology (2013).

Acknowledgements

Matt Kloosterman, Francis Poulin Natural Sciences and Engineering Research Council of Canada Ontario Graduate Scholarships Program