%\documentclass[]{article}%proc
\documentclass[12pt,english]{article}
\usepackage{geometry}
\geometry{verbose,letterpaper,tmargin=2.54cm,bmargin=2.54cm,lmargin=2.54cm,rmargin=2.54cm} %can cancel
\usepackage{setspace} %can cancel
\usepackage{lineno}
\linenumbers
\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\renewcommand{\thecorollary}{\thesection.\arabic{corollary}}
\renewcommand{\thelemma}{\thesection.\arabic{lemma}}
\renewcommand{\thetheorem}{\thesection.\arabic{theorem}}
\newtheorem{proposition}{Proposition}
\renewcommand{\theproposition}{\arabic{proposition}}
\renewcommand{\theequation}{\arabic{equation}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\ba}{\begin{eqnarray*}}
\newcommand{\ea}{\end{eqnarray*}}
\newcommand{\dsize}{\displaystyle}
\usepackage{graphicx}
\usepackage{epsfig}
\usepackage{color}
\usepackage{natbib}
\bibpunct{(}{)}{;}{a}{}{}
%this is cool
\newcommand{\nagy}[1]{\textcolor{blue}{\textit{#1}}}
%\doublespacing
\begin{document}
\begin{spacing}{1.5} %can cancel
%\begin{flushleft}
\title{Modeling the population dynamics of lemon sharks}
\author{Easton R. White$^{1,2,\dag }$, John D. Nagy$^{2,3}$ and Samuel H. Gruber$^{4}$}
\maketitle
{\footnotesize 1. School of Life Sciences, Arizona State University
2. Department of Life Sciences, Scottsdale Community College
3. School of Mathematical and Statistical Sciences, Arizona State University
4. Bimini Biological Field Station, Bimini, Bahamas
$^\dag$ Current address: Department of Biology, University of Victoria, BC, Canada}
\bigskip
* Corresponding author: E.R. White, Email: eastonrwhite@gmail.com
%\tableofcontents
\begin{abstract}
Long-lived marine megavertebrates (e.g. sharks, turtles, mammals, and seabirds) are inherently vulnerable to anthropogenic mortality. Although some mathematical models have been applied successfully to manage these animals, more detailed treatments are often needed to assess potential drivers of population dynamics. In particular, factors such as age-structure, density-dependent feedbacks on reproduction, and demographic stochasticity are important for understanding population trends, but are often difficult to assess. Lemon sharks (\emph{Negaprion brevirostris}) have a pelagic adult phase that makes them logistically difficult to study. However, juveniles use coastal nursery areas where their densities can be high. Thus, we use a stage-structured, Markov-chain stochastic model to describe lemon shark population dynamics from a 17-year longitudinal dataset at a coastal nursery area at Bimini, Bahamas. We found that the interaction between delayed breeding and demographic stochasticity accounts for 33 to 49\% of the variance. Demographic stochasticity contributed all random effects in this model, suggesting that the existence of unmodeled environmental factors may be driving the majority of interannual population fluctuations. In addition, we are able to use our model to estimate the natural mortality rate of older age classes of lemon sharks that are difficult to study. Further, we use our model to examine what effect the length of a time series plays on deciphering ecological patterns. We find that--- even with a relatively long time series---our sampling still misses important rare events. Our approach can be used more broadly to infer population dynamics of other large vertebrates in which age structure and demographic stochasticity are important.
\end{abstract}
%% main text
\section{Introduction}
Many large marine megavertebrates (e.g. sharks, turtles, mammals, seabirds) are particularly vulnerable to anthropogenic mortality due to their complex life history characteristics, including long lifespans, delayed maturity, low fecundity, and extended migrations \citep{Fujiwara2001, Baum2003, Lewison2004, Ward-Paige2012,Senko2014}. These animals often act as ecological keystones, and their removal can lead to considerable ecosystem changes such as cascading ecological effects on lower trophic levels \citep{Lewison2004, Myers2007, Heithaus2008,Wirsing2008,Baum2009,Ferretti2010,Heithaus2010}. For example, as predators, sharks not only regulate their own prey populations but also those of species deeper in the food web \citep{Myers2007,Baum2009,Ferretti2010,Heithaus2010} and see recent review by \citep{Heupel2014}. Given their importance to ecosystem stability and the multiple anthropogenic threats they face \citep{Kyne2012,Worm2013}, it is imperative that we develop a better understanding of shark population dynamics, particularly to identify primary drivers of annual population variation.
Physiologically structured population models \citep{Crouse1987,Caswell2001,MorrisDoak2002,BrauerCastillo-Chavez2012} that incorporate delayed breeding \citep{Gourley2004, Wang2009}, density-dependent mechanisms \citep{NeubertCaswell2000,Caswell2001}, demographic stochasticity \citep{MorrisDoak2002,OvaskainenMeerson2010, Jenouvrier2012,Mills2012,McCarthyPossingham2012}, or some combination of these processes, have been applied to many ecological systems to answer questions related to population dynamics, conservation, and management. In examining shark populations, physiologically structured discrete demographic models have been used to study overfishing and population viability, calculate specific demographic parameters, and predict population dynamics \citep{Hoenig1990,Cortes1998,Gallucci2006,Gedamke2007,Mollet2002,Knip2010,Cortes2002a,Beerkircher2002,Booth2011,Tsai2010,Forrest2009,Dulvy2010}, also see review \citep{Cortes2007}. Although these demographic models are useful, they typically have three key shortcomings \citep{Cortes2007}: (1) they typically include assumptions that are biologically unrealistic, including density independent, deterministic mechanisms, which makes application of model outputs to real data difficult to accept; (2) given that these models are deterministic, they cannot capture demographic stochastic events, which are likely to be important drivers of interannual population fluctuations; and (3) parameterization and validation of such models from data are often logistically difficult and require long-term field operations \citep{Lewison2004,McCauley2012}.
The third challenge has been met by a longitudinal field study of lemon sharks (\emph{Negaprion brevirostris}) at Bimini, Bahamas. Data from this study includes an annual population census of juvenile lemon sharks (ages 0-2 years) from 1996 to the present. The number of juveniles in our study population (see Methods) typically fluctuates between 50 and 100 sharks, although the complete range is estimated to be between about 35 and 150 (Fig.\ \ref{biminisharks}), which illustrates the significance of annual fluctuations in the juvenile age class of this lemon shark population. Also, fecundity and mortality rates have been estimated precisely using mark-recapture and genetic methods \citep{Gruber2001,Feldheim2002, Feldheim2004}.
The causes of annual variation in population size remain unclear for many species, and we are unaware of any previous studies that have assessed these causes in detail for lemon sharks. Furthermore, little is known regarding mortality rates of both the larger juveniles (ages 3-11), who leave the nursery site around age three, and the adults (ages 12+) who mature at approximately 12 years of age \citep{Brown1988,Kessel2010}.
Here we present a mathematical model detailing annual fluctuations in a juvenile lemon shark population at Bimini, Bahamas. The model is physiologically structured, with age class as the (discrete) structuring variable. Hoenig and Gruber (1990) also used a physiologically structured model in the form of a Leslie matrix. Unlike our model, their work only deals with deterministic models. We make births and deaths stochastic, but fix environmental parameters. We parameterize the model using estimates of fecundity and juvenile mortality rate estimates obtained from the field study described above and apply inverse pattern-oriented techniques to fit the model to data \citep{Wiegand2003,Grimm2005,Hartig2011,Anadon2012}.
We show that demographic stochasticity in the model predicts only 33\% to 49\% of the observed variance. Therefore, we predict that another source of stochasticity, probably environmental, accounts for at least half of the variance observed in annual population fluctuations in this population of lemon sharks. In addition, our use of inverse pattern-oriented modeling allows us to estimate the unknown mortality rates of subadults and adults, which illustrates the utility of this type of modeling. Further, such modeling can be used to study how sample size (length of the time series) affects estimated population parameters and dynamics. We find that time series with 15 consecutive years may be too short to capture critical, but rare, stochastic events.
\section{Methods}
\subsection{Study Site and Field Data}
This study builds on field work conducted in Bimini Lagoon, Bimini, Bahamas (25$^o$44N, 79$^o$16W). The Biminis are located approximately 86 km east of Miami, Florida and provide habitat for numerous species of fish, arthropods, birds, and mollusks \citep{Jennings2012}. Of the three lemon shark nursery sites (as defined by \citep{Heupel2007}) in Bimini, our study focuses on the most northerly one, known as the North Sound. Between 1996 and 2012, standardized gillnet methods were used to capture juvenile lemon sharks within 45 days of parturition (Fig. \ref{biminisharks}). For a more detailed treatment of the gillnetting protocols and yearly censuses, see \citep{Manire1993,Gruber2001,Feldheim2004}; and \citep{DiBattista2011}.
In addition to population censuses, genetic analyses from tissue samples were used to reconstruct family pedigrees \citep{Feldheim2002, Feldheim2004,DiBattista2011} from which we estimate per-female annual fecundity in the Bimini population (Fig. \ref{littersize}). Reproductive-age female lemon sharks (ages 12+) show strong philopatry to their natal nursery sites, with about 45\% returning to a given nursery area to reproduce every other year \citep{Feldheim2002}. Newborn and juvenile sharks (ages 0-2 years old) stay in these protected, mangrove fringed nursery areas \citep{Morrissey1993}. In addition, there appears to be very little dispersal among nursery sites in this region, so the population of juvenile lemon sharks in the North Sound is essentially closed \citep{Gruber2001}. At about age 3, lemon sharks enter their subadult phase (ages 3-11), begin to leave the lagoon area and move to deeper waters \citep{Morrissey1993,Franks2007,Newman2010}.
Our model is constructed to capture this natural history in such a way that key model parameters, such as mortality rates of different age groups, can be estimated from field data. This allows us to use inverse pattern-oriented methods to estimate other life history parameters that are otherwise difficult or presently impossible to measure directly.
\subsection{Model}
We model the Bimini lemon shark population as an age-structured, Markov-chain stochastic process. We choose this formalism due to the complexity of the lemon shark's life cycle---in particular the delay in breeding to the 12th year---and because breeding populations at nursery sites in any given year appear to be too small to be buffered from fluctuations due to demographic stochasticity. Since the maximum age for lemon sharks is thought to be 25 years \citep{Cortes1998,Gedamke2007}, we assume a maximum of 26 age classes (including the 0$^{th}$ age class).
Let $\mathbf{x}(n)$ be the shark population vector at census time $n$; that is, its elements, $x_a(n),\, a \in \{0,1,\ldots,25\},\, n \in \{0,1,2,\ldots\}$, represent the number of lemon sharks of age $a$ in the North Sound population, including all animals born to and breeding in the North Sound nursery, whether they are in the nursery or open ocean at census $n$. Age class 0 represents sharks born the year of the census. To match the timing of the actual Bimini census, we assume that this census occurs just after reproduction (i.e., pups are born April/early May and are sampled late May/June).
\subsubsection{Fecundity}
We assume an equal sex ratio and that females only reproduce every other year after their 11th year of life \citep{Feldheim2002}. Let $R$ be a random variable taking on values in $\{0,1,2,\ldots\}$ with probability density $\{p_0, p_1, p_2, \ldots\}$. We interpret $R$ as the number of offspring born to a particular breeding female, and $p_i$ as the probability that a female gives birth to $i$ pups. Let the number of breeding females in year $n$ be denoted $b_n$; that is,
\begin{equation}
b_n := \frac{1}{4} \sum_{a=12}^{25} x_a(n);
\end{equation}
the coefficient of $1/4$ follows from the assumptions of equal sex ratio and biennial breeding with the further assumption that, for each age class, the proportions of females breeding in even and odd numbered years are equal. We assume that all breeding females have the same reproductive potential regardless of age class, time or population density. Therefore, the set $\{R_i; i \in \{1,2,\ldots,b_n\} \}$ is a collection of independent, identically distributed random variables, and $R_i$ is the reproductive output of the $i$th female in year $n$. Therefore,
\begin{equation}
B(n) = \sum_{i=1}^{b_n} R_i
\end{equation}
is the total reproductive output of the population in year $n$. Note that the dependence of $B$ on $n$ comes only through the number of breeding females in year $n$, not through $R$.
The probability density, $\{p_0, p_1, p_2, \ldots\}$, for $R$ can be obtained from a variety of assumptions. We consider two possibilities. In some simulations, we obtain this density from data; in particular, each $p_i$ is set to the observed frequency of females producing $i$ pups (Fig. \ref{littersize}), with the convention that $p_j = 0$ for all $j >18$. In the second case, we assume that all $R_i$s are Poisson-distributed with fixed mean $\lambda$. In this case, the probability density for the total population fecundity in year $n$ becomes
\begin{equation} \label{PBn}
\Pr(\{B(n) = j\}) = e^{-b_n \lambda} \frac{(b_n \lambda)^j}{j!} \mbox{ for all } j \in \{0,1,2,\ldots\}.
\end{equation}
A Poisson distribution is often assumed to be a good fit for a birth process. We explicitly test that assumption here.
\subsubsection{Mortality}
We assume that the probability of mortality is evenly distributed across all individuals in a given age class; therefore, within an age class the number of sharks that die between censuses is distributed binomially. Generally speaking, the parameter of that distribution---the probability that a given shark dies---could potentially depend on population size. However, in the case of lemon sharks, we only have evidence for density-dependent mortality in the first age class \citep{Gruber2001,Gedamke2007}. There is insufficient evidence to support either density-dependent or -independent mortality assumptions in other age classes; indeed, very little is known about lemon sharks once they leave their nursery area. Therefore, as a first approximation we chose density-independent mortality for all age-classes above the first.
In this first age class, the probability that a shark pup dies between birth (age class 0) and its second census (i.e. dies in age class 1) is a generally increasing function of the size of its cohort ($x_0(n)$) in the lagoon in that year \citep{Gruber2001,Gedamke2007}. This type of density-dependent mortality may be a result of reduced prey resources (although the population does not appear to be prey-limited in any way), predation from large barracudas, predation from other shark species, or cannibalism, which has been documented for this population \citep{Morrissey1993,Guttridge2012}. We model this cohort-density dependence with a generalized Michaelis-Menten function (equivalent to a Hill function):
\begin{equation}\label{death_rate}
\hat{\mu}(x_0(n))= \frac{x_0(n)^h}{k^h + x_0(n)^h},
\end{equation}
with (constant) Hill and shape parameters $h\geq 1$ and $k > 0$, respectively. Let $M_0(n)$ be a random variable representing the number of sharks born in year $n$ that die between their first and second censuses. Then $M_0(n)$ has probability distribution
\begin{equation}
\Pr(\{M_0(n) = m\}) = {x_0(n) \choose m} \hat{\mu}^m (1-\hat{\mu})^{x_0-m}, \hspace{2ex} m \in \{0,1,\ldots,x_0(n)\},
\end{equation}
where we suppress the dependence of $\hat{\mu}$ on $x_0(n)$ for clarity.
As a first approximation, we assume that no age classes except the first have density dependent mortality. We further assume that the probability of mortality for any shark in age classes 1 or higher is invariant across individuals regardless of age (this assumption could be relaxed with our model structure). We denote this constant probability as $\mu$ and define $M_a(n)$ to be a random variable representing number of deaths in age class $a \in \{1,2,\ldots,25\}$. Then
\begin{equation} \label{adult_mort_pdf}
\Pr(\{M_a(n) = m \}) = {x_a(n) \choose m} \mu^m (1-\mu)^{x_a - m}, \hspace{2ex} m \in \{0,1,\ldots,x_a(n)\}.
\end{equation}
It is important to note that many of these assumptions can be relaxed without altering the form of our model (see below). For example, here we assume no fishing mortality because there is no shark fishery in Bimini. One could easily incorporate such an assumption into $\mu$, and even make $\mu$ age-class- and (or) density-dependent with fairly obvious alterations to the probability distribution for mortality equation (\ref{adult_mort_pdf}) which have no effect on the overall model form.
\subsubsection{Model Form and Parameterization}
The development above generates a model with the following form:
\begin{equation} \label{model}
\left \{
\begin{array}{rcl}
x_0(n+1) & = & B(n),\\
x_a(n+1) & = & 1-M_{a-1}(n),\hspace{2ex} a \in \{ 1,2,\ldots,25 \},\\
\mathbf{x}(0) & = & \mathbf{x}_0,\\
n & \in & \{0,1,\ldots\},
\end{array}
\right .
\end{equation}
where $\mathbf{x}_0$ is the initial age distribution.
Application of model (\ref{model}) to the lemon shark population requires field estimates of fecundity and mortality. Starting with the former, as noted above we can estimate the probability distribution directly from data (Fig. \ref{littersize}), or we can assume that per-female reproductive output is Poisson-distributed with mean $\lambda$. Data from Bimini over the last 20 years suggests that $\lambda \approx 6.1$ pups per female \citep{Feldheim2002}, although this value is somewhat lower than that used in previous modeling studies (perhaps because of the high mortality of pups between birth and our sampling season; in this case the 6.1 pups per female simply represents the number of sharks that make it past that interim period) \citep{Hoenig1990,Gedamke2007}.
Less is known about mortality in this species. The first-year mortality function, equation (\ref{death_rate}), requires two parameters: the Hill parameter ($h$) and the shape parameter ($k$), whereas non-first-year mortality only requires an estimate of mean per-shark probability of mortality, $\mu$. Because of the lack of data, we compare model output to population data from the Bimini study to define a range of potential values for these parameters using a sensitivity analysis similar to that in \citep{Hartig2011}. We describe this method in the next section.
\subsection{Simulations and analysis}
For clarity of exposition we will refer to sharks from ages 0 to 2 as juveniles, from ages 3 to 11 as subadults and above age 12 as adults. We partition the traditionally-defined juvenile class into two groups to make connections with the Bimini study---``juveniles" as defined above are the animals actually caught each year in the Bimini nursery census. In particular, we evaluate the model by comparing its behavior to the Bimini nursery census data. Model (\ref{model}) was implemented and all analyses were conducted using the open-source computing language R \citep{Rcitation}.
We used an inverse pattern-oriented technique to quantitatively compare simulations and data \citep{Wiegand2003,Grimm2005,Hartig2011,Anadon2012}. In this approach, one compares actual means and variances from data with the distribution of means and variances predicted by model simulations. We explored 9000 distinct parameter combinations ($\lambda$: range 1-15; $k$: range 0-200; $\mu$: range 0-0.35; $h$ fixed at 1). Note that, although $\lambda$ is known from Fig. \ref{littersize}, it is still of interest to examine a range of values for $\lambda$ to evaluate the type of compensatory responses generated by variations in $\lambda$. We fixed $h=1$ in equation (\ref{death_rate}) because this gave the best fit (least sum of a squares) to the relationship between mortality and density obtained by \citep{Gedamke2007}. Likewise, we also limited $k$ below 200 because larger values produce a linear mortality rate that greatly underestimates that measured by \citep{Gedamke2007} and \citep{Gruber2001}. For each of the 9000 combinations, model simulations were repeated 100 times with the same initial conditions. Simulations were run for 300 time steps (in ``years") or until the population went extinct. Among-year mean and variance for juveniles were recorded for each of the 100 simulations, which provided an estimate of parametric distributions of simulation means and variances. To assess how well any given parameter combination represented the field data, we determined if both the mean and variance of the field data set fell within the middle 95 percent of the distribution of means and variances generated from the simulations of a particular parameter combination. The fit between model and data was deemed ``good" if (i) mean size of the juvenile population size from the data fell within the middle 95\% of the distribution of mean juvenile population sizes from the simulations; (ii) variance in juvenile population size from the data fell within the middle 95\% of its distribution from the simulations; and (iii) the population remained extant after 300 years in each of the 100 runs. By eliminating parameter combinations that fail to satisfy any of these three criteria, we effectively constrain possible values for the unknown parameters (see Fig. \ref{alldata}).
\subsubsection{Testing the effect of sampling length}
We generated the mean and variance distributions by drawing random 17-year-long samples out of each of the 100 trials; that is, we sampled a randomly-chosen sequence of 17 consecutive years---after initial transient dynamics have settled down---from each simulation as an analogue to the 17 consecutive years of field data at hand. However, a question arises regarding how well a 17 year data set represents centuries of ecological dynamics. To assess this, we compare data to various sized samples from simulations, including complete (300 year) simulations, and find that interannual variance in population size is strongly affected by the duration of sample run.
\section{Results}
The Bimini nursery data suggest that on average about 77 juvenile sharks inhabit the lagoon at census time, on average, with interannual variance, $s^2$, of 498 (Fig.\ \ref{biminisharks}). With default parameters ($\lambda=6.1$, $\mu=0.15$, $k=100$) we can match mean juvenile population sizes from the simulation to the actual mean; however, at default the mean of variances from simulation runs, denoted $s_{sim}^2$, is 176 ($n=100$ trials), which represents only 35 percent of $s^2$ in the actual data set.
We compared field estimates of annual mean and variance of population size to annual mean and variance of simulated population sizes generated by each of the 9000 parameter combinations. By the methods described in the simulations and analysis section, we determined that most of the 9000 parameter combinations were a poor fit to actual data. Therefore, the volume of the possible parameter space that admits dynamics having any chance of representing the actual Bimini population is greatly constrained (Fig. \ref{alldata}), even when $\lambda$ is allowed to vary from 1-15. Note that in Fig. \ref{alldata}, lemon shark births were Possion distributed according to equation (\ref{PBn}). Using the observed distribution of births per female (Fig. \ref{littersize}) further tightens our constraints on $k$ and $\mu$ (Fig. \ref{halfsatvsmort}).
Our model also places tight constraints on adult mortality ($\mu= 0.14-0.17$). Values for $\mu$ greater than 0.17 drive the population to extinction. This range for $\mu$ also tends to agree with the indirect methods for estimating mortality given by \citep{Pauly1980,Hoenig1983} and \citep{Jensen1996} and summarized in Table 2. These indirect methods place the mortality rate between 0.086-0.179. Interestingly, the half-saturation value ($k$) is much less constrained; good fits can be obtained for any $k>100$ (Fig. \ref{halfsatvsmort}).
It is important to note that even for parameter combinations that qualified as good fits, all greatly underestimated annual variance observed in the actual data.
In general, our model dynamics were robust with respect to the two assumed distributions of per-female fecundity. Specifically, both produced very similar regions of parameter combinations that matched Bimini (Fig.\ \ref{halfsatvsmort}), with the nuanced exceptions noted above. In addition, however, the actual distribution of litter sizes consistently generated a higher variance in annual population size than did the Poisson distribution. Therefore, although the Poisson distribution is a reasonable choice to use when the actual distribution is not available, one needs to be aware that it tends to underestimate variance.
To test the effect of sample size on characterization of population dynamics from real data, we compared sequences of years of various lengths as described in the methods (Fig.\ \ref{vt}). We find that variance is a generally increasing function of sample length and appears to approach an asymptote which represents our estimate of the parametric variance. This asymptote is considerably higher than the variance typically obtained using 17-year sample runs (vertical red line in Fig.\ \ref{vt}). Nevertheless, this estimated parametric variance is still well below the observed variance of 498 sharks$^2$ (horizontal green line in Fig.\ \ref{vt}).
\section{Discussion}
Lemon sharks have complex life histories---they delay breeding for over a decade, mature in an environment (nursery lagoons at the Bimini site) vastly different from their adult habitat (open ocean) and when mature breed every other year. Also, nursery populations are typically not large enough (order $10^2$ at most) to buffer demographic stochasticity; indeed, demographic stochasticity can dominate dynamics in patchy systems with sizes orders of magnitude larger than this one \citep{McKaneNewman2004,McKaneNewman2005,McKaneetal2007}.
Therefore, generalized, deterministic population models can hope to elucidate only the broadest outlines of lemon shark population dynamics and should be interpreted only in the ``ensemble average" sense \citep{McKaneNewman2004,vanKampen1992}. That is, deterministic models at best provide an expectation or mean behavior for an infinite number of Bimini's lemon shark populations. Although this abstract notion of an ensemble mean is sensible and provides some insight about expected behavior of the population, under a suitable definition of ``expected," that insight is limited because such models provide no measure of the fluctuations about this ensemble average one can expect to see in any real instance \citep{OvaskainenMeerson2010}).
We addressed this shortcoming by developing a model of the lemon shark population at Bimini incorporating both demographic stochasticity and age structure. Despite the added realism, the model remains relatively simple. Parameters requiring estimates include the probability distribution for the number of pups born to breeding females in a give year, or just the mean number of pups per female if one assumes a Poisson distribution, two parameters characterizing density-dependent mortality in the first age class, and the probability(-ies) of mortality for all individuals in all other age classes, which here we assume to be invariant across individuals. We obtain the reproductive parameters directly from a field study of the Bimini lemon shark nursery (Fig.\ \ref{littersize}), and we use an exhaustive parameter search to obtain bounds on the other parameters.
\subsection{Fecundity assumptions}
As far as we know, this is the first elasmobranch study to use actual litter sizes derived from genetic data to parameterize a mathematical or computational model's birth function. Typically such studies rely on an assumed distribution, of which Poisson is usually thought to be a good first estimate. Since we have a data-driven, realistic estimate of the distribution of litter sizes, we can explore the consequences of making the Poisson assumption, and we find that, in this instance, the Poisson assumption appears reasonable. If researchers only have mean number of pups per female, which is often the case without genetic maternity data, the Poisson distribution approximates the actual distribution remarkably well (Fig.\ \ref{littersize}). Also, the regions of viability within the parameter space for both actual and Poisson distributions are comparable (Fig. \ref{halfsatvsmort}). However, the simulations we ran using the Poisson assumption consistently exhibited lower interannual variance than did simulations using the actual reproductive data. Therefore, although the Poisson assumption generates estimates of adult mortality that are essentially identical to those produced using the data, it should be used with care when modeling to assess population viability and fluctuations.
\subsection{Mortality assumptions}
In the present study we model density-dependent mortality in the first age class as a nonstationary Bernoulli process; that is, the probability of mortality is a generally increasing function of the number of sharks in this age class. This assumption is justified based on field data \citep{Gruber2001,Gedamke2007}. We represented this density-dependent mortality using a Hill function (equivalently a Michaelis-Menten form) from phenomenological considerations only---we hypothesize a monotonically increasing function, which the Hill function exhibits, with the added benefits of relative simplicity and plasticity. Nevertheless, the parameters of this function have biological meaning. In particular, the shape parameter, $k$, measures how sensitive sharks are to competition from same-age conspecifics. As such, although it quantifies an important biological response, in general it will be very difficult to estimate accurately in the field. Importantly, this model demonstrates that practical estimates need not be very precise. The model matches data for a wide range of this parameter's values; in fact, for $k > 100$ or so, the model fit is largely unaffected (Fig. \ref{halfsatvsmort}). Therefore, this model is robust for parameterizations of $k$.
The bounds we found for subsdult and adult mortality rates are remarkably narrow (Figs.\ \ref{alldata} \& \ref{halfsatvsmort}). In fact, a key prediction from this model is that the probability of mortality for any individual in any age class, after the first, lies between 0.14 and 0.17. Indeed, above 0.17, populations invariably die out very rapidly. These mortality estimates compare favorably to estimates obtained using the techniques in Table 2. We emphasize that the technique applied here, known as inverse pattern-orientated modeling, has seldom been used to estimate parameters for any shark population. If subadult and adult mortality were to increase by only a few percentage points, the model predicts rapid extinction. Therefore, we suggest that any added fishing pressure (there is currently no fishing pressures for lemon sharks at Bimini) to this population would threaten its sustainability. This result agrees with others in suggesting that long-lived species with low fecundity, like lemon sharks, would not be able to handle added fishing mortality in adult age classes \citep{Ward-Paige2012}.
One difficulty this model faces is a lack of information about mortality in subadults and adults. As a first approximation, we assume a fixed probability of mortality in all age classes after the first. However, we recognize the tentative nature of this assumption. In particular, it is challenging to relate age to mortality rate of untagged adult sharks. \citep{Peterson1984} suggested that one method to overcome this problem is to construct a function that maps shark mortality rate to mass. In this case age can then be related to mass with conversion methods such as those used in \citep{Beerkircher2002}. This method needs modifications in situations where mortality is density-dependent, and it relies on assumptions of von Bertalanffy growth and reliable estimate of mass, length, and age from catch data. However, where applicable, this technique may only be needed for the first few age classes, or at least until the age at which individual sharks are large enough avoid being preyed upon by natural predators \citep{Cortes1998}. However, in this study such a technique is unavailable because size data of adult sharks is limited to length only; length-to-weight standards exist only for juvenile age classes (Gruber, unpublished data).
\subsection{Effect of sampling length}
As one might expect, we found that longer sampling lengths (i.e. longer time series or more years sampled) from simulated populations represent the entire simulation better than smaller samples do. The quality of the representation, as measured by a comparison of sample variance to variance in the entire simulation, appears to approach some parametric value asymptotically. In this study, that asymptotic approach requires samples measured in units of centuries (Fig.\ \ref{vt}). Apparently, random walks to exceptionally large or small population sizes, caused only by demographic stochasticity, occur on time scales on the order of hundreds of years. This observation further supports our prediction that environmental stochasticity generates much of the variance observed in the 17-year dataset we study here. It also calls into question the generality of conclusions about demographic stochasticity drawn from samples even decades long. If such forces strongly influence population dynamics of species with similar life histories, modeling will be required to correctly characterize the dynamics; studies relying solely on statistical assessments of data at hand are likely to miss significant dynamical processes.
\subsection{Environmental stochasticity}
Our primary method for comparing simulation output to data focuses on matching means and variances. Within the portion of parameter space that admitted reasonable fits to the data, our simulations consistently matched the mean population size in the Bimini nursery (Fig.\ \ref{biminisharks}); however, simulations regularly generated variances distinctly lower than that seen in the actual data set, even when data variance fell with the middle 95\% of the distribution of simulation variance. Specifically, simulations on average account for approximately 33-49\% of the data variance. The variance in our models is generated entirely by demographic stochasticity and any instabilities caused by delayed breeding and age structure. What, then, caused the missing variance? We postulate it may have been a combination of variations in prey abundance, environmental stochasticity, including weather patterns and global climate change, habitat loss, and effects of fishing.
Whatever this environmental stochasticity is, we predict that its effects are relatively sparse, even though we underestimate the actual variance by a considerable amount. This prediction follows from a close inspection of the data from North Sound (Fig.\ \ref{biminisharks}). It appears that the high and low points in this time series (2005 and 2008, respectively) may be outliers. Removing these points from the data set reduces $s^2$ to 164.12 but leaves the mean at 75.13, both of which are in almost exact agreement with simulation mean and variance at default parameter settings. Therefore, the discrepancy between our simulations and the data appear to be driven by only two events in the 17-year data run.
The next step with this model is to incorporate environmental stochasticity so that more accurate assessments of population viability can be made \citep{OvaskainenMeerson2010}. However, such modifications will be difficult because stochastic environmental effects include an enormous array of possibilities. In addition, predator-prey dynamics, including juveniles as prey for both conspecifics and other species, should be modeled. Cannibalism, which has been documented in this species \citep{Morrissey1993,Guttridge2012}, needs careful attention because it can have very drastic effects on population dynamics \citep{Dennis2001,Ziemba2000}.
\subsection{Implications of model}
Our modeling approach is applicable to other shark populations as well as other megavertebrates. The age structure and stochastic birth and death rates can be easily altered to fit a different population. This type of model is also useful when applied to species that exhibit distinct stages in their lifecycle that can then be ordered into specific physiological classes \citep{Crouse1987,Caswell2001}. The lemon shark life cycle is naturally broken into three primary stages, but this structuring can be relaxed to fit life cycles in many other species. In addition, this model is ideally suited to study populations where basic data on annual population size is available; all the technique requires are estimates from the data of mean population size and interannual variance from at least some age class or classes. However, even this requirement can be relaxed. Any measurement made in the field that can be mapped to a variable in the model could be used to determine which combinations of parameters are a good fit \citep{Hartig2011}. For example, if the only field data available is an estimate of subadult mortality, researchers can match distributions of that variable to the field data and thereby constrain the biologically relevant parameter space using the same techniques we employ here. The ability of this technique to predict bounds for parameters that are not easily estimated in the field has important implications for management and conservation, generating as it does predictions about which parameters and life cycle stages may be most sensitive to anthropogenic impacts such as overfishing or bycatch. For lemon sharks at Bimini, we show it is essential to include density-dependent mortality in the first age class and to incorporate delayed breeding to predict even basic population dynamics. We also show that adult lemon sharks must have a mortality rate below 0.17 in order for the population to remain viable. Although we have a relatively long data set (17 consecutive years), longer time series may be required to capture important, rare stochastic events. These types of events, whether they be environmental or demographic, seem to be the primary factor in driving the fluctuations in the population size of juvenile sharks in Bimini.
\section{Acknowledgements}
We would like to thank all the volunteers and staff who have helped collect data at the Bimini Biological Field Station. We are also thankful to M. Braynen and C. Higgs, Directors of the Bahamas Department of Fisheries, for issuing a scientific permit in support of the field research. We would like to thank Jos\'{e} Anad\'{o}n and Jesse Senko for their helpful comments in preparing this manuscript.
%\bibliography{/Users/easton2/Desktop/White_bib}{}
\bibliography{White_bib}{}
%\bibliographystyle{plain}
\bibliographystyle{plainnat}
\newpage
\section{Tables}
\begin{table}
\begin{minipage}[c][\textheight]{\textwidth}
Table 1. Notation and interpretations of model parameters, their default values, ranges and sources for the lemon shark (\emph{Negaprion brevirostris}). \\[1ex]
\begin{tabular}{c p{4.2cm} c l p{4.8cm}}\hline
Parameter & Meaning & Default & Range & Source \\ \hline \hline
$ \lambda$ & Pups born per female & 6.1 & 1-18 & \citep{Feldheim2002, Feldheim2004}\\
$h$ & Juvenile mortality Hill parameter & 1 &NA & This paper \\
$k$ & Juvenile mortality shape parameter & 100 & 0-200 & This paper \\
$t_{max}$ & Maximum age for adult & 25 & 20-35 & \citep{Cortes2002a,Hoenig1990}\\
$x_m$ & Age at maturity & 12 & NA & \citep{Cortes2002a,Hoenig1990}\\
$\mu$ & Mortality rate for all animals above age one & 0.15 & 0.05-0.30 & This paper \\ \hline
\end{tabular}
\end{minipage}
\end{table}
\begin{table}
\begin{minipage}[c][\textheight]{\textwidth}
Table 2. Indirect methods used to calculate mortality rates. Here $M$ and $Z$ represent natural and total mortality, respectively. Similar analysis as \citep{Heupel2002} and \citep{Knip2012}. \\[1ex]
\begin{tabular}{l p{7.8cm} c }\hline
Method & Relationship & Value \\ \hline \hline
Hoenig (1983) (fish) & $ln(Z)= 1.46 - 1.01 \: ln(t_{max})$ & 0.167 \\
Hoenig (1983) (cetacean) & $ln(Z)= 0.941 - 0.873 \: ln(t_{max})$ & 0.154 \\
Hoenig (1983) (combined) & $ln(Z)= 1.44 - 0.982 \: ln(t_{max})$ & 0.179 \\
Pauly (1980) & $log(M) = -0.0066 - 0.279 \: log(L_\infty) + 0.6543 \: log(K) + 0.4634 \: log(T)$ & 0.140 \\
Jensen (1996) (age) & $M= 1.65/x_m$ & 0.138 \\
Jensen (1996) (growth) & $M= 1.5 \: K$ & 0.086 \\
Jensen (1996) (Pauly) & $M= 1.6 \: K$ & 0.091 \\
\hline
\end{tabular}
\textbf{Note:} Life history parameters are based on \citep{Brown1988}. $K$, body growth parameter (0.057); $L_\infty$, maximum theoretical length (317.65 cm); $x_m$, age at maturity (12 years); $t_{max}$, maximum age (25); T, mean temperature ($27.1\,^{\circ}\mathrm{C}$, \citep{Newman2007})
\end{minipage}
\end{table}
%
\newpage
\section{Figures}
\begin{figure}[htbp]%function to make graph float
\centering
\includegraphics[width=8.64cm,height=5.76cm] {biminisharks4.pdf}% requires the graphicx package
\caption{Juvenile population data from the past 17 censuses in the North Sound of Bimini.}
\label{biminisharks}
\end{figure}
\begin{figure}[htbp]%function to make graph float
\centering
\includegraphics[width=7cm,height=7cm] {littersize(after1995andNS)withdiscretecurve.pdf}
\caption{Distribution of litter size per female lemon shark in North Bimini. Grey bars: data from \citep{Feldheim2002, Feldheim2004}, from 1996 to 2010 (n=264). Red curve: discrete Poisson distribution, $\Pr\{N=i\}= {e^{-\lambda}} {\frac{\lambda^i}{i!}}$, with $\lambda$ equal to the mean of the litter size distribution depicted by the grey bars.}
\label{littersize}
\end{figure}
\begin{figure}[htbp]%function to make graph float
\centering
\includegraphics[width=10cm,height=10cm] {9000paramcombo(100trials_300years).pdf}% requires the graphicx package
\caption{Region of parameter space in which simulations exhibited a ``good" fit to the data of the lemon shark population based on criteria described in the main text. Each filled circle represents one of the 9000 parameter combinations that met the criteria of a good representation. The change in color represents degree of half saturation value, with red indicating smaller values of $k$.}
\label{alldata}
\end{figure}
\begin{figure}[htbp]%function to make graph float
\centering
\includegraphics[width=12cm,height=7cm] {goodrep_actual_vs_poisson(0-400).pdf}% requires the graphicx package
\caption{Left: Half-saturation value versus the mortality rate for subadults and adults for series of combinations utilizing the actual distribution of litter sizes for fecundity rate. Right: Same as left but uses a Poisson distribution for fecundity rates. Both pictures represent cases when $\lambda$ was set at 6.1 for the Poisson distribution which is equivalent to the average of the actual distribution of litter sizes. }
\label{halfsatvsmort}
\end{figure}
\begin{figure}[htbp]%function to make graph float
\centering
\includegraphics[width=10.1cm,height=7cm] {var_vs_samplinglength2.pdf}
\caption{Simulation variance as a function of sample size. The sample at Bimini is a total of 17 years (indicated by the vertical red line). The green line represents the variance in the actual population size ($s^2$=498). }
\label{vt}
\end{figure}
\end{spacing}
\end{document}