Exemplary Final Term Paper 1

Author: Bradley Potter

Introduction

This paper looks at the mathematical model of oyster reef populations in terms of reef volume. The model was created in Mathematica 8.0 for Students through the cooperative effort between my partner and I. The model was based on an original work on oyster reef bistability by William Jordan-Cooley, et al. After full development of the original model, an extension was added based on other biological references and research to acquire a deeper and more complete understanding of the system. The complete Mathematica file is attached here for your reference.

Significance of the Problem

The native oyster population in the Chesapeake Bay has historically been a robust species, tolerating a wide range of salt concentrations and water temperatures, but it has still fallen victim to pollution, disease, and over-harvesting. As a result, the population size has dropped to as low as 1% of its prior levels. Human activity has played a significant role in the degradation of these reefs: a combination of pollution; over-harvesting; destruction of the reefs by fisherman which lowers the reef profile; and work on the land increasing sediment flow into the bay has slowed the reproduction and increased the death rate (choking out the population by covering it with sediment). At greater depths, the water flow rate is also slower, which translates into higher sediment deposition rates, which combined with lower profile oysters causes severe damage to the oyster population. This is a major problem not only for the oyster industry and for conservation of the species, but the native oyster species also plays an important role in the construction of the underwater ecosystem on which many other species depend.

Many restoration efforts have been made, especially starting in the 1990s, driven largely by the formation of the Oyster Recovery Partnership in 1993 whose primary goal is the recovery of the Chesapeake Bay oyster population by providing knowledge, resources, funding, and personnel to support the recovery. Unfortunately, many of these efforts have found only minimal success. An experiment performed by the Army Corps of Engineers in 2004 used oyster reefs that were constructed at two different heights, and results suggested that the taller reefs favored successful restoration far more than the short, low-relief reefs. This experiment provides the fundamental basis for the model reconstruction we will pursue.

Background and Hypothesis

Before getting into the details of the model and system, it is important to have an understanding of what an oyster reef even is. Oysters are a bivalve mollusc, small sea-dwelling creatures inside of a hard shell. Oysters grow in larger collections called reefs, which are essentially just piles of live oyster shells, dead oyster shells, and sea sediment which collect to form a dynamic underwater biological system. These reefs are frequently harvested for human consumption, but are also struggling for species survival, pitting conservationists against harvesters; however, research has suggested that leaving some stable, unharvested reefs in tact will provide an important habitat for the prosperity of oysters as well as other marine life which dwell in these reefs, ultimately benefiting both the conservationists and the harvesters. These reefs are frequently categorized by their vertical height, called relief. High relief reefs have a large height, while low relief reefs are short. With this basic understanding of the structure of oyster reefs, the details of the system can now be discussed.

The stability of natural habitats often show great sensitivity to minor changes in their surroundings. This is the phenomenon of bi-stability. Sometimes changing the conditions of an ecosystem slowly within a critical parameter range can result in a drastic change, yet returning to the original conditions externally does not force the system to revert back to its previous state. Other mollusk populations have exhibited this bi-stability, including the horse mussel in New Zealand and the blue mussel along the Atlantic Coast of North America. Combining this awareness of bi-stability with the experiment by the Army Corps of Engineers, we will be developing a model of oyster populations and analyzing these mathematical models for bi-stable equilibrium points with varying initial conditions and parameters. The hypothesis we wish to test using this model is that reefs of higher relief are more likely to have stable, living populations of oysters due to an offset of heavy sedimentation and improved disease resistance relative to lower height reefs which will degrade and die in just a few years.

Reproduction of the model from the original paper was successful. There was a slight quantitative discrepancy between the reproduction and the original, but after significant analysis it was determined that there was an error in the given parameter values. The qualitative analysis holds and matches. Based on the model we created, it was shown that our hypothesis is accurate: reefs which start with a higher relief will have a better chance of long-term health and stability. This informs the restoration efforts that in an effort to successfully recover the oyster population, conservationists should concern themselves with constructing tall, high relief reefs, rather than broad, low-relief reefs, using a mix of live and dead oyster shells for the reef's structure.

After the model presented by Jordan-Cooley, et al. was successfully reproduced and analyzed with respect to our hypothesis, a variation of the model was produced to incorporate seasonal weather changes. Previous research had been done developing a periodic function for the temperature over the course of a year using a sine function The temperature of water was also related to the growth rate of an oyster by using a power function of the temperature By combining the temperature function and the growth rate function, we were able to integrate periodic temperature oscillations into the model. This ultimately showed that with seasonal variation, the oyster populations required even more acutely optimized parameters and initial conditions to prosper. From a restoration perspective, this means one would need an even higher relief reef and would need to minimize the sediment deposition rate as much as possible.

Model Description

The initial model was based directly on the equations provided by Jordan-Cooley, et al. The model was adjusted slightly to accommodate the seasonal temperature variations for the model extension. Each of these models is described below.

Original Model

State Variables and Parameters

The oyster reef system has three state variables which are used to define the system. These variables are described in Table 1, below.

Table 1 - The state variables, their units, and their descriptions.

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

In the initial set of equations, there are eleven parameters which are used. These parameters are given below in Table 2.

Table 2 - The parameters, including units, description, and typical range of values.

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 below 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 Input 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.

Table 3 - Parameter values used for initial model simulation.

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 do not match the response shown in the original paper. We did a detailed check of our equations and parameters relative to the original paper and found no differences. I then contacted the corresponding author, Junping Shi, to ask for his input. Ultimately it was determined that the error was somewhere in the list of parameter values given in the original paper, but that the essence of the paper really lay in the qualitative understanding of the system. Since our qualitative responses matched, we decided to move on and accept the quantitative differences.

The final portion of the paper we had to reproduce before expanding upon the original paper was to create bifurcation diagrams of each of the state variables. By solving the system of equations for equilibrium points given a specific set of parameters, we were able to get three equilibrium points for each state variable. Adjusting the sediment deposition rate, parameter C, we were able to repeat this process and get another equilibrium point. By using Mathematica to repeat this process on a large scale, we were able to find and plot the equilibrium points versus the parameter C to reproduce bifurcation diagrams. This will be shown in the results section to follow.

Extended Model

As described in the introduction, we used prior research on the influence of seasonal temperature variations on oyster growth as a basis for a slight modification of our model. The new model incorporated a growth rate scaling function to force oscillations of the growth rate between 0.7 and 1.3, the standard range given in Table 2 above.

State Variables and Parameters

The state variables describing the system do not change with the addition of a growth rate scaling term. The system is still defined by the volume of live oysters (o), dead oysters (b), and sediment (s), as described in Table 1 above. The eleven parameters given in Table 2 above also still apply to the revised system. However, there are two additional parameters for the new system. These are given below in Table 4.

Table 4 - Additional parameters needed for revised model.

Parameter	Units	Description	Value
a	1 / Celsius	Pre-exponential coefficient	0.00032915
z	unitless	Exponent to the power function	2.36

Assumptions

The extension to the model relies on a few key assumptions and simplifications. First and foremost, it uses actual numerical data for the year-round temperatures of the North Inlet Ecosystem in Georgetown, South Carolina. It is therefore assumed that the seasonal swings in temperature in Chesapeake Bay are similar in magnitude and shape to those in South Carolina. Given the proximity of Virginia (where Chesapeake Bay is located) to South Carolina, this assumption is fair. Furthermore, it assumes that the year-to-year changes in temperature do not differ substantially, and therefore can be modeled using a periodic function with a period of one year. This assumption is also reasonable. Lastly, the model assumes that all other factors influencing growth rate are optimized so that when temperature is at a max, the growth rate is at its max (1.3) and when temperature is lowest, the growth rate is at its lowest (0.7). This simplification is reasonable. It is not clear exactly what other factors may influence growth rate, nor how one would model such a dependency.

Model Input and Descriptions

Most of the equations remain unchanged. The three secondary functions, equations 4, 5, and 6 above, are all exactly the same as they were for the original model. The differential equations for the dead oyster volume and sediment volume, equations 2 and 3 above, are also unchanged. The differential equation for the live oyster population, however, has a scaling factor in the first term, as shown below in Equation 7.

Equation 7:

$\frac{dO}{dt} = T(t) 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 )$

The new term, T(t), is the scaling function based on the annual temperature fluctuations. This scaling function is defined by Equation 8, below.

Equation 8:

$T(t) = a Temp(t)^z$

In this function, parameters a and z are the parameters defined in Table 4 above. These values are selected by fitting the expression to experimental data and were adjusted slightly so that when T(t) was multiplied by parameter r, the product oscillated between 0.7 and 1.3. The function Temp(t) is the function which directly relates temperature to time over the course of a year. This relationship is define by Equation 9, below.

Equation 9:

$Temp(t) = 23.795 + \left\vert 7.14 \sin( \pi t) \right\vert$

Once again, the values of this equation, as well as the absolute value used to the right of the plus sign, are all based on fitting the equation to experimental data. These equations, when combined with the rest of the original system, exhibit the oscillatory behavior expected, and influence the basins of attraction of the bifurcation diagram.

Simulation Method

The simulation of the revised model was done in the same way as the original model. NDSolve was used in Mathematica to numerically solve the system of differential equations. The solution to the system was plotted, showing each of the state variables versus time. In this case, there were no pre-existing figures to compare the results to, but the plots showed roughly the same qualitative behavior, just with oscillations in the long-term behavior and changes in the range of sustainable parameters. The results are detailed below.

One aspect of the simulation that was distinctly different was the equilibrium point analysis. Because the solutions had continuous oscillations over time, there was no true equilibrium "point", but rather a limit cycle. This meant that the method used previously for analyzing the long-term stability and equilibrium values would not work. Instead, when creating a plot which allowed for manipulation of the parameters, we included a horizontal line whose position could also be manipulated. By adjusting the line so that it fell in the middle of the limit cycle, an approximate average value for the limit cycle, we could describe the system in a quantitative fashion with reasonable accuracy.

Results

In recreating the model of oyster reef populations based on the paper by William Jordan-Cooley, et al. we used the numerical simulation and bifurcation analysis of the system to try to match our results to the original work. Once we had the model up and running as described in the Model Description, we then moved to performing numerical simulations using Mathematica with parameter values given in the original paper. The details of this computation and simulation are provided in the Mathematica notebook attached at the beginning of this paper and again here for convenience. Deeper biological insight regarding these results will be provided in the Discussion session to follow.

Numerical Simulation

Our Results

Using NDSolve to solve the system of equations and applying the parameters shown in Table 5 below, we were able to create plots of the state variables (O = live oyster volume, B = dead oyster volume, and S = sediment volume) versus time.

Table 5 - Parameter values taken from table 3 in the original paper

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

Using the parameters in Table 1 and initial conditions of O(0) = 0.01 and S(0) = 0.01, Figures 1-4 were produced. In each figure, the first image (left) shows the numerical simulation result I obtained. The second image (right) shows the numerical simulation result from the original paper using the same conditions.

Figure 1 - B(0) = 0.2

Figure 2 - B(0) = 0.1

Figure 3 - B(0) = 0.12

Figure 4 - B(0) = 0.11

Comparison to Original Paper

Clearly the plots in Figures 1-4 above show that our specific simulations do not match the values from the original paper. This was a point of major concern. First we addressed the issue by meticulously walking through our code, comparing all of our equations to the original paper and comparing our code against one another. Still no discrepancies could be found, so the last likely cause for error that came to mind was that the parameter values defined in the original paper were not exactly accurate. Even just small adjustments in the parameters make a substantial difference. I contacted the corresponding author from the original paper, Junping Shi, to try to clarify the issue. Mr. Shi responded saying that they emphasized creating a qualitatively accurate system, so while he could not confirm the parameter set used for the original analysis, he said that if I could reveal the same qualitative behavior in my system to the original paper, then that was sufficient. Therefore I made a minor adjustment to the growth rate, increasing r from 1 to 1.2. This minor change revealed the qualitative bistability described in the original paper. The new numerical simulations are given in Figures 5-7. The new complete parameter set is given in Table 6, below.

Table 6 - Parameter values adjusted to reveal bistability

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

Figure 5 - B(0) = 0.2

Figure 6 - B(0) = 0.168

Figure 7 - B(0) = 0.16

Analysis and Connection to Hypothesis

The plots in Figures 5-7 now show the bistable nature of the system. All it took was a change of r from 1 to 1.2, indicative of just how delicate the system really is. While the numerical values do not match exactly, these plots are qualitatively similar to the numerical simulations from the original paper.

These plots provide information necessary to make a conclusion regarding our hypothesis. The hypothesis was that a taller, higher-relief reef would thrive more than a short reef because it would be more resilient to the sediment deposition. The volume can be seen as the height times the area of the sea floor covered. Therefore, a higher volume reef covering the same unit area will be a taller reef. This corresponds to the initial condition B(0) being larger. Figure 5 shows that for the parameters specified in Table 2 and B(0) = 0.2, the live oysters will ultimately thrive and achieve a steady state equilibrium at a substantial total volume. Yet as the reef height decreases, meaning B(0) decreases, the system crosses a critical threshold and the live oysters die out completely leaving nothing but sediment. Therefore, this simulation suggests that our hypothesis is correct. Reefs that start taller favor the long-term livelihood of the oyster population.

Bifurcation Analysis

Our Results

With the numerical simulations worked out, we then moved on to look at the long term behavior of the system relative to the specific parameter C, the maximum sediment deposition rate. First we analyzed the original model, prior to the extension. In the original paper, the authors created bifurcation diagrams to graphically display the equilibrium behavior of the system. Figures 8-11 show the bifurcation diagrams of each state variable relative to parameter C. The original paper still claims to be using the parameter set from Table 1 above. However, based on our analysis of the numerical simulations, we have adjusted to parameter set to elicit the desired qualitative behavior. The parameter set we used is shown in Table 7 below. *Note that once again our plots are on the left and the plots from the original paper are presented on the right.

Table 7 - Parameter set used to create the bifurcation diagrams in Figures 8-11

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

Figure 8 - Bifurcation diagrams of Live Oyster Volume

Figure 9 - Bifurcation diagrams of Dead Oyster Volume

Figure 10 - Bifurcation diagrams of Sediment Volume

Figure 11 - Bifurcation diagram of Sediment Volume with adjusted scale

Comparison to Original Paper

Using our slightly adjusted parameter set, we end up with bifurcation diagrams which match those of the original paper quite precisely. This further confirms my belief that there was just a small error in their published parameter set which is to blame for our discrepancies. Qualitatively, our system matches the original paper very well. Our bifurcation diagrams were adjusted to show the stability of each region, and we also included the biologically obvious equilibrium of complete death of the oysters in which case sediment just accumulates in their wake.

These findings show that at a specific parameter set, increasing the height of the system, which is to say increasing the initial values of live and/or dead oyster shells, will push the system closer to long-term stability, supporting the findings from the numerical simulation. What's more, it also provides insight into the role of sediment deposition rate. As the sediment deposition rate increases, the conditions under which the system will ultimately prosper become more and more restrictive, eventually reaching a point where no matter what the initial conditions are, there is simply no stable equilibrium other than a dead reef. This means that it is extremely important from a conservation perspective to monitor the sediment being pushed into the bay by human activity on the shore.

Revised Model

After complete analysis of the original model as described by William Jordan-Cooley, et al. we began simulation and analysis of the model with the seasonal temperature changes incorporated. Using this model and the parameters given in Tables 4 and 6, numerical simulations were completed. These simulations are shown below in Figures 12, 13, and 14 below. In each figure, O(0) = S(0) = 0.01, as before.

Figure 12 - Extended model with B(0) = 0.2

Figure 13 - Extended model with B(0) = 0.275

Figure 14 - Extended model with B(0) = 0.28

The basic shape of these responses are the same as for the original model, but there are a few key differences. First, it is clear to see that the system is oscillating, in time as a result of the oscillations in temperature. These oscillations have several repercussions. First, the basin of attraction for the living (non-zero) equilibrium has shrunk. With all of the parameter values exactly the same, the original system's critical value of B(0) was at 0.167452. However, the modified system with seasonal oscillations has a much higher critical value at B(0) = 0.277519. This means that the initial volume of dead shells has to be 1.66 times as large. A difference of this magnitude will be substantial in determining the well-being of a reef.

In addition to the biological ramifications, the new system does not truly have an "equilibrium point" for the surviving system. Instead, with the infinite oscillation, the system ends up in a limit cycle. Therefore, the method of finding equilibrium points and plotting them to develop a bifurcation diagram does not work. Instead the system's long-term behavior must be determined by using the numerical simulation. The stable equilibrium points can be determined using the the long-term solutions to the numerical simulation. Unfortunately, due to the somewhat complex three-dimensional nature of the system, the location of the saddle points, which are the most interesting aspect of the equilibrium diagrams, cannot be determined by simply using the numerical simulation.

Discussion

Model and Hypothesis

The results (shown in the results section) indicate that the larger the initial volume of dead oyster shells, the more sustainable an oyster reef will be. Specifically, there is a critical volume below which the oyster reef will die and above which the oyster reef will go to a stable, living equilibrium. This bistability can be seen in the numerical simulations as the initial condition, B(0), is varied. A very small perturbation can cause a drastic, qualitative shift in the long-term behavior of the system. The same bistable nature can be seen in the bifurcation diagrams, where changes both in the initial conditions of any of the state variables or changes in the parameter value C (sediment deposition rate) will cross from one basin of attraction to another. This hyper-sensitivity only applies to the boundary region between the basins of attraction. At extreme ends of the conditions, the system will either be doomed to failure or definitively prosperous. The effort, then, would be to take systems which are small and developing, systems which may be at risk of death, and try to force these systems toward a more secure, stable regime by cleaning away sediment, adding oysters to the reef, or changing any other environmental factors available for that system (isolate or protect from predators, for example). The hypothesis was that a taller, or higher-relief, reef would be more likely to stay alive long term. The volume of the reef is related to the height of the reef by dividing the volume over the area on the sea floor which the oyster reef inhabits. So larger volumes (over roughly the same area) translate to higher-relief reefs. Therefore, our model results suggest that our hypothesis was correct: high-relief reefs are more likely to reach a stable, living equilibrium.

This analysis, on a qualitative level, applies to both the model as described in the original paper as well as the adjusted model which accounts for seasonal temperature changes. However, it can be seen that as the temperature fluctuates, the net effect is to put more strain on the oyster reef system, tending to pull the system toward the stable equilibrium at which both live and dead oyster shells go to zero. This means that in regions where the climate is more unpredictable, or if there is a particularly cold winter one year, restoration and conservation efforts may need to go to greater lengths, building up restored reefs as high as possible. This also has implications with global climate shift. If the average temperatures begin to decrease in the bay, the oyster growth rate will decrease, and the oyster population may rapidly approach endangerment.

Model Limitations

The model that we have created based on that given by Jordan-Cooley, et al. has some inherent limitations. First of all, the numbers that we have are not directly real-world values. For example, the volumes of sediment, live oysters, and dead oysters show the response of the system for hypothetical values, but to apply this model to a specific, real oyster reef, all of the parameters and initial conditions would have to be measured and scaled to fit the system. What we learn from our model is really more of a qualitative analysis. Adjusting a parameter here will have this effect, or changing an initial condition in this way would have this effect. Using this knowledge, we could determine an ideal scenario, such as maximizing the initial reef height, minimizing sediment deposition rate, etc., and then try to impose these optimal conditions on a real system to the best of our abilities in an effort to conserve or reconstruct oyster reefs.

Furthermore, the model that we have does make several assumptions, as detailed in the model description, which are all biologically founded, but as a real system is never ideal, there will inevitably be some loss of accuracy and precision with these assumptions.

Lastly, there are some additional details and factors which the model ignores. Some examples include a more detailed geometric description of each individual reef (something more hill-shaped rather than a rectangular prism perpendicular to the sea floor)as well as the geometry and influence of neighboring reefs. The model treats each reef as a completely independent, free-standing system, which in a bay with a substantial concentration of reefs may not be entirely accurate. There are also many other factors which may influence the oyster prosperity to varying degrees, including the seasonal climate changes, salinity or chemical composition of the water (including pollution), and harvesting for human consumption.

Discrepancies With Original Paper

As shown and described in the results section, the specific numerical results obtained using the exact equations and parameters of the original paper do not quantitatively agree. After extensive checking of our work and email contact with the corresponding author, we have determined that the likely source of the discrepancy is in the parameter values given in Table 3 of the original paper. Unfortunately, this table of parameter values is used for all other analysis in the original work, making it impossible to duplicate exactly what the authors describe. However, we have results that match quantitatively with everything described in the paper, including both numerical simulations and bifurcation analysis. Since, as described in the Model Limitations above and confirmed by the corresponding author, the purpose of the model is to gain a qualitative understand of the effects of certain conditions, the quantitative discrepancy is not a major concern. We have noted the discrepancy and researched possible causes, and have determined that our model is still adequate for analysis of the system.

Another point worth noting when analyzing the differences between our simulations and those of the original paper is that when our system dies, the live oyster volume goes to zero, followed by the dead oysters as those shells decay. At that point, there are no oysters remaining, and the sediment reaches a stable point at a very high value. However, looking at the original paper's numerical simulations (see Figures 2 and 4 above), there is a qualitative difference. The live oyster shells go to zero, however the dead oyster shells reach a stable equilibrium at a non-zero volume. Biologically this does not make sense because the oyster shells would decay with time, and with no live oyster population to die and feed the dead oyster volume, the dead oyster volume should go to zero. Mathematically this makes sense as well. See Equation 2. The first two terms in this equation have O(t), which goes to zero. Therefore these terms go to zero. This just leaves the linear degradation term. The dead oyster volume should decrease and approach 0. Both biologically and mathematically, the dead system should become a pile of sediment and nothing else, as my model shows. The original model seems to have a small volume of dead shells remaining at all times. While I cannot find any clear source for this error, I am certain that this is indicative of an error in the original paper, and may be the true source of discrepancy in our quantitative analysis (as opposed to a mistake with the parameter values). I have not contacted the corresponding author yet, but may in the near future to see if he has any explanation.

Relationship to Literature

Hypothesis

As mentioned, our model shows that high-relief reefs should be more stable and sustainable long-term than low-relief reefs, which will tend to die out. This has been confirmed experimentally in some of the restoration efforts of the past decade. It was found in a study of reconstructed reefs, some high-relief and some low-relief, in the Great Wicomico River (Virginia) that "oyster density was fourfold greater on high-relief than on low-relief reefs." This agrees with our findings through the model reconstruction and confirms our hypothesis.

Model

In the construction of the model, several assumptions were made and each term in the equations was meant to represent a specific biological phenomenon. This is all detailed on the Model Description page. The basis for the model comes from a rough biological understanding and the detailed research that already exists in the field. The logistic growth of live oysters was described by Schulte, et al. The degradation of dead oyster shells, resulting in a decrease in the volume of dead oyster shells, was described by Smith et al. The influence of reef height in the water column and filtration by the live oyster volume on the sediment volume, as well as the influence this sediment volume has on the live oyster population, was described by Lenihan, et al. and Jordan These sources all provide a background in the literature for the assumptions and model construction used in the original paper.

Conclusion and Future Work

The models constructed in this paper provide important insight into the restoration and conservation efforts of oyster populations. These results may also indicate an important balance necessary for the long-term stability of the oyster harvesting industry for human consumption. Ultimately it seems that it is important to maximize the height of live and dead oysters in a developing reef and to minimize the sediment volume present in order to support long-term health of the oyster reef. Furthermore, the bifurcation analysis and manipulation of the numerical analysis show the important role other environmental factors play, including the sediment deposition rate. While this is not as easily controlled, near-shore construction or manufacturing which could pollute the water may be very detrimental to the survival of oyster reefs. The influence of temperature was also seen in the extended model. Major climate shifts, or perhaps even thermal pollution of the water by local industry, can also have substantial impacts on the well-being of oyster reefs. A decrease in water temperature of only a few degrees may shift an oyster reef from one basin of attraction to another, an issue of life or death.

Moving forward, there is still plenty more to learn about the dynamics of oyster populations. All of the analysis done within this paper has treated each oyster reef as a completely independent system. In reality, there are likely to be many oyster reefs in close proximity to one another, especially in a relatively small body of water like the Chesapeake Bay. A study of the dynamics between neighboring reefs may be very revealing. For example, if a few large reefs are maintained as no-harvesting zones, can those reefs provide the stability, protection, and offspring necessary to sustain other reefs nearby, which can be harvested more aggressively without permanent damage? Performing bifurcation analysis of some of the other potentially variable parameters may provide a better understanding of how important parameters like the mortality due to predation and disease or the rate of sediment erosion influence the health of an oyster reef? Once that information is detailed, then restoration efforts can try to accommodate those factors, isolating new reefs from predators or putting up barriers to minimize sediment erosion.

As more is learned about the dynamics of oyster systems, efforts to preserve the species can be more precisely oriented to optimize resources, with the ultimate goal of establishing a substantial and robust population which can sustain itself.

Appendix - Mathematica File

All of the programming in mathematica is compressed into a single file. Annotations within the file describe what each section of the code does. The code can be found at the link below:

Bsp6_Final_File