MODELING THE DAILY ACTIVITIES OF BREEDING COLONIAL SEABIRDS: DYNAMIC OCCUPANCY PATTERNS IN MULTIPLE HABITAT PATCHES

We constructed differential equation models for the diurnal abundance and distribution of breeding glaucous-winged gulls (Larus glaucescens) as they moved among nesting and non-nesting habitat patches. We used time scale techniques to reduce the differential equations to algebraic equations and connected the models to field data. The models explained the data as a function of abiotic environmental variables with R(2 )=0.57. A primary goal of this study is to demonstrate the utility of a methodology that can be used by ecologists and wildlife managers to understand and predict daily activity patterns in breeding seabirds.

1. Introduction.Ecologists study the distributions, abundances, behaviors, and interactions of living organisms.Identification of environmental determinants that elicit temporal patterns in these systems is important to both community ecologists and resource managers.Dynamical systems theory offers powerful but underutilized tools for identifying temporal patterns in data and for identifying environmental determinants associated with those patterns [11].
Differential equations have been used to model the diurnal movements and behaviors of harbor seals (Phoca vitulina) [6], herring and great black-backed gulls (L.argentatus and L. marinus) [7], and glaucous-winged gulls (Larus glaucescens) [3], [8], [9], [10], as well as the consequences of cowbird parasitism on avian communities [5] and the dynamics of waterfowl [16].A number of these studies show that a variety of environmental cues play crucial roles in determining diurnal movements and behaviors.
In a study of glaucous-winged gulls during chick-rearing season, Damania et al. [3] used a system of differential equations to model the diurnal movement among habitat patches at Protection Island National Wildlife Refuge, Washington.A standard compartmental modeling approach was used to track the occupancy of three focal habitats adjacent to a large breeding colony.The model compartments consisted of a jetty used for loafing; a small marina used for bathing and drinking; a beach used for loafing, bathing, and drinking; and an "other" compartment consisting of all other locations, including the nesting colony itself.The nesting colony, however, comprises a crucial component of the seabird ecosystem during the breeding season.In this study we substantially revise the model of Damania et al. to include nesting colony dynamics.The resulting differential equation model tracks glaucous-winged gull occupancies in the colony, jetty, marina, and"other" locations.Unlike the colony, which is occupied exclusively by nesting gulls, noncolony habitats are occupied by both nesting and non-nesting gulls and census counts cannot differentiate between these two categories of birds.Hence, the model must keep track not only of the numbers of nesting gulls in each habitat, but also of the numbers of non-nesting gulls in the jetty, marina, and "other" locations.This necessitates a seven-compartment model.Furthermore, the model must account for the fact that at least one of each nesting pair attends its territory at all times during chick-rearing season, as well as the fact that non-nesting gulls are promptly driven from the colony.
A primary goal of this study is to demonstrate the utility of a methodology that can be used by ecologists and wildlife managers to understand and predict daily activity patterns in breeding seabirds.Realistic models of dynamic field systems are difficult to connect to data because of computational constraints and limited data sets.We suggest that our methodology may be useful for solving some of these problems.55'W), Jefferson County, Washington, which contains the largest marine bird breeding colony in the state.Violet Point, a gravel spit on the southeast end of the island, supported about 2900 pairs of nesting glaucous-winged gulls during the data collection period.Three well-defined habitats were designated for study: a sample area within the colony, a small marina in the middle of the larger colony, and a rock jetty adjacent to the larger colony (Fig. 1).All areas, including the marina, are closed to the public and experienced little human disturbance.The gulls nest in short and tall grass habitats on the colony, engage in loafing behaviors (sleep, preen and rest) on the jetty, and use the marina for drinking, bathing, preening, and floating.

Data collection.
Hourly counts of the colony sample area, marina, and jetty were taken during daylight hours, 05:00-20:00 Pacific Standard Time (PST), 28 June to 16 July 2004 during chick-rearing season, from an observation point at least 100 m from the closest focal habitat.Observations were made using a 20-60x spotting scope from the observation point atop a 33m high bluff that borders the west end of Violet Point (Fig. 1); the presence of observers did not appear to affect the system in any way.A comprehensive census was taken for the marina and jetty.Each census of the colony sample area was multiplied by a scaling factor to represent a census of the entire colony.The scaling factor (41.93) was computed in the following way. 1) The maximum occupancy of the entire colony was calculated by multiplying the 2004 colony-wide nest count by two (2914 x 2 = 5828).This relies on the fact that each nest is associated with two birds and that non-nesting birds are promptly driven from the colony.2) The scaling factor for the colony was computed as the ratio of the maximum occupancy of the entire colony to the maximum occupancy observed within the colony sample area during the data collection period (139); that is, (5828 / 139 ≈ 41.93).This assumes that the dynamic occupancy of the sample area remained approximately proportional to the dynamic occupancy of the entire colony.
Hourly tide heights and solar elevations were obtained from the National Oceanic and Atmospheric Administration (NOAA).A weather station located 2 m above site elevation on the northwest end of Violet Point tracked a large number of environmental conditions on the colony, including humidity, wind speed and direction, temperature, barometric pressure, heat index, THSW index, rainfall, and solar radiation.The heat index and THSW index are computed from temperature and relative humidity as measures of how hot the air feels; THSW also includes the effects of solar radiation and wind speed on perceived temperature [17].
2.3.Model assumptions.Mathematical models are based on parsimonious sets of simplifying assumptions that capture the main dynamics of the system.Our models were based on eight assumptions: A1) The total number of gulls K that used the focal habitats was assumed constant over the data collection period (28 June to 16 July 2004).Although K varies during the year from essentially zero during the winter to many thousands at the height of the breeding season, the value is relatively stable during chick-rearing season.The value of K was estimated in the following way: During the 2004 data collection period, occupancy data were collected for two additional habitat patches-a sample beach area, and the water north and south of the spit extending approximating 200 m out from the beach.We assumed that all the gulls returned from remote feeding locations back to the colony, jetty, marina, beach, and water by dusk.Thus, K was calculated as maximum of the summed occupancies for these five habitats (7556, occurring at 1700 hours on day 193), where the colony sample occupancy was scaled up by a factor of 41.93 as explained under "Data Collection" above and the beach sample occupancy was scaled up by a factor of 5.5.The beach scaling factor was estimated from the total length of the beach and the lengths of sample areas of various densities using the method detailed in [14].
A2) The total number of gulls K in the system consisted of gulls without nests in the colony ("non-nesting gulls") and gulls with nests in the colony ("nesting gulls").That is, K = K 1 + K 2 , where K 1 and K 2 denote the number of non-nesting and nesting gulls, respectively.From the 2004 colony-wide nest count mentioned under "Data Collection" above, K 2 = 5828.From A1, A3) Colony occupancy at any given time was between K 2 /2 and K 2 .This assumption was based on the facts that at least one of a pair of nesting gulls remains at the nest at any given time during chick-rearing season and that non-nesting gulls are immediately chased away by colony residents.Since territoriality is not a predominant behavior in the other two study habitats (jetty and marina), no such maximum or minimum occupancy thresholds were assumed for those.
A4) The per-capita transition rates r ij from habitat j to habitat i were assumed proportional to powers of nine abiotic nondimensionalized environmental variables: time of day (hour), tide height (tide), solar elevation (sun), humidity (hum), wind speed on the colony (wind), temperature (temp), barometric pressure (bar), heat index (heat), and THSW index (thsw), where each environmental variable x was nondimensionalized and scaled so that 1 ≤ x ≤ 2. That is, where A5) The per-capita transition rates for nesting and non-nesting gulls into and away from the jetty, marina, and "other" were assumed to be the same.
A6) The system was assumed to recover rapidly after disturbance.That is, the environmental variables were considered essentially constant during system recovery.This is based on our observations over many years at the Protection Island colony that occupancies (and behaviors) essentially recover within 15 minutes after shortterm "point disturbances" by eagles and humans [3], [10], [11].
A7) The main source of noise in the census data was assumed to be demographic rather than environmental stochasticity.This assumption was motivated by the fact that major environmental correlates were incorporated explicitly into the model.Demographic stochasticity in this context arises from independent random binary choices of individual gulls as they arrive in or depart from a habitat.This is analogous to a stochastic birth-death process at the population level [8].
A8) Hourly residual model errors were assumed to be uncorrelated in time.That is, it was assumed that a stochastic event affecting the census at one hour would not affect the census an hour later, because of rapid recovery of the system postperturbation, assumption A6.Furthermore, at any given time t the covariances of the residuals between habitats were assumed small relative to the variances.That is, it was assumed that stochastic events mainly affected single habitats.We can describe the dynamics of the colony by the differential equation Here r ij is the per-capita rate at which gulls move from habitat j to habitat i, and hence the positive terms r 12 x 2 , r 15 x 5 , and r 17 x 7 are the inflow rates into the colony from the jetty, marina, and "other" location, respectively.By assumption A3, the factor (x 1 − K 2 /2) is the number of birds that are eligible to leave the colony.Thus, the negative terms r 21 ( ) and are the outflow rates from the colony to the jetty, marina, and "other" location, respectively.By assumption A2, the number of nesting gulls in the "other" location is x 7 = K 2 − x 1 − x 2 − x 5 .Thus, equation ( 2) can be written Given that the number of non-nesting gulls in the "other" location is by assumption A2, one can construct an equation for each of the remaining focal compartments.The model is the system of differential equations in which r 25 = r 34 , r 43 = r 52 , r 63 = r 72 , r 27 = r 36 , r 64 = r 75 , and r 46 = r 57 by assumption A5.Because of assumption A6, one can apply the multiple time scale analysis methods explained in detail in [11] to obtain an approximation to the dynamic steady state of model (4).The analysis involves the appearance of a small parameter multiplying each derivative on the "slow" time scale.When the left hand sides of equations ( 4) thus are assumed small (essentially zero), the values of x 1 , x 2 , • • • , x 5 can be approximated by solving the resulting nonhomogeneous linear system of algebraic equations using Cramer's Rule: Here A is the 5×5 matrix of coefficients on x 1 , x 2 , • • • , x 5 from the right hand sides of the differential equations in model ( 4), and each matrix A i is obtained from matrix A by replacing the i th column of matrix A with the vector of nonhomogeneous terms 2.5.Stochastic model.Demographic noise is approximately additive on the square root scale.That is, the variance of the residual model errors (the deviations of model predictions from observations) is stabilized by a square root transformation [4], [6], [8], [15].Given this fact and the assumption of demographic noise (A7), we can therefore model the square-root transformed residuals as where t denotes the discrete hours of data collection, the ε i (t) are drawn from a standard normal random distribution, and the σ i are positive constants.Thus, a stochastic model for the system is where the quantity inside the square is taken to be zero if negative.The ε i (t) are uncorrelated from hour to hour by assumption A8.Since the same parameters α ij appear in each of the x i (t), the random variables ε 1 (t), ε 2 (t), • • • , ε 5 (t) would in general be expected to covary.In this study, however, the covariances were assumed small relative to the variances and were taken to be zero, assumption A8.

Model selection.
Given assumption A5, model (4) has 12 per-capita flow rates r ij to be determined from data.By assumption A4, each of the 12 per-capita flow rates involves one coefficient parameter and nine exponent parameters, for a total of 120 model parameters.Numerical estimation of this many parameters with maximum-likelihood (ML) methods is computationally difficult and would require a large data set.To circumvent this problem, we designed a four-step procedure: 1) We first supposed that each of the exponents 1) had value -1, 0, or 1.This had two consequences.First, it reduced the number of ML parameters to the 12 coefficients, yet allowed each environmental factor to have a negative, zero, or positive effect on the flow rate.Second, it created a large number of alternative models in the form of model ( 4).Thus, the heavy computational burden was shifted from parameterization to model selection.
2) Next we eliminated most of the alternative models.Exhaustive model selection would have involved estimating ML parameters for each of the alternative models and choosing the model(s) that best explained the data.This approach is neither feasible nor desirable for such a large number of models; most such alternatives should be eliminated by other means before applying model selection techniques [1].We eliminated most of the alternative models in two ways.First, we used the biologists' knowledge of the system and statistical investigations of correlations between habitat occupancies and environmental variables to discard some of the possibilities.Second, we followed the "step-wise" method introduced by Damania et al. [3].In this approach, one pursues the (much easier) task of determining the best flow rates for each possible two-compartment model (in this case, colony and "other", jetty and "other", marina and "other"), then uses this information to determine the best flow rates for each possible three-compartment model; in this way one finally identifies a relatively small subset of alternative models for the larger system.
3) We then selected the "best" of the remaining alternative models to be the one with the smallest (fitted) sum of squared residuals after ML parameterization (see next section), and we discarded the others.Because each of the alternative models had the same number of ML parameters, it was not necessary to use an informationtheoretic model selection index such as the Akaike Information Criterion (AIC) [1].Note that the selected model had the property that each of the exponents on the environmental variables had value -1, 0, or 1.
4) Finally, we adjusted the exponents on the environmental variables in the selected model.In particular, exponents with value -1 (or 1) were decreased (or increased) by integer units until the "best" integer exponents were obtained.This was done by changing the value of the exponents by hand in the computer program and recalculating the sum of squared residuals until the minimal value was obtained.2.7.Model parameterization.ML parameters were estimated from hourly census data.The census counts of the jetty and marina did not differentiate between nesting and non-nesting birds.Thus, hourly model predictions for nesting and nonnesting birds were summed to produce total jetty and marina hourly predictions (x 2 + x 3 and x 4 + x 5 ) that were then compared to the hourly observations of jetty and marina occupancies.The residual error at time t for the i th habitat, computed on the square root scale, is given by Because the residuals were assumed uncorrelated across habitats (assumption A8) the ML parameter estimates are equivalent to the nonlinear least squares parameter estimates, obtained by minimizing the sum of the squared residuals for the three censused habitats (colony, jetty, and marina) as a function of the vector θ of parameters, where q is the number of data points [3].The minimizer θ is the vector of ML parameter estimates for the model.We produced the model predictions using MATLAB and minimized RSS(θ) using the Nelder-Mead Algorithm, a convenient downhill method for numerically finding minima of functions [13].
2.8.Goodness-of-fit.The goodness-of-fit was measured by a generalized R 2 where "mean i " denotes the sample mean of the square roots of the observations for the i th habitat.The R 2 estimates the proportion of observed variability explained by the model and thus provides a measure of the accuracy of the model's predictions.The higher the R 2 value, the better the model fit, with R 2 = 1 indicating a perfect fit.
3. Results and discussion.The selected model is shown in Table 1.Predictions and observations are shown in Figure 3A-D.The R 2 was highest for the jetty and lowest for the marina.The overall R 2 value was approximately 0.57.Colony counts and predictions were high in early morning, lower at midday, and highest in late evening (Fig. 3A).The model flow rates in Table 1 suggest that nesting gulls leave the colony when the sun is high and the tide is low.Low tides present the most opportune times for shoreline and mudflat feeding, and during the breeding season, low tides are most likely to occur between mid-morning and midafternoon.When conditions were reversed, patterns of flow were reversed.Daily low counts, but not predictions, decreased over the course of the data collection period; this was likely due to typical accumulation of nest failures, a source of midday attrition not included in our model.Table 1.Parameters and environmental variables for per capita transition rates of glaucous-winged gulls between colony (C), jetty (J), marina (M), and 'other' (O) on Protection Island.Colony R 2 = 0.47; jetty R 2 = 0.81; marina R 2 = 0.44; overall R 2 = 0.57.Daily jetty counts and predictions followed an occupancy pattern similar to that of the colony.Midday lows on the jetty, however, were essentially zero (Fig. 3B), reflecting the absence of territory defense in this habitat.Both nesting and nonnesting gulls use the jetty as a loafing site for preening, resting, and sleeping.Table 1 suggests that birds move to the jetty when the tide is high and the sun is low, conditions most likely to occur toward evening when feeding sites are not exposed.Nesting birds move from the jetty to the colony late in the day, suggesting that they use the jetty as a resting place between feeding sites and territorial activities on the colony.
Marina counts were relatively low until the last five days of data collection, when counts peaked at 141 at 12:00 hr on Day 193 (Fig. 3C).Marina counts were lowest in the morning, higher in the afternoon and early evening, and low again in late evening.On hot days, nesting gulls moved from the colony to the marina water (Table 1).They also moved here from the jetty on hot days and from "other" locations when the hour was late.Gulls unload significant amounts of heat from their feet [12]; thus, simply floating in the cool marina water serves an important thermoregulatory function.Gulls also use the calm marina waters for drinking and bathing.
"Observations" and predictions for the "other" habitat were generated by subtracting the sum of the observed/predicted occupancies for the three focal habitats from K, the total number of gulls in the system (Fig. 3D)."Observed" and predicted occupancies for the "other" habitat category decreased from mid-morning to evening.
Several comments and caveats are in order.First, it is important to note that the connection between "driving" environmental factors in the model and the behavioral dynamics is that of mathematical implication rather than scientific causation.Environmental "determinants" identified in this study are correlative and may or may not be causative [11].The identification of environmental determinants, however, narrows the search for cues that elicit behavior.
Second, although time scale analysis greatly simplifies parameterization by reducing computation time, and renders the methodology accessible to ecologists and wildlife managers, the resulting steady state models generally present an inverse problem that is not well posed [11].That is, parameters cannot in general be determined uniquely because various combinations of inflow and outflow rates can give rise to the same steady-state dynamics [6], [11].In particular, a given set of parameter estimates may not support the assumption of rapid transients in the original differential equations.Although this is immaterial to the usefulness of the steady state model as a tool to predict occupancies, it is important to remember that the parameter estimates themselves have no biological meaning except relative to each other in the context of the steady-state model.They should not be used in the differential equation models and cannot explain transient dynamics.In order to model transient dynamics, one must collect data densely over short time spans following disturbance of the system.A detailed discussion of this is given in [11].Furthermore, the environmental variables in the flow rates cannot in general be determined uniquely, for similar reasons.Extensive experience using the method suggests, however, that the modeling procedure robustly identifies environmental determinants.Most alternative models do not parameterize (one or more parameters approaches zero or infinity); the ones that do typically have flow rates based on similar combinations of environmental variables.Third, our procedure for identifying the 120 model parameters was a nonstandard combination of ML estimation, model selection techniques, and trial-and-error computations.The final model chosen (Table 1) does not necessarily represent the 120-parameter maximal likelihood incarnation of equation ( 4).Realistic models of field systems are difficult to connect to data because of computational constraints and limited data sets.We suggest that our procedure, although somewhat ad hoc, may be useful for solving some of these problems.
Fourth, an overall R 2 of 0.57, although quite high for dynamic field data from a complicated system, indicates a substantial amount of stochasticity relative to the deterministic model.Departures of model predictions from data are due to at least four sources of variability: 1) Observational error.Count inaccuracies occur as a result of poor light conditions, fog, or concealing vegetation; uncertainty as to whether birds near habitat boundaries are in or out of the habitat; and inflow and outflow during counts.2) Environmental stochasticity.Strong winds over the Strait and rain suppress gull movement; bald eagle flyovers and predation events disrupt "normal" activities of gulls; and human disturbances occur such as boats entering the marina.3) Demographic stochasticity.Gulls, like many vertebrates, exhibit considerable individual variation in behavior.This variation is due to innate differences in "temperament", different histories, and varying physiological states.4) Error in modeling assumptions.Each of the modeling assumptions is to some degree overly simplistic, leading to departures of model predictions from the data.For example, assumption A8 is not completely valid; some prolonged disturbances, such as deer walking through the colony, can last longer than one hour.Indeed, post hoc diagnostic analyses of residuals do show some correlation between temporally subsequent pairs of residuals in all three focal habitats (ρ ≤ 0.45).
Henson et al. [9] provide more detail on the caveats associated with this type of analysis, including the observation that several of the environmental variables are correlated (such as temperature and relative humidity) or related through calculated formulas (such as temperature and THSW).4. Summary.Gulls exhibit considerable between-individual variation in decision making.The choice to enter or leave a habitat by one gull may be influenced by behaviors of other gulls [2] or may be made independently.Given the complexity of colony life, individual decision-making seems particularly pronounced in nesting gulls compared to loafing and swimming gulls.Nesting gulls engage in a wide variety of activities including courtship, copulation, nest-building, egg-laying, incubation, provisioning young, thermoregulation, territory defense, and response to predators.Some of these activities involve entering and leaving the colony in response to a wide variety of contingencies such as approach by predators, the drive to access water to drink and thermoregulate during hot periods, and the exploitation of ephemeral food sources near the colony.Moreover, because nesting gulls enter all habitats but non-nesting gulls enter the colony relatively infrequently, a more complex system of equations was needed to model the system than that used by Damania et al. [3].
One goal of this study was to identify the environmental factors that drive movement of glaucous-winged gulls among three focal habitats in or adjacent to a breeding colony in the Strait of Juan de Fuca, Washington.The results of this study are consistent with those of earlier studies of Protection Island gulls that identify environmental conditions correlated with behaviors and habitat occupancies.Tide, solar elevation, temperature, and humidity are to a large degree responsible for shaping the daily movement patterns and behaviors of these birds [3], [8], [9], [10], [14].Another goal was to demonstrate the utility of a methodology that can be used by ecologists and wildlife managers to understand and predict daily activity patterns in breeding seabirds.This study suggests that mathematical modeling can play an important role in the establishment of priorities, goals, and policies in wildlife management.

2. 4 .
Deterministic model.As shown in Figure2, let x 1 = number of gulls in the colony x 2 = number of nesting gulls on the jetty x 3 = number of non-nesting gulls on the jetty x 4 = number of non-nesting gulls in the marina x 5 = number of nesting gulls in the marina x 6 = number of non-nesting gulls in the "other" location x 7 = number of nesting gulls in the "other" location

Figure 2 .
Figure 2. System of habitat patches, with separate compartments for nesting (N) and non-nesting (NN) glaucous-winged gulls on Protection Island. 181

Figure 3 .
Figure 3. Model predictions (continuous curve) and hourly observations (open circles) of numbers of glaucous-winged gulls in each habitat.Each graph shows one day of data, with the day of the year indicated in the upper left-hand corner.A. Colony, B. Jetty, C. Marina, and D.Other.The data points in B and C designated by a star exceeded the maximum value on the vertical axis.'Observations'/predictions for 'Other' were generated by subtracting the sum of the observed/predicted occupancies for the three focal habitats from K (the total number of gulls in the system).