Exemplary Model Description Draft 1

Author: Bradley Potter

On this page I will describe the mathematical construction of the model for oyster reef development and bistability as developed in “Bistability in a Differential Equation Model of Oyster Reef Height and Sediment Accumulation” by William C.Jordan-Cooley, Romuald N. Lipcius, Leah B. Shaw, Jian Shen, and Junping Shi. This model consists of three state variables, three differential equations, and three functions which fall within the differential equations. How these equations were constructed, what assumptions were made, and what the role of the parameters are will all be addressed below.

State Variables and Parameters

State Variables

The three state variables are given in the table below:

State Variable	Units	Description
O(t)	Cubic meters	Live oyster volume
B(t)	Cubic meters	Dead oyster volume
S(t)	Cubic meters	Sediment volume on reef

Parameters

There are also eleven parameters which will be used in the formulation of the model. The parameters, along with their units, a brief description, and a typical range of values, are all given below:

Parameter	Units	Description	Typical Range
r	1 / year	Rate of growth	0.7 - 1.3
K	Cubic meters	Oyster carrying capacity	0.1 - 0.3
µ	1 / year	Death rate due to predation and disease	0.2 - 0.6
ε	1 /year	Death rate due to suffocation from sediment	0.94
γ	1 / year	Rate of shell degradation	0.5 - 0.9
F0	1 / year	Max rate of sediment filtration	1
C	Cubic meters / year	Max rate of sediment deposition	0.04 - 0.08
y0	Year / cubic meters	Amount of sediment at location of max filtration	0.02
β	1 / cubic meters	Rate of erosion of sediment	0.02 - 0.04
h	1 / cubic meters	Scaling factor	10 - 30
η	1 / cubic meters	Rate of decay of sediment deposition	3.33

Assumptions

The model described above is based on several key assumptions. First it is assumed that oysters above the sediment line will grow at a logistic rate while dying at a linear rate. Oysters which are below the sediment line are essentially suffocated or choked out by the sediment, which causes the oyster population to die at rate proportional to the total volume of live oysters below the sediment line. The function, f(d), ranges from 0 to 1 and represents the fraction of the oyster volume that is above the sediment line (1 means all of the oysters are above the sediment, while 0 means that the entire reef is buried). Once shells die, they degrade at a rate proportional to the total volume of dead oyster shells. Sediment gets deposited onto the reef at a decreasing rate with increasing reef height because water flow rates at lower depths are slower and will allow for more rapid settling of the sediment. Sediment deposition rate also decreases exponentially relative to the oysters' sediment filtration rate. Sediment on the reef is eroded away at a rate proportional to the volume of sediment, and it gets filtered by the oysters at a rate proportional to the total volume of live oysters. This filtration rate is assumed to have a maximum value for an optimal sediment concentration. It is assumed that there is an exponential decay of sediment concentration with increasing height within the body of water. Lastly, the filtration rate is assumed to be a Ricker-type function of the concentration of sediment in the water, meaning that it has a specific rate of change and a limiting, maximum filtration rate which we defined as F0.

Model Inputs and Descriptions

Below I describe each of the six equations used in the model, three differential equations and three functions used within the differential equations. These six equations act as the mathematical model of the oyster reef system.

Equation 1:

$\frac{dO}{dt} = r O(t) f(d) \left ( 1 - \frac{O(t)}{k} \right ) - \mu f(d) O(t) - \epsilon O(t) \left ( 1 - f(d) \right )$

This equation gives the rate of change of the live oyster population with respect to time. The first set of terms gives the logistic growth of the oyster population which is limited by the carrying capacity parameter k and has a growth rate parameter of r. Again, the function f(d) is an expression for the fraction of the live oysters which are above the sediment level. The total oyster volume, by definition and mathematics, cannot exceed the carrying capacity k. The next term in this equation represents the death of oysters which are not covered by sediment. These oysters die as a result of predation and disease, with a proportionality constant equal to the parameter mu. It is assumed that only oysters above the sediment die of these causes, while oysters covered by sediment only die due to suffocation. This is represented in the third term in the equation, which has a rate parameter epsilon for the mortality rate of oysters that are buried under and suffocated by the sediment.

Equation 2:

$\frac{dB}{dt} = \mu f(d) O(t) + \epsilon O(t) \left ( 1 - f(d) \right ) - \gamma B(t)$

This equation is a representation of the rate of change of the dead oyster population with respect to time. The increase in dead oyster volume is a direct result of the death of living oysters. Therefore, the first two terms in this equation are just the two death terms from equation 1, however here these terms are positive because a decrease in living oysters translates into an increase in dead oysters. The last term accounts for the degradation of the dead oyster shells. This degradation rate is proportional to the volume of dead oyster shells with a proportionality constant equal to the parameter gamma.

Equation 3:

$\frac{dS}{dt} = - \beta S(t) + C g \exp \left ( \frac{-F O(t)}{C g} \right )$

This equation represents the change in sediment volume on the oyster reef with respect to time. The first term represents the natural erosion of the sediment caused by ocean currents, which decreases the amount of sediment left on the reef at a rate proportional to the total volume of sediment with a proportionality constant equal to beta. The second term represents the deposition of sediment onto the reef. Sediment is deposited more rapidly when water flow rates are slower, which is the case at lower reef heights (deepest depth, meaning the bay floor). The deposition rate in the absence of a reef is equal to C, the maximum deposition rate. The argument to g, defined below, is the height of the reef, or O(t) + B(t), which in the absence of a reef is zero. Therefore g(0) = 1, and g(infinity) = 0. Biodeposition, meaning the waste products of other aquatic life that fall to the floor, are simply incorporated into the other parameters and is not explicitly included in the system of equations. The filtration of sediment (F) by the live oyster population is accounted for by the exponential term, where the filtration is proportional to the volume of live oysters and inversely proportional to the sediment concentration which is represented by Cg(x) at height x = O + B in the reef. The filtration rate is limited by parameter y0 due to clogging of the oysters' gills at excessive sediment concentrations. This is shown by F(y) below.

Function 1:

$f(d) = \frac{1}{1 + \exp (-h d)} \qquad \longrightarrow \qquad d = \frac{O(t)}{2} + B(t) - S(t)$

This function shows the fraction of live oysters which are above the sediment layer. The argument is a function of the three state variables, essentially a difference between the height (or volume) of the oysters, both living and dead, and the height (or volume) of the sediment. There is a scaling parameter, h, which defines the exact sigmoidal shape of the function, so that as h increases the transition from f(d) = nearly 0 to f(d) = nearly 1 happens more rapidly. The key transition point takes place when d is near zero. Below zero the function equals essentially 0, and above 0 the function equals essentially 1.

Function 2:

$g(x) = \exp (- \eta x) \qquad \longrightarrow \qquad x = O(t) + B(t)$

This function takes as its argument the height of the oyster reef, O + B, and is an important part of equation 2 above as it scales the sediment deposition rate. The sediment deposition rate varies with height, due largely to the change in water flow rate (which is inversely related to the deposition rate). The function g(x) follows exponential decay with increasing reef height, meaning that as the reef is taller the deposition rate of sediment on the reef decreases exponentially with an exponential scaling parameter eta.

Function 3:

$F(y) = \frac{F_0}{y_0} y \exp \left ( \frac{y_0 - y}{y_0} \right ) \qquad \longrightarrow \qquad y = C g(x)$

Lastly, this function give the filtration rate of the oysters. There is a maximum in the filtration rate, F0, which takes place at the concentration of sediment y0. As y exceeds y0, the gills of the oysters become clogged and cannot filter any faster. However, when y is less than y0, the oysters' gills are not at full capacity and could filter more sediment if it were present. This is the so-called Ricker model or Ricker function relating filtration rate to sediment concentration y = Cg.

Simulation Method

The system of equations described above was solved and simulated using Mathematica. Initially the three functions, f(d), g(x), and F(y), were plotted as d was varied between -1 and 1, x was varied between 0 and 1, and y was varied between 0 and 0.1. This step was included so that we could reproduce the graphs in figure 2 of the original paper, which confirmed that the functions were set up correctly.

Once the functions were up and running, we went about trying to perform a numerical simulation of the system. This is a key factor in analyzing the stability of the system for various parameter values and initial conditions. Using a nest of the Manipulate, NDSolve, Plot, Evaluate, and Show function in Mathematica, we created a numerical simulation which plotted the response of the three state variables, O, B, and S, over time from 0 to 50 years. This method was successful in reproducing such plots, and is a useful tool because the parameters and the initial value of B could be varied to see how the time-response would change. In an effort to reproduce the specific plots shown in figure 4 of the original paper, we used the following set of parameters (given in table 3 of the original paper) with initial B values of 0.2, 0.1, 0.12, and 0.11, respectively.

Parameter	r	K	µ	ε	γ	η	y0	F0	β	h	C
Value	1	0.3	0.4	0.94	0.7	3.33	0.02	1	0.01	20	0.02

Unfortunately, upon plotting the time-response using our manipulate function and the values given above, our system exhibited the bistability shown in the plots of figure 4, however the exact values at which the long-term stability changes does not match the response shown in the original paper. We are still trying to determine what the cause of this may be, though my current belief is that the parameter values given above are not what was ultimately used by the original research group in creating the plots shown in figure 4. I will update these results when possible.

The final portion of the paper we must reproduce before expanding upon the original paper is to create bifurcation diagrams of each of the state variables. The paper gives bifurcation diagrams of the state variables while varying parameter C, getting a saddle-node bifurcation near C = 0.078. We have not completed this portion of the model yet, due largely to the delays encountered in figuring out why our time-response plots don't match those provided. We want to make sure our model is set up correctly and functioning as expected before progressing too far into the model.