N(utrient)-P(hytoplankton)-Z(ooplankton)-D(etritus) Model

A toy interactive model of ocean ecosystem dynamics

Riley X. Brady

[email protected]


References:

  1. J.L. Sarmiento and N. Gruber (2006). Ocean Biogeochemical Dynamics. Chapter 4: "Production."

  2. A.M. Edwards (2001). Adding Detritus to a Nutrient–Phytoplankton–Zooplankton Model :A Dynamical-Systems Approach. J Plankton Res.


Summary

We can create a reduced model of a complex lower-trophic ocean ecosystem with a few differential equations. Here, we choose to model four variables: an arbitrary nutrient concentration (generally thought of as a macronutrient such as nitrate), a phytoplankton concentration (maybe a diatom), a zooplankton concentration (with option to change the type of zooplankter), and a detritus concentration (waste products).

In reality, we are just modeling a finite reservoir of nitrate, and considering how it gets redistributed around the ecosystem, given a few initial conditions and parameter settings. In other words, we aren't explicitly modeling phytoplankton cell count or biomass, but rather tracking where the nitrate goes as it is incorporated into organic matter via photosynthesis, or consumed by zooplankton.

Differential equations (DE's) are complex things to deal with. In a model like this, we have four DE's interacting with one another, because the rate of change of the given population (nutrient, phytoplankton, zooplankton, or detritus) is dependent on the current state of the other three populations. Thus, it is a lot easier to discretize a model into time steps, and reduce our DE's into algebraic equations that may be solved in reference to the current state of the system.

Here, I use an explicit time-differencing scheme (forward Euler method) to model this simple ocean ecosystem.


Differential Equations Contributing to the Model

The four raw DE's are as follows (where N is our nutrients, P is our phytoplankton, Z is our zooplankton, and D is our detritus):

$\frac{dN}{dt} = -V_{m}\left(\frac{N}{K_{N}+N}\right)f(I_{0})P~+~\alpha R_{m}\left(1-e^{-\lambda P}\right)Z~+~\epsilon P~+~gZ~+~\phi D$

$\frac{dP}{dt} = V_{m}\left(\frac{N}{K_{N} + N}\right)f(I_{0})P~-~R_{m}\left(1 - e^{-\lambda P}\right)Z~-~\epsilon P~-~rP$

$\frac{dZ}{dt} = \beta R_{m}\left(1-e^{-\lambda P}\right)Z~-~gZ$

$\frac{dD}{dt} = rP~+~(1 - \alpha - \beta)R_{m}\left(1 - e^{-\lambda P}\right)Z~-~\phi D$


Terms

Bulk Terms

If you look closely at each DE, you note that these are simply source (+) minus sink (-) equations. This simple model only has a few nitrogen exchange processes:

$V_{m}\left(\frac{N}{K_{N} + N}\right)f(I_{0})P$ : Phytoplankton grazing term. How much inorganic nitrogen are they taking up?

$R_{m}\left(1-e^{-\lambda P}\right)Z$ : Zooplankton grazing term. How much nitrogen are they taking up after consuming phytoplankton and releasing some as waste? The $\beta$ coefficient represents the proportion taken up into zooplankton organic matter; the $\alpha$ coefficient is that which is dissolved back into nutrients (perhaps from urine); and (1-$\alpha$-$\beta$) is that which is excreted as fecal pellets (to the detritus compartment).

$\epsilon P$ : How much nitrogen is being returned to the pool from phytoplankton death?

$gZ$ : How much nitrogen is being returned to the pool from zooplankton death?

$rP$ : How much nitrogen is being respired by phytoplankton into detritus?

Phytoplankton Terms

$V_{m}$ : Maximum growth rate of an individual plankter (div per day). This value is dependent on temperature, via a lab-derived equation for diatoms: $V_{m} = a\cdot b^{T}$, with $a=0.6d^{-1}$, $b=1.066$, and $c=1(degC)^{-1}$.

$K_{N}$ : Half-saturation constant for nitrogen uptake ($\mu molNl^{-1}$). This is the nitrogen concentration at which the phytoplankton growth rate is at half its maximum value.

$f_{0}$ : Light intensity (0 to 1 weighting function). This is a simple parameterization of a more complex hyperbolic term that uses a similar term to $K_{N}$.

$\epsilon$ : Phytoplankton death rate (cells per day).

$r$ : Respiration rate (per day).

Zooplankton Terms

$R_{m}$ : Maximum grazing rate of zooplankton on phytoplankton (cells per day).

$\lambda$ : Grazing constant ($\mu molNl^{-1}$).

$\beta$ : Proportion of assimilated nitrogen by zooplankton. In other words, when they graze upon a phytoplankter, how efficient are they at taking up the nitrogen? (dimensionless)

$\alpha$ : Proportion of nitrogen taken up by zooplankton that returns to the environment as dissolved nutrients (urine?).

$g$ : Zooplankton death rate (critters per day).

Detritus Terms

$\phi$ : The remineralization rate of detritus back into dissolved nutrients (per day).


Fixed Values/Initial Conditions

Parameter Symbol Default Value
Ambient Temperature T 15 degC
Half-saturation constant for $N$ uptake $K_{N}$ 1$\mu$mol per L
Maximum Grazing Rate R$_{m}$ 1 $d^{-1}$
Zooplankton Death Rate g 0.2$d^{-1}$
Zooplankton Grazing Constant $\lambda$ 0.2$\mu$mol per L
Phytoplankton Death Rate $\epsilon$ 0.1$d^{-1}$
Proportional Light Intensity f$_{0}$ 0.25
Zooplankton Dissolved Excretion Fraction $\alpha$ 0.3
Zooplankton Assimilation Efficiency $\beta$ 0.6
Phytoplankton Respiration Rate r 0.15
Detritus Remineralization Rate $\phi$ 0.4 $d^{-1}$

Tools

This model was built using Python 3 and visualized using Bokeh.

In [1]:
# Outside packages
import numpy as np

# Bokeh packages
from bokeh.io import output_notebook, show, gridplot, output_file, save
from bokeh.layouts import column, widgetbox
from bokeh.models import CustomJS, ColumnDataSource, Slider, FixedTicker
from bokeh.models.widgets import Slider, Dropdown, RadioButtonGroup
from bokeh.plotting import figure
from bokeh.charts import Area, Bar

# Set up colors (from ColorBrewer discrete colors)
cmap = ["#bebada", "#8dd3c7", "#fb8072", "#e5d8bd"]
In [2]:
# Allow Bokeh to be utilized inline with Jupyter.
output_notebook()
Loading BokehJS ...

Default Model View


We first need to compute the model in Python with some basic parameterizations that we know work. This will serve as the default view when the user opens the page.

In [3]:
# Here we set up the default parameters/coefficients. 
DT = 1 # Time Step (in days)
NUM_STEPS = 150 # Number of time steps to be computed and plotted

# Temperature-Dependent Growth Rate
a = 0.6
b = 1.066
c = 1
T = 15
Vm = a * b**(c*T) # Maximum growth rate (per day)

# Other parameters
Kn = 1    # Half-saturation constant for nitrogen uptake (umolN per l)
Rm = 1    # Maximum grazing rate (per day)
g = 0.2  # Zooplankton death rate (per day)
lambda_Z = 0.2  # Grazing constant (umolN per l)
epsilon = 0.1  # Phyto death rate (per day)
f = 0.25 # Light intensity (assumed constant)

# Detritus-related stuff.
alpha = 0.3 # Fraction of zoo. uptake that goes immediately to dissolved nutrients.
beta = 0.6  # Assimilation efficiency of zooplankton.
r = 0.15 # Respiration rate.
phi = 0.4 # Remineralization rate of detritus.
In [4]:
# Set Initial Conditions (umol per L)
N_0 = 4 
P_0 = 2.5 
Z_0 = 1.5
D_0 = 0
In [5]:
# Initialize Arrays
N = np.empty(NUM_STEPS, dtype="float")
P = np.empty(NUM_STEPS, dtype="float")
Z = np.empty(NUM_STEPS, dtype="float")
D = np.empty(NUM_STEPS, dtype="float")

# Insert Initial Values
N[0] = N_0
P[0] = P_0
Z[0] = Z_0
D[0] = D_0

Compute Simulation in Python

In [6]:
# Here we use the Euler forward method to solve for t+1 and reference t. 
for idx in np.arange(1, NUM_STEPS, 1):
    t = idx - 1
    
    # Common terms for simpler code
    gamma_N   = N[t] / (Kn + N[t])
    zoo_graze = Rm * (1 - np.exp(-lambda_Z * P[t])) * Z[t]
    
    # Equation calculations
    N[idx] = DT * (-Vm*gamma_N*f*P[t] + alpha*zoo_graze + epsilon*P[t] + g*Z[t] + phi*D[t]) + N[t] 
    P[idx] = DT * (Vm*gamma_N*f*P[t] - zoo_graze - epsilon*P[t] - r*P[t]) + P[t]
    Z[idx] = DT * (beta*zoo_graze - g*Z[t]) + Z[t]  
    D[idx] = DT * (r*P[t] + (1-alpha-beta)*zoo_graze - phi*D[t]) + D[t]

Set up Bokeh Data Structure

In [7]:
x = np.arange(1, NUM_STEPS + 1, 1)
N = N
P = P
Z = Z
D = D

# Bokeh likes reading data via its own version of dictionaries.
# I also prime this dictionary with additional variables for 
# multiple plots, which makes the custom JS interaction a lot easier to manage.
source = ColumnDataSource(data = {
        'x'    : x,
        'N'    : N,
        'P'    : P,
        'Z'    : Z,
        'D'    : D,
        'Psum' : N + P, # These sum variables can be removed... are used for a stacked bar plot.
        'Zsum' : N + P + Z,
        'Dsum' : N + P + Z + D,
    })
In [8]:
# Functions for plotting
def plotlines(plot, x, y, source, legend, line_width=3, line_alpha=0.75,
             color='black'):
    plot.line(x, y, source=source, line_width=line_width,
             line_alpha=line_alpha, color=color, legend=legend)

Standard Visualization (before interaction)

In [9]:
# TIME SERIES PLOT
plot = figure(plot_width=900, plot_height=300,
             toolbar_location="right", tools = "save",
             x_range=(1, NUM_STEPS + 1), title="N-P-Z-D Time Series",
             webgl=True)

# Plot data
plotlines(plot, 'x', 'N', source, "Nutrients", color=cmap[0])
plotlines(plot, 'x', 'P', source, "Phytoplankton", color=cmap[1])
plotlines(plot, 'x', 'Z', source, "Zooplankton", color=cmap[2])
plotlines(plot, 'x', 'D', source, "Detritus", color=cmap[3])


# Plot aesthetics
# Title
plot.title.align           = "center"
plot.title.text_font_style = "normal"
plot.title.text_font_size  = "14pt"

# Axes
plot.yaxis.axis_label                 = 'Concentration (umolN per L)'
plot.yaxis.axis_label_text_font_style = "normal"
plot.yaxis.axis_label_text_font_size  = "10pt"
plot.xaxis.axis_label                 = 'Model Days'
plot.xaxis.axis_label_text_font_style = "normal"
plot.xaxis.axis_label_text_font_size  = "10pt"

# Grid
plot.ygrid.grid_line_alpha       = 0.2
plot.xgrid.grid_line_alpha       = 0
plot.xgrid.minor_grid_line_color = 'grey'
plot.xgrid.minor_grid_line_alpha = 0.2
In [10]:
# STACKED BAR PLOT
plot2 = figure(plot_width=900, plot_height=300,
            toolbar_location="right", tools = "save",
            x_range=(1, NUM_STEPS + 1), y_range=plot.y_range, 
            title="Nutrient Distribution", webgl=True)

# Plot data
plot2.vbar('x', width=0.2, bottom=0, top='N', source=source, color=cmap[0])
plot2.vbar('x', width=0.2, bottom='N', top='Psum', source=source, color=cmap[1])
plot2.vbar('x', width=0.2, bottom='Psum', top='Zsum', source=source, color=cmap[2])
plot2.vbar('x', width=0.2, bottom='Zsum', top='Dsum', source=source, color=cmap[3])

# Aesthetics
plot2.y_range.start = 0

# Title
plot2.title.align           = "center"
plot2.title.text_font_style = "normal"
plot2.title.text_font_size   = "14pt"

# Axes
plot2.yaxis.axis_label                 = 'Concentration (umolN per L)'
plot2.yaxis.axis_label_text_font_style = "normal"
plot2.yaxis.axis_label_text_font_size  = "10pt"
plot2.xaxis.axis_label                 = 'Model Days'
plot2.xaxis.axis_label_text_font_style = "normal"
plot2.xaxis.axis_label_text_font_size  = "10pt"

Javascript Interaction

In [11]:
callback = CustomJS(args=dict(source=source), code="""
    // Ingest main model data for modification
    var data = source.get('data');
    x    = data['x'];
    N    = data['N'];
    P    = data['P'];
    Z    = data['Z'];
    D    = data['D'];
    // Again, these are optional and are to serve a stacked bar plot visualization.
    Psum = data['Psum'];
    Zsum = data['Zsum'];
    Dsum = data['Dsum'];
    
    // Parameters
    var Kn = 1; // Anything that is fixed and unable to be modified just gets the same value as the default view.
    var g = z_death.get('value'); // Anything that will be a slider or button gets this JS call.
    var lambda_Z = 0.2;
    var epsilon = p_death.get('value');
    var f = light.get('value');
    var dt = 1;
    var T = temperature.get('value');
    
    // Detritus-Related Parameters
    var alpha = 0.3;
    var beta = 0.6;
    var r = 0.15;
    var phi = 0.4;
    
    // Calculate Maximum Grazing Rate
    var a  = 0.6;
    var b  = 1.066;
    var Vm = a * Math.pow(b, T);
    
    // Zooplankton Species (impacts grazing rate)
    // Need to use IF statements for the button widget.
    var entry = zooSpecies.get('active');
    if (entry === 0) {
        var Rm = 1.6;
    } else if (entry === 1) {
        var Rm = 1.8;
    } else if (entry === 2) {
        var Rm = 1;
    } else {
        var Rm = 2;
    }


    // Initial Conditions with modifications allowed
    var N_0 = nut.get('value');
    var P_0 = phyto.get('value');
    var Z_0 = zoo.get('value');
    var D_0 = det.get('value');
    
    // Insert Initial Values for model
    N[0]    = N_0;
    P[0]    = P_0;
    Z[0]    = Z_0;
    D[0]    = D_0;
    Psum[0] = N_0 + P_0;
    Zsum[0] = N_0 + P_0 + Z_0;
    Dsum[0] = N_0 + P_0 + Z_0 + D_0;
    
    // Run Model
    for (i = 1; i < x.length; i++) {
         t = i - 1;

        // Common terms
        gamma_N   = N[t] / (Kn + N[t])
        zoo_graze = Rm * (1 - Math.exp(-lambda_Z * P[t])) * Z[t]
    
        // Equation calculations for model
        N[i] = dt * (-Vm*gamma_N*f*P[t] + alpha*zoo_graze + epsilon*P[t] + g*Z[t] + phi*D[t]) + N[t];
        P[i] = dt * (Vm*gamma_N*f*P[t] - zoo_graze - epsilon*P[t] - r*P[t]) + P[t];
        Z[i] = dt * (beta*zoo_graze - g*Z[t]) + Z[t];
        D[i] = dt * (r*P[t] + (1-alpha-beta)*zoo_graze - phi*D[t]) + D[t];
        
        // Sum Variables
        Psum[i] = N[i] + P[i];
        Zsum[i] = N[i] + P[i] + Z[i];
        Dsum[i] = N[i] + P[i] + Z[i] + D[i];
    }
    source.trigger('change');
""")

Building Sliders and Buttons

In [12]:
# Nutrient Initial Conditions
nut_slider = Slider(start = 0, end = 10, value = 4, step = 0.5, title = "Initial Nutrient Concentration", callback=callback)
callback.args["nut"] = nut_slider # The quoted "nut" is what is called in the JS interaction.

# Phytoplankton Initial Conditions
phyto_slider = Slider(start = 0, end = 10, value = 2.5, step = 0.5, title = "Initial Phytoplankton Concentration", callback=callback)
callback.args["phyto"] = phyto_slider

# Zooplankton Initial Conditions
zoo_slider = Slider(start = 0, end = 10, value = 1.5, step = 0.5, title = "Initial Zooplankton Concentration", callback=callback)
callback.args["zoo"] = zoo_slider

# Detritus Initial Conditions
det_slider = Slider(start = 0, end = 10, value = 0, step = 0.5, title = "Initial Detritus Concentration", callback=callback)
callback.args["det"] = det_slider

# Ambient Temperature
temp_slider = Slider(start = 0, end = 25, value = 15, step = 1, title = "Water Temperature (degC)", callback=callback)
callback.args["temperature"] = temp_slider

# Phytoplankton Death Rate
pdeath_slider = Slider(start = 0, end = 0.5, value = 0.1, step = 0.05, title = "Phytoplankton Natural Death Rate (per day)", callback=callback)
callback.args["p_death"] = pdeath_slider

# Zooplankton Death Rate
zdeath_slider = Slider(start = 0, end = 0.5, value = 0.2, step = 0.05, title = "Zooplankton Natural Death Rate (per day)", callback=callback)
callback.args["z_death"] = zdeath_slider

# Light Intensity
light_slider = Slider(start = 0, end = 1, value = 0.25, step = 0.05, title = "Proportional Light Intensity", callback=callback)
callback.args["light"] = light_slider

# Zooplankton Species (for grazing)
zoo_species = RadioButtonGroup(labels=["Cladoceran", "Copepod", "Mysid", "Rotifer"], active=2, 
                               callback=callback)
callback.args["zooSpecies"] = zoo_species

Final Layout and Saving the Plot

In [13]:
layout = gridplot([[nut_slider, phyto_slider, zoo_slider], [det_slider, pdeath_slider, zdeath_slider], 
                   [light_slider, temp_slider, zoo_species], 
                   [plot], [plot2]])

# Option 1: Display in the Notebook
show(layout)

# Option 2: Save to HTML to run in browser
# output_file('NPZD-Model.html')
# save(layout)