Exemplary Results Draft 1

Author: Bradley Potter

Modeling Results

In recreating the model of oyster reef populations based on the paper "Bistability in a differential equation model of oyster reef height and sediment accumulation" 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 will be provided in the Mathematica notebook file with the final paper.

Numerical Simulation

Our Results

Using NDSolve to solve the system of equations and applying the parameters shown in Table 1 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 1 - 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.2	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 2.

Table 2 - 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.

What do these results say about our hypothesis? The hypothesis was that a taller, higher-profile 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. 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 3 below. *Note that once again our plots are on the left and the plots from the original paper are presented on the right.

Table 3 - Parameter set we 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 do not relate directly to the hypothesis regarding the height of the reef relative to its long-term success. However, the bifurcation analysis may help provide a direction for the model extension to follow.

Natural Extensions of Model

Having seen the influence of changing B(0) and C in the system, it naturally leads to questions like "what about the other parameters?" and "what about the other initial conditions?" These are very valid questions which, with minor adjustments to the Mathematica code, can be addressed. Perhaps we will find that if we can collect a sufficiently large volume of living oysters and we add them to a system to increase O(0), then the system will flourish. Investigating these concepts can provide the basis for new conservation and restoration methods.