Using time series data to assess recent population dynamics of Bald Eagles in the northeast United States

We explored the use of time series data to study the recent demographics of Bald Eagles in the northeast United States. We merged a stage-structured population matrix model (PMM) into a combinatorial optimization algorithm (COA) to estimate stage-wise abundances, survival rates, and other dynamics of the eagle population from 1990 to 2018. The PMM includes biannual periods (breeding and nonbreeding), density dependence, and stage-wise migration and dispersal. We demonstrated theoretically that the information contained in 17 yr of time series data of the reproductive stage would have contained all the information for the entire population. We then tested the performance of the COA in estimating transient and asymptotic population dynamics using simulated data. Finally, we assessed the COA under differing modeling conditions (without stage-wise dispersal and in the absence of biological assumptions). We found that time series data of the reproductive class, the skeleton of the PMM, and COA are sole ingredients to reliably estimate stage-wise abundances and survival, but additional biological assumptions are helpful in honing precision of the estimates. The algorithm is a tool to assess population dynamics when demographic data are limited to counts of reproductive adults.


INTRODUCTION
Population matrix models (PMMs) are used to assess population viability, understand population dynamics, and provide guidance for management (Caswell 2001). From the PMM arise important demographic measures of shortand long-term viability (Caswell 2001). Given the information on population performance that may be extracted from the model, the PMM has been applied to thousands of biotic systems around the globe (Salguero-G omez et al. 2014;2016).
The structure of a PMM can be adapted to the characteristics of any life history (e.g., see Salguero-G omez et al. 2014Salguero-G omez et al. , 2016, but applied PMMs are frequently limited by reliable information on individual parameter values (Gallardo et al. 2019). Estimation of parameters is a key endeavor in ecology, and most of the current methods rely on extensive field effort. Limitations in funding, time, or personnel have stymied the adaptation of a PMM to countless species whose management may otherwise benefit from the tool.
We explored the use of a combinatorial optimization algorithm (COA) to attain population projections, estimates of stage-wise survival rates, matrix elements, migration, and dispersal, and functions of these parameters (eigenvalues, eigenvectors, sensitivities, elasticities, damping ratios, cumulative, stochastic, and transient growth rates, harmonic and arithmetic means, inertia, reactivity, maximum amplifications and attenuations, and life tables), all from annual counts of the reproductive stage.
The use of a COA in PMMs is an emerging computational technique in ecology. Any PMM can be rewritten in a series of difference equations, where combinations of the elements of the PMM are the coefficients. For example, let an arbitrary 3 9 3 matrix contain nine elements, where a ij is the element in the ith row and jth column. Let a 3-stage vector consist of juveniles, subadults, and adults at time t (J t ; S t ; and A t ). Then, the equation can be rewritten as the series of difference equations.
J tþ1 ¼ a 11 J t þ a 12 S t þ a 13 A t ; (2) S tþ1 ¼ a 21 J t þ a 22 S t þ a 23 A t ; and A tþ1 ¼ a 31 J t þ a 32 S t þ a 33 A t : Eqs. 2-4 reveal that the matrix elements are mixed with the stage-wise time series counts to project abundances at time t þ 1. Traditionally, the set of a ij s and vector of J t ; S t ; and A t are used to calculate J tþ1 ; S tþ1 ; and A tþ1 . But Eqs. 2-4 suggested that information on J t ; S t ; and A t and J tþ1 ; S tþ1 ; and A tþ1 may be used to calculate the set of a ij s.
Indeed, Hanley et al. (unpublished manuscript) tested this idea with the COA to estimate sets of parameter values for Northern Spotted Owl (hereafter "NSO"; Noon and Biles 1990) from adult time series data. The COA successfully recovered simulated population properties in a blind test given nothing but simulated adult time series data (Hanley et al., unpublished manuscript). This work suggested that observed changes in time step abundances contain information regarding the underlying matrix parameters, even if those parameters may not be directly observable in the field. Using a PMM in combination with a COA, modern computing power, and empirical sequential counts of breeding adults, important demographic parameters may be estimated with otherwise minimal (or unobtainable) field data.
We explored the use of adult time series data and the COA in estimating demographic parameters for Bald Eagles (Haliaeetus leucocephalus) in the northeast United States (NE) during period of population recovery. Like the NSO, the model for Bald Eagles contained three stages of long-lived individuals, but unlike the NSO, the population of eagles was assumed to contain birds that spent time away from their breeding area and to have fecund birds that failed to reproduce each year (Buehler et al. 1991a). The COA initially estimated parameters in a plummeting population of NSO, while the focus herein was on estimating parameters of an expanding population. Given these discrepancies in system-scale dynamics, it was unclear whether the COA would perform as well for eagles as it did for NSO.
We hypothesized that the adult time series data of breeding adults (hereafter "data") could be integrated into a symbolic, skeletal structure of a PMM and estimate underlying parameter values. Given the mathematics in Hanley et al. (unpublished manuscript), we saw no reason that the COA could not be adapted to satisfy the mathematics of eagles, but we noted that additional complexities in the structure of the PMM would likely mandate a longer sequence of data than that needed for NSO.
Our study objectives were fivefold: (1) to derive a PMM that reflected the biology of Bald Eagles, (2) to identify the length of data necessary to estimate all model parameters, (3) to assess the use of a COA for parameter estimation, (4) to assess the conditions under which the COA could produce reliable estimates, and (5) to assess the robustness of COA to unmet biological assumptions.

The population matrix model
The Lefkovitch matrix (Lefkovitch 1965) contained three stages: a juvenile stage (hatchling; hereafter "SI"), a combined immature and nonreproducing adult stage (SII), and a reproducing adult stage (SIII). We assumed an equal sex ratio (Millsap et al. 2014, Gallardo et al. 2019) and modeled females only, but the modeling of males could be conducted without issue.
The life history was categorized by biannual periods of roughly equal length (USFWS 2019a). The window of the breeding period varies by latitude (USFWS 2019a), but we assumed that it occurred between January and June (USFWS 2019a). The window can be shifted without issue. We assumed that the breeding period was approximately six months in length (35 d incubation, 82 d in nest, 70 d until independence;Buehler 2000) and that SIs fledged during the following non-breeding period. We further assumed that SII lasted a minimum of four years (Hunt et al. 2009, USFWS 2019a and that SIIIs remained in SIII (and reproduced) or regressed to SII (and failed to reproduce; Hansen and Hodges 1985) each year. Equations with differing spans of SII are in Appendix S1: Part X.
Migration and dispersal were modeled during the non-breeding periods, but location during the breeding period determined residency (Hansen and Hodges 1985). Birds that were present in the NE all year were year-round residents (hereafter "YR"), while birds that migrated during the non-breeding period were part-year residents (PY). We treated the data as counts of mating pairs that only bred once per year, with possible replacement clutches only if first nesting attempts failed (Wood andCollopy 1993, Buehler 2000). We had no way to verify this assumption given our data (see Appendix S1: Part XI).
These life history characteristics were the mathematical building blocks of the PMM (de Kroon et al. 2000). The PMM for the breeding period was 0 0 a 13 0 a 22 a 23 0 a 32 a 33 2 4 3 5 ; (5) where a 13 was the fertility of SIII, a 22 and a 33 were probabilities that SIIs and SIIIs would survive but fail to transition,a 32 was the probability that SIIs would survive and transition into SIIIs, and a 23 was the probably that SIIIs would survive and transition into SIIs. No transition from SI to SII occurred during the breeding period due to the parameterization of the PMM. Survival differed across the year. The PMM for YR during the non-breeding period was where a was the probability that SIs would survive and fledge, b and e were the probabilities that SIIs and SIIIs would survive but fail to transition,d was the probability that SIIs would survive and transition into SIIIs, and c was the probably that SIIIs would survive but regress to SIIs. A similar matrix was used to represent the PY during the non-breeding period (Table 1). Migration depends on age and location, but eagles tend to return to the same territory year after year (Buehler 2000, Hunt et al. 2009, USFWS 2019a. Departure is often more synchronous than arrival (Buehler et al. 1991b), but we modeled migrations as grouped events with no mortality and instantaneous transit. We included parameters for migration and dispersal (Table 1).
Matrix elements represented the compound probability of survival and transition (Doak et al. 1994). Stage-wise survival was constrained by unity, while the biology governed the probability of the transition in the ith stage (x i ; i ¼ 1; 2; 3Þ. SIs perished or fledged by fall, so x 1 was equal to one. Eagles reached sexual maturity at five years (Hunt et al. 2009), which we defined as the onset of SIII. Successful maturation into SIII required x 2 =1/6. Derivation of the density dependence parameter (x 3 ) mandated that we assume a life expectancy, so we used the published average life span (25 yr; Buehler 2000, Lutmerding and Love 2017, USFWS 2019a). The expected value of the death at the onset of the 26th year required where S iB was survival in the ith stage during the breeding period. We did not assess potential biases that may have arisen from the use of this life expectancy, but Appendix S1: Part I Eqs. S32-S61 and Part X Eqs. S112-S157 may be modified to do so. Letting t represent time (yr) and J t S t A t ½ 0 be the vector of abundances in SI (J t ), SII (S t ), and SIII (A t ) in the breeding period, then abundances in year t þ 1 were fully captured by the matrix skeleton: See Appendix S1: Part I Eqs. S1-S27 for the derivation of this complex PMM.

The difference equations
The matrix skeleton (Eq. 8) can be rewritten in an abbreviated form: where E, F, K, and I are defined as in Table 2.
Equations 10, 11 imply and By substitution, Eq. 11 became Continuing with additional substitutions, we arrived at the nth-order difference equation (17) where A tþ1 is the count of SIII at time t þ 1,S tÀn is the count of SII at time t À n, A t is the count of SIII at time t, and A tÀn is the count of SIII at time t À n. The derivation is in Appendix S1: Part I (Eqs. S62-S90).

The data
The data consisted of 29 consecutive years of annual counts of SIIIs, taken during the breeding period by agencies and societies in the NE (Smith and Clark 2018, New  The data were assumed to represent a census of the segment of eagles in the NE that made themselves available for observation. We did not have information regarding sampling designs nor the error of estimates in the counts from each state. We therefore limited our inference to that segment of the population that are observable and assumed that all data were accurate.

The combinatorial optimization algorithm
The goal of the COA was to estimate matrix elements from SIII time series data. The COA searched the parameter space at random for combinatorial values of E; F; K, and I (hereafter "sets") in the life history skeleton that could have generated identical (or similar) SIII trajectories (e.g., the A tþ1 ; A t ; A tÀ1 ) to reality. The sets were linked to PPM parameters via Table 2.
Parameter values were bounded by their traditional ranges in the PMM (Table 1). For example, stage-wise survival is limited to unity, and was therefore bounded between zero and one. The algorithm therefore checked candidate sets where stage-wise values satisfied this criterion. We further narrowed the algorithmic search to the region of the parameter space with (1) fledgling survival between 0.63 (long movements; Wood 1992) and 1.0 (short movements; Buehler et al. 1991a, Hunt et al. 2009), (2) fertility between 0 and 1.5 (Bent 1937, Buehler 2000, (3) the probability of reaching SIII (breeding age) between 0.35 and 0.65, (4) life expectancy >10 yr, and (5) initial conditions of SII being at most 3 9 SIII. Detailed explanation is in Appendix S1: Part II.
The COA searched for combinations of E; F; K, and I that could have produced the data. In all, it estimated five matrix parameters for the breeding period, 10 matrix parameters for the nonbreeding period, four migration and dispersal parameters, and abundances of all stages through time. From these parameters, the calculation of other population characteristics followed. Technical details of the COA are in Appendix S1: Part II.

The reliability of the COA
We assessed whether the COA successfully estimated underlying parameter ranges from simulated SIII trajectories. We prepared a mock PMM according to the complex structure in Eqs. 5-8, and then used a random number generator to initialize the exact values in a mock set of underlying parameters. We used the mock set to project a time series of SIII eagles across 29 yr. We furnished this mock SIII trajectory to the COA, ran the COA with a randomly generated Table 2. The notation for the coefficients in the difference equations (sets) and their equivalence in population matrix model parameter notation.

Difference equation
Parameter equivalence set of initial conditions, and then compared the COA estimates against the underlying mock parameter values (Validation 1, hereafter "V1").
Comparisons were conducted for annual, biannual, biannual stage-structured abundances, migration and dispersal parameters, survival rates, asymptotic quantities (eigenvalues, stable stage structure, reproductive values, sensitives, elasticities, and damping ratios), transient quantities (cumulative, transient, and stochastic growth rates, reactivities, inertia, and maximum amplification and attenuation), and genetic quantities (harmonic and arithmetic means, theoretical loss of genetic diversity). Further information regarding these population characteristics is found in Appendix S1: Part IX.

The robustness of the algorithm to unmet assumptions
We assessed the performance of the COA when parameters values were left unbounded. We repeated the validation procedure with randomly generated parameter sets, but systematically relaxed each of the eagle-specific biological assumptions (Table 3). Validation 2 (V2) removed the assumed bounds on fledgling survival, Validation 3 (V3) removed the assumed bounds on fertility, Validation 4 (V4) removed the assumed bounds on the probability of reaching breeding age, Validation 5 (V5) removed the assumed lower bound on life expectancy, Validation 6 (V6) removed the assumed initial proportion of SIIs, and Validation 7 (V7) omitted all five of these assumptions (also see Appendix S1: Part IV).
We validated the use of the COA with a simple PMM. The simple PPM assumed that stage-wise matrix parameters were identical across the year (e.g., the value of a in the breeding period was equal to the value of a* in the non-breeding period; Appendix S1: Part V, Tables S6-S7). We repeated the assessment of the COA with this simple PMM given all five biological constraints (Validation 8, V8), and with no constraints (Validation 9, V9; Table 3).
We validated the use of the COA in a closed system. We repeated the assessment of the COA when the complex PMM contained no dispersal and all five biological constraints (Validation 10, V10), and when the complex PMM contained no dispersal nor biological constraints (Validation 11, V11). We then repeated the assessment of the COA when the simple PMM contained no dispersal and all five biological constraints (Validation 12, V12), and again with no biological constraints (Validation 13, V13).
Finally, we validated the use of the COA when the age at first breeding was four (Validation 14, V14), six (Validation 15, V15), and three (Validation 16, V16). In V14-V16, all five biological assumptions were assumed to exist (Table 3). Derivations, software, and further details of V1-V14 are in Appendix S1: Parts III-X.
We prepared an interactive software application to display the results of the COA using real data. The EaglePOPd software contains the results of the COA for the NE, plus additional narratives regarding population dynamics and important references. EaglePOPd is available at https://cwhl.vet.cornell.edu/tools/eaglepopd.

RESULTS
The difference equation suggested that complete information on system-level population dynamics is contained in 17 yr of SIII data. When at least 17 yr of sequential data exist, the S tÀ17 term in Eq. 17 (S tÀ17 being the counts of SII at time t À 17) tends to zero, leaving where A tþ1 , A t , and A tÀ17 are the counts of SIIIs at times t þ 1, t, and t À 17, respectively, while K; E; F; and I are the coefficients that contain the matrix elements (see Table 2). Appendix S1: Part I shows the full derivation (Eqs. S62-S90 and Figs. S4-S6). The COA accurately estimated simulated parameter ranges for many demographic properties given a randomly selected example model (Tables 3, 4). In Table 4, the validation received a Y if the estimated range included the true value of the mock parameter from the simulated example. A full catalog of comparisons for all 78 demographic parameters for each of the 16 validations of the simulated model is in Appendix S1: Parts III-X, Figs. S7-S423.
The COA successfully predicted the ranges of stage-wise abundances, SI and SII matrix elements, stochastic growth rates, and the biannual survival rates for the simulated model under ❖ www.esajournals.org most modeling conditions, and under the relaxation of several biological assumptions (Table 4). The COA in combination with the complex PMM accurately estimated the ranges of annual and biannual abundances for all stages across all 58 biannual time periods. As well, the COA successfully estimated survival of all three life stages in both the breeding and non-breeding periods.
The COA was robust to violations in several biological assumptions but performed much poorer by comparison when fertility was unbounded (V3; Table 4). Simplification of the structure of the PMM (V8-9, V12-13) improved estimation of movement parameters (migration and dispersal), but decreased accuracy in the estimation of abundances. The omission of dispersal (V10-13) caused the COA to perform poorer by comparison in the estimation of annual and biannual abundances, but improved estimation of the stable stage proportions. Shifts in the age at first breeding to four or six (V14-15) did not significantly alter performance as compared to the age of first breeding of five (V1), but a shift in the age of first breeding to three caused poorer overall performance.
The COA exhibited marginal overall performance in predicting characteristics that were functions of matrix elements, such as eigenvalues, eigenvectors, and damping ratios. As well, the COA was marginal in recovering transient characteristics for the example model (Table 4).
We identified the algorithm specifications that produced the most accurate estimate for each parameter. We generated estimates of the resulting dynamics given the data in the NE. See EaglePOPd (https://cwhl.vet.cornell.edu/tools/ eaglepopd).

DISCUSSION
It is long known that using an equation with three variables, the unknown value of the third may be calculated given information on the other two. We applied this idea to investigate whether PMM elements may be estimated using information on abundances at times t and t þ 1 only. Further, sequential substitutions in the set of difference equations suggested that matrix elements may be estimated given the information on the reproductive stage only.
Time series data of the reproductive stage can be thought of condensed information regarding the sum of system-level signatures of the myriad ways individuals can replace their forebearers. For example, the COA in the NSO model required a time series of four years, because the full completion of generational replacement occurred in four years or less (Noon and Biles 1990). Here, we parameterized the life history of eagles using a complex matrix skeleton that allowed birds to reproduce (or not) across their average 25-yr life The Y indicates yes. PMM, population matrix model.
❖ www.esajournals.org span. Because individuals could enact countless permutations of transitions up to and including reproduction, more data were necessary to capture the total signatures of the permutations. The equations revealed that at least 17 consecutive annual counts of reproductive eagles constituted the minimum amount of data necessary to capture the myriad of life cycle permutations for Bald Eagles.
Comparison of the COA estimates with the underlying mock parameters revealed that the COA was capable of correctly estimating parameter ranges for many important demographic traits given nothing but time series abundances of reproductive adults. For example, in both the simple and complex PMM structures, the COA correctly identified abundances of SIs and SIIs in 58 consecutive biannual time periods from 29 yr of data on SIIIs (Table 4). Given the random selection of the experimental set, the probability that these successful estimations were driven by chance was Table 4. The results of the combinatorial optimization algorithm (COA) in estimating parameters under differing modeling conditions.
Notes: A Y indicates that the COA accurately captured the true value of the parameter in its output range. Model outputs are i, the annual counts; ii.a-ii.c, biannual counts for SI, SII, and SIII; iii.a-iii.b, incoming migration proportions for SII and SII; iv.a-iv.b, outgoing migration proportions for SII and SIII; v.a-v.b, incoming dispersal proportions for SII and SIII; vi.a-vi.b, outgoing dispersal proportions for SII and SIII; vii.a-vii.b, the spring survival for SII and SIII; viii.a-viii.b, the fall survival for YR SI, SII, and SIII; and ix.a-ix.c, the fall survival for PY SI, SII, and SIII. A full catalog of figures that compare algorithm ranges and the underlying mock parameters is in Appendix S1: Parts III-X, specifically Figs. S7-S423. PY, part-year. Ellipses signify not applicable, parameters were not tested in that validation. ð Þ 58 ð Þ= 2 58 À Á 2 58 À Á ¼ 4:04 Â 10 À32 . Similar exceptional performance existed in estimating the ranges of stage-wise survival (for both the breeding and non-breeding periods), most matrix elements (a 13 ; a 22 ; a 23 ; a; b; c; d; aÃ; bÃ; cÃ; dÃ), and for the stochastic growth rate.
The COA incorrectly estimated SIII stasis (a 33 ) in nearly all scenarios. Incorrect estimation of a 33 cascaded into incorrect estimates of the population properties that are derived from the matrix elements. For example, the asymptotic growth rate, sensitivities, elasticities, stable stage distributions, and reproductive values rely on accurate estimations of each matrix element, so incorrect values of a 33 propagated into incorrect estimates in each of these other quantities. Should the COA accurately estimate all matrix elements, improvement in the accuracy of each of these quantities would follow.
It is of immediate interest for the broader use of this tool to assess whether an alteration in the number of algorithmic runs, the probability of selecting a worse trajectory (in the optimization), or an alteration in the selection of the initial conditions would improve the performance of matrix estimation. We note, however, that this assessment was based on the performance of the COA in one scenario, and the COA may perform differentially given a differing scenario. We provide software that may be modified to tailor, assess, and hone the COA to specific applied scenarios.
While it is generally insufficient to assess the overall performance of a tool based on its performance in one example, we point out that the mock scenario consisted of a random set of coefficients that were selected by a random number generator. Table 4 therefore illustrates the performance of the COA in estimating population parameters when no attempt was made to control the nature of the simulation. In lieu of repeating these validations ad infinitum (to generate a probability distribution of the performance of the COA under all plausible scenarios), we instead provide software that may be modified to assess performance under specific scenarios.
The empirical time series data were collected by professional biologists in the NE from 1990-2018. We assumed that the data accurately reflected the abundances of SIIIs in the region during that time. However, habitats are not equally accessible for observation, nesting behavior could have been missed, environmental fluctuations exist in the size of natural populations (see, e.g., Saether et al. 2000, 2004, Coulson et al. 2001, Stenseth et al. 2003, Berryman and Lima 2006, or heterogeneities in sampling methodologies between biologists, agencies, or organizations may have introduced error in the empirical trajectories. We had no way to know if and to what extent the empirical time series deviated from reality in the NE, and we did not test how such derivations could have biased the results. Instead, we provided detailed documentation regarding our data, the derivations, the software, and each step in our methodology, and encourage its use for deriving equations to incorporate such sources of error into the COA (see Appendix S1: Parts I-XI).
Theory and simulation reveal that ranges of parameter values are attainable by the COA, but exact point estimates of each parameter are not. To see why, consider the example 8 = U9W, where U and W are natural numbers. It is fact that 8 = 492 = 891 = 294 = 198, so we can infer that U 2 (4, 8, 2, 1) and W would follow deterministically. Without additional information, we could say with certainty that 1 ≤ U ≤ 8, but the exact value of U would remain unknown. If, for example, W was known to be an odd number, then we could infer the exact value of U (8). Thus, the purpose of using published information to constrain values of Bald Eagle parameters was to pare candidate sets with the hope that the remaining combinations of parameter would carry higher information content for conservation.
Biological assumptions were therefore helpful in increasing information content, but they needed to be used cautiously. For example, in V1-complete with all five biological constraints-the COA estimated 400-2400 SIs and 0-6000 SIIs would exist in 2018 (see Appendix S1: Figs S7-S423). In V7, however, the COA estimated 200-2500 SIs and 0-17,000 SIIs. Both COAs successfully captured mock SI and SII trajectories for 2018, but the constraints in V1 made the estimates more precise by comparison with those in V7. Indeed, such constraints generally yielded a higher precision in estimates, but too many restrictive assumptions could result in inaccuracy (Table 4). Alternatively, too few restrictions often resulted in uninformative results (where the probability of survival was calculated to be between zero and one, for example).
Overall, the reliability and robustness of the COA in producing accurate and precise ❖ www.esajournals.org parameter estimates should be assessed on a case-by-case basis, preferably with the inclusion of only the best available information for the species. We point out that the values used herein to restrict the parameter space have undergone vigorous peer review prior to its use.
Despite several consequential technical decisions and the overall inability to generate familiar point estimates, even ranges of parameter values may be informative in advancing conservation. First, the COA seized emerging techniques in computational biology to analyze data that have long been collected but may otherwise be unanalyzed. Second, the COA aids in understanding dynamics of population processes that are currently not well understood using current field techniques. Third, the COA identifies the plausible limits of life history plasticity contained in a life cycle. Fourth, the COA allows the study of historical conditions, from which researchers may compare changes population dynamics over time. And finally, this technique allows interdisciplinary teams to compare population-scale dynamics of subpopulations in the presence and absence of veterinary diseases.
We used both theory and experimentation to show that the use of adult time series data in conjunction with a PMM and a COA constitutes an alternative emerging method to estimate viability properties for Bald Eagles in the NE (and beyond).