Population Limits with Multiple Habitats and Mechanisms: I. Age-Independent Periodic Breeders

1. Introduction

Fretwell (1972) has discussed many of the factors influencing the selection and dispersal among various habitats for populations in seasonal environments. A number of studies of neotropical migrant birds, such as reviewed by Rappole et al. (1983), Rappole (1995) and Sherry and Holmes (1995), have suggested that limited and deteriorating resources in the wintering habitats have contributed to population declines for several species. Sutherland (1996) has provided an estimate of the decline in equilibrium populations from reduction of the breeding or non-breeding habitats based on the density dependencies of the respective birth and death rates. Seasonal population dynamics models addressing the problem of variable habitats are customarily based on the classical logistic differential equation of Verhulst (1845,1847),

$${\displaystyle \frac{dx}{dt}}=rx\left(1-x\right),$$

(1a)

describing the density-limited rate of growth of a population. Here, $x=N/K$ is the population density, N, relative to the carrying capacity, K, of the habitat, and r is the ”intrinsic” growth rate, usually associated with the difference between the birth and death rates, $b-d$. The $F\left(x\right)=\left(1-x\right)$ factor represents the density dependence of r, moderating otherwise exponential Malthusian growth. This expression has the analytical solution (Pearl 1927, Maynard Smith 1968, Pielou 1969, Acheson 1997),

$$x=1/\left[1-(1-x_0^{-1})e^{-rt}\right],$$

(1b)

where $x_0=N_0/K$ is the initial population density. For $x_0\ll 1$, this solution starts out with nearly exponential growth, but eventually saturates exactly at the carrying capacity, $x=1$, exhibiting the well-known ”sigmoidal” shape. Pearl (1927) demonstrated the utility of the logistic curve in fitting the growth of laboratory-cultured yeast, but similar logistic projections of the US population (saturating at ~200 million in year ~2000) from historic census data (Pearl, Reed and Kish 1940) have not been realized. Though extrapolation of any model fit to limited, noisy data is always risky, such failure likely results from inadequacies of the model.

The Verhulst logistic, Equation $\left(1\right)$, is a deterministic, single-species, single-mechanism, single-isolated-habitat, age-independent, spatially-homogeneous, instantaneous-response, first-order density-dependent model. These assumptions may restrict its applicability, and the literature is replete with extensions for statistical fluctuations of the parameters, predator-prey and competitive species interactions, age-specific birth and death rates, spatial diffusion, delayed response, and nonlinear density dependencies (see texts: Lotka 1956, Maynard-Smith 1968, 1974, Pielou 1969, Royama 1992, Krebs 1994, Kot 2001, Caswell 2006, and references therein). In this work, we primarily consider generalizations for multiple habitats and mechanisms encountered by migratory or resident seasonal breeders. We also examine threshold effects and possible instabilities due to the delayed response inherent in substituting discrete difference equations for differential equations, such as Equation $\left(1a\right)$. In a companion paper (Pine and Rappole 2008), we discuss the requirements for the validity of overall density-dependent population dynamics models for periodic breeders that exhibit age-structured birth and death rates.

2. Generalized Single-Habitat Logistic Models

The Verhulst logistic model, Equation $\left(1\right)$, has only two parameters, r and K, available to fit data, and K is just a scale factor. For increased flexibility, Richards (1959) and Gilpin and Ayala (1973) introduced an additional ”asymmetry parameter”, p, to the density dependence, in the form, $F\left(x\right)=\left(1-x^p\right)$. Herein, we refer to this p parameter as the ”saturation power”, since it affects the abruptness of the population saturation, as shown below. However, this modification, as well as the standard $p=1$ case, may have some difficulties for populations that exceed the carrying capacity because of the negative-going $\left(1-x^p\right)$ density-dependent factor. In particular, if the birth rate is less than the death rate, then extinction is predicted for $x_0<1$, as expected, but anomalous growth occurs for $x_0>1$. Indeed, for $r<0$ and $x_0>1$, the denominator in Equation $\left(1b\right)$ vanishes at $t=r^{-1}ln(1-x_0^{-1})$, leading to a population singularity (Kuno 1991, Royama 1992). For this reason, the range of validity is often restricted to $0\le x\le 1$. However, for multiple sequential habitats as considered here for periodic breeders, the population out of one habitat may exceed the carrying capacity of the next. To eliminate such anomalous growth for non-breeding seasons, we modify the logistic by simply allowing for separate density dependencies for the birth and death rates,

$${\displaystyle \frac{dx}{dt}}=x\left[b\left(x\right)-d\left(x\right)\right]=x\left[b_{0}F\left(x\right)-d_{0}G\left(x\right)\right].$$

(2)

Here, the birth-rate density-dependent factor, $F\left(x\right)$, is constructed to be unity for $x=0$ and to vanish or be greatly reduced at the carrying capacity, whereas the death-rate density-dependent factor, $G\left(x\right)$, increases with density from unity at $x=0$. Separation of the density dependencies of the birth and death rates is not a new idea of course (Williamson 1972, Sutherland 1996), but it is critical to the periodic breeder problem. Fretwell (1972) avoided the anomaly for ”short generation” species by ignoring the seasonal change in the intrinsic growth rate, r. Kuno (1991) and Royama (1992) invoked a negative effective K for $r<0$, which is unphysical, but reflects the dominance of $G\left(x\right)$ over $F\left(x\right)$ when deaths exceed births.

Some illustrative functional forms for $F\left(x\right)$ and $G\left(x\right)$ are considered here,

$$F\left(x<q\right)=\left[1-{\left(x/q\right)}^p\right],F\left(x\ge q\right)=0,$$

(3a)

$$F\left(x\right)=exp\left[-{\left(x/q\right)}^p\right],$$

(3b)

$$F\left(x\right)=1/\left[1+{\left(x/q\right)}^p\right],$$

(3c)

$$F\left(x\le q\right)=1,F\left(x>q\right)={\left(q/x\right)}^p,$$

(3d)

and

$$G\left(x\right)=1/\delta \left(x\right),\delta \left(x\right)=\left[1-{\left(x/q\right)}^p\right]\mathrm{or}\delta \left(x\right)=ϵ<1\mathrm{if}\delta \left(x\right)<ϵ$$

(4a)

$$G\left(x\right)=exp\left[{\left(x/q\right)}^p\right],$$

(4b)

$$G\left(x\right)=\left[1+{\left(x/q\right)}^p\right],$$

(4c)

$$G\left(x\le q\right)=1,G\left(x>q\right)=1+p[1-\left(q/x\right)].$$

(4d)

In these expressions, q is the ratio of any particular birth or death carrying capacity to an arbitrary K (usually the minimum value), and p is again the saturation power. Equation $\left(3a\right)$ is basically the Richards-Gilpin-Ayala expression, with the physical provision that the birth rates cannot be negative. Equations $\left(3b,c\right)$ are continuous positive functions for all $x\ge 0$, and have been studied previously in discrete-time, single-habitat population dynamics models (Ricker 1954, May and Oster 1976, Royama 1992, Leslie 1948, Beverton and Holt 1957, Thomas et al. 1980). Equation $\left(3d\right)$ is a hyperbola for $x>q$, piecewise continuous with unity for $x\le q$, assuring the physical requirement that $0\le F\left(x\right)\le 1$. Functions $\left(4a-c\right)$ for $G\left(x\right)$ are essentially the reciprocals of the corresponding functions $\left(3a-c\right)$ for $F\left(x\right)$, with prevention of a vanishing or negative denominator in $\left(4a\right)$. Although Equations $\left(3d\right)$ and $\left(4d\right)$ appear to be artificial constructs, they actually have a fairly straightforward biological interpretation. For example, if there are $K/2$ nest sites (assuming mating pairs) and any excess population is non-productive, but benign, then $\left(3d\right)$ with $p=1$ simply represents the reduction of the per capita birth rates. Higher values of p represent a faster reduction of birth rates, perhaps due to competition for the available nest sites or territories. If the excess population leaves in search of a better habitat, then $G\left(x\right)$can be represented by Equation $\left(4d\right)$ with p interpreted as the ratio of the emigration rate to death rate. Note that emigration in this case is equivalent to death for the habitat of interest. By contrast, immigration cannot be treated similarly to birth, as it need not depend on the existing population. The other functions may represent a distribution of quality or availability of nest sites or resources (Royama 1992). The $\left(1-x\right)$ factor in the Verhulst model is sometimes referred to as the ”available space” or ”unutilized opportunity” (Krebs 1994), and in the analogous ”epidemic” model (Acheson 1997) as the ”uninfected”.

Equation $\left(2\right)$ can be integrated for the population density at any time, provided an initial value, $x_0$. In general, the integration can be carried out numerically, very accurately using the fourth-order Runge-Kutta method (Press et al. 1992, Acheson 1997). For the special case when $F\left(x\right)=\left[1-{\left(x/q\right)}^p\right]$ and $G\left(x\right)=1$, and when the parameters are constant for a given time, t, there is an analytical solution to Equation $\left(2\right)$ in the form,

$$x=q/{\left[a-\left(a-{\left(x_0/q\right)}^{-p}\right)exp\left(-prt\right)\right]}^{1/p},$$

(5)

where $r=b_0-d_0$ and $a=b_0/r$. This solution follows by integrating the partial fraction $1/\left[x\left(1-a{x}^p\right)\right]=1/x+a{x}^{p-1}/\left(1-a{x}^p\right)$ after substituting $y=1-a{x}^p$ and $dy/dx=-ap{x}^{p-1}$. The standard analytical solution to the Verhulst logistic, Equation $\left(1b\right)$, is a special case of Equation $\left(5\right)$ for $a=p=q=1$. In the lower panel of Figure 1, we plot the time development of the population from Equation $\left(5\right)$ for several values of p with $q=1$ and $x_0=0.1$ and the corresponding density-dependent $b_{0}F\left(x\right)$ and $d_0$ shown in the upper panel. Note that the saturations occur where the net growth rates vanish at the crossings of the respective birth and death rates in the upper panel, below the $x=N/K=1$ nominal carrying capacity in this case. As p increases, the saturation approaches the carrying capacity, K, with a more abrupt inflection, and the growth portion more closely follows the density-independent exponential Malthusian curve.

Numerical integration for a birth density dependence given by Equation $\left(3d\right)$ and death by Equation $\left(4b\right)$ yields the results shown in Figure 2. In this case, the saturations occur at the crossings above the $x=1$ nominal carrying capacity, but again an increase in the saturation power, p, more closely approaches both the Malthusian growth and the $x=1$ limits with a more abrupt transition. For the $F\left(x\right)$ factors given by Equations $\left(3b\right)$ or $\left(3c\right)$, saturation may occur above or below the carrying capacity, depending on whether $d_{0}G\left(x\right)$ crosses the birth rates below or above $b_0/e$ or $b_0/2$, respectively. The saturation value is then the ”effective” or ”combined” carrying capacity, which depends not only on K, but on contributions from the distinct birth and death factors, F and G.

As indicated above, there may be more than one density-dependent mechanism operating simultaneously in the same habitat, controlling either the birth or the death rates. For example, birth rates may be limited by nesting sites and territorial behavior; death rates by food, water, shelter, predators, and communicable disease. Each of these mechanisms may be characterized by a different functional form, such as Equations $\left(3a-d\right)$ or $\left(4a-d\right)$, with specific carrying capacities and saturation powers. The overall effective density dependence is then given by the product of the various functions of the operative mechanisms. As might be expected, the mechanism with the minimum carrying capacity dominates, but the other mechanisms generally reduce the saturation population and may inhibit the growth rate. A reduced saturation value is evident from the upper panels of Figure 1 and Figure 2 since additional limiting resources increase the upward (downward) curvature of the death (birth) rates, thereby lowering the crossing density. Slower growth may be inferred from Figure 1 where it is seen that the population deviates from the Malthusian curve well below saturation for any saturation power, the anticipation being greater for lower p. For the mechanisms represented by Equation $\left(3d\right)$ shown in Figure 2 or by Equation $\left(4d\right)$, there is no deviation from Malthusian growth for $x<q$, where $F\left(x\right)$ or $G\left(x\right)=1$, so there is no anticipation, and any other simultaneous mechanism would be manifest only very near saturation.

3. Multiple-Habitat Models

In the above discussion, the parameters representing the birth and death rates, carrying capacities and saturation powers have been assumed to be constant in time. In general, though, they may vary explicitly with time, perhaps randomly due to weather, gradually due to habitat deterioration, or periodically due to seasonal breeding or migration. If the time dependence is known or is predictable, then it may be incorporated easily into the numerical integration of Equation $\left(2\right)$ (Press et al. 1992, Acheson 1997).

In the case of periodic breeders, the time dependence of the parameters generally may be expressed as a Fourier series. For example, Skellam (1966) has discussed a ”periodic normal distribution” of the form, $b\left(t,T\right)/b_0=1+2{\sum }_{n=1}^{\infty }exp\left(-n^{2}T/\tau \right)cos\left(2\pi nt/\tau \right)$, which has a single peak at $t=0$ of width scaled by T in a total period of $\tau$. However, here we take the parameters to be a sequential series of steps, being constant throughout the seasonal occupation of a particular habitat. This latter approach also permits a seasonal change of the density dependence functional form, e.g. Equations $\left(3\right)$ and $\left(4\right)$, if desired. In this case, we give a seasonal index label, m, to the parameters, $b_m$, $d_m$, $q_m$, $p_m$, and functions, $F_m$ and $G_m$, and integrate Equation $\left(2\right)$ over each season of duration ${\tau }_m$ separately, with an initial seasonal value obtained from the final value of the previous season (using $x_0$ for the first season). If there are M seasons in the total period time, $\tau$ (e.g. one year), then $\tau ={\sum }_{m=1}^M{\tau }_m$, and the m label is cyclical, $M+m=m$. Here, we refer to a ”season” as the dwell time in a habitat of substantially constant parameters, such habitats being geographically separate for migratory species, but which may coincide or overlap for resident species. We assume that the entire population either migrates or not, in order to avoid the complications associated with partial migration (Kaitala et al. 1993).

In Figure 3, we show the time development of the population using the generalized logistic, Equation $\left(2\right)$, for the case of two sequential habitats with ${\tau }_1={\tau }_2$, but only one breeding season. Here we have taken the birth rate to have the controlling density dependence (smallest carrying capacity) with the functions in Equations $\left(3a-d\right)$ as noted; Equation $\left(4b\right)$ is chosen for the density dependence of the death rates in both seasons. The sawtooth curves given by the solid lines represent the seasonal variations of the populations. The equilibrium population, $x_{eq}$, which is the average population over the cycle once saturation or steady state has been achieved, is given by

$$\sum _{m=1}^M\left[b_m{F}_m\left(x_{eq}\right)-d_m{G}_m\left(x_{eq}\right)\right]{\tau }_m=0.$$

(6)

This is a generalization for multiple habitats of arbitrary duration of the graphical analysis for a single season shown in the upper panels of Figure 1 and Figure 2 and for two seasons of equal duration with breeding in just one given by Sutherland (1996). Here, we see that saturation occurs over a wide range below and above the nominal carrying capacity, depending on the form of $F\left(x\right)$. If the density-dependent mechanism controlling the growth is dominated by the death rates, then the functional form for the birth rates makes little or no difference. If the death rate carrying capacity is lower in the breeding season, the seasonal population variations are proportionally smaller than shown in Figure 3, whereas if death in the non-breeding season dominates, the variations are larger.

4. Discrete Models

The above age-independent differential equation models imply that the progeny have the same fertility and mortality at birth as their progenitors. This is unrealistic for most complex organisms, and Fretwell (1972) has also considered the equilibrium populations for ”long-generation” species where the newborns are infertile during their birth season. In this case and for many species that abandon their eggs, the birth rate depends more on the adult population at the beginning of the breeding season than the total density at birth. In such circumstances, it may be appropriate to replace the derivative in Equation $\left(2\right)$ with a seasonal difference relation, $dx/dt\to (x_{m+1,n}-x_{m,n})/{\tau }_m$, yielding an iterative expression,

$$x_{m+1,n}=x_{m,n}\left[b_m{F}_m\left(x_{m,n}\right){\tau }_m+s_m\left(x_{m,n}\right)\right],$$

(7a)

where $s_m$ represents the seasonal survival,

$$s_m\left(x_{m,n}\right)=1-d_m{G}_m\left(x_{m,n}\right){\tau }_m.$$

(7b)

Here, $x_{m,n}$, is the density at the beginning of season m in breeding cycle n. Again, if there are M seasons per cycle, then $x_{M+1,n}=x_{1,n+1}$. The density-dependent factors, $F_m$ and $G_m$, may still be given by Equations $\left(3\right)$ and $\left(4\right)$, by substituting the seasonal carrying capacities, $q_m$, and saturation powers, $p_m$, and the $x_{m,n}$ for x. This approximation simplifies and speeds up the computations dramatically over the integration of the differential Equation $\left(2\right)$. Of course, this model also implies that the death rate depends on the initial seasonal population, rather than the existing density, which is more difficult to justify biologically. Also, Equation $\left(7b\right)$ may lead to unphysical negative populations if $d_m{G}_m\left(x_{m,n}\right){\tau }_m>1+b_m{F}_m\left(x_{m,n}\right){\tau }_m$. These difficulties may be alleviated if we replace Equation $\left(7b\right)$ with

$$s_m\left(x_{m,n}\right)=exp\left[-d_m{G}_m\left(x_{m,n}\right){\tau }_m\right],$$

(7c)

which represents the probability of survival and is always positive. Equation $\left(7c\right)$ approximates $\left(7b\right)$ if the probability of dying during the season is small, that is if $d_m{G}_m\left(x_{m,n}\right){\tau }_m\ll 1$. Moreover, this approximation is compatible with the well-known discrete age-structured single-habitat population model of Leslie (1945,1948), either in the limit of age-independent birth and death rates or in the circumstance of a stable age distribution. The generalizations for multiple habitats under which this correspondence holds are examined in a following paper (Pine and Rappole 2008).

In Figure 3, we superimpose the population dynamics for the discrete model, Equations $\left(7a,c\right)$, for comparison with the continuous differential model, Equation $\left(2\right)$, for the same seasonal parameters and functions. The unconnected symbols show the discrete results. Generally, the discrete points lag the changes in the comparable continuous curve, which tends to slow the early growth, but which may overshoot and exceed the saturation values somewhat. These lags result from the inherent delay in births and deaths in response to initial seasonal populations. If we make the same comparison for the single habitat case, shown in Figure 4, we find that saturation values for the discrete models are slightly higher than for the differential models owing to the approximation of Equation $\left(7c\right)$. However, such smooth monotonic behavior in a single habitat is not necessarily the consequence of this discrete model.

We show in Figure 5 the dynamics predicted for a four times larger birth rate with the same density factors as in Figure 4. Here we obtain steady, complex oscillation when $F\left(x\right)$ is given by Equation $\left(3a\right)$, damped oscillation for Equation $\left(3c\right)$, and aperiodic fluctuations for Equation $\left(3d\right)$. We note that the seasonal variations observed for multiple habitats, as seen for example in Figure 3, have a completely different physical origin than the ”inductive” oscillations seen in Figure 5. These inductive fluctuations are a consequence of overshoot due to the delayed response in the discrete model, and have been well studied for non-overlapping generations in a single habitat (q.v. May 1973, 1976, May and Oster 1976, Krebs 1994, Acheson 1997, Weisstein 1999). Here we see similar behavior for overlapping generations (May and Oster 1976, Kot 2001), with a significant qualitative dependence on the functional form of the density dependence.

The single-habitat fluctuations seen in the upper trace of Figure 5 for $F\left(x>1\right)={\left(q/x\right)}^p$ are not entirely chaotic since the displacements above and below the differential saturation value show a high degree of symmetry. These fluctuations require, not only a high birth rate, but the strong density dependence of the death rate. A constant death rate for the same high birth rate produces monotonic saturation. Also, such inductive instabilities may be repressed for multiple sequential habitats because of the shorter delays involved with seasonal updates. In fact, a large number of very short seasons would better approximate the continuous differential model in Equation $\left(2\right)$. Here we are ignoring explicit delayed density regulation effects (Maynard Smith 1968, 1974, Royama 1992, Kot 2001) which may arise, for example, for species with gestation, egg incubation, or seed dormancy longer than the cycle period. Arbitrary delay due to age-structured fertility is considered in the following paper (Pine and Rappole 2008).

5. Metamorphic and Spawning Species with Non-Overlapping Generations

Equations $\left(7a,c\right)$ can also represent discrete non-overlapping generations for metamorphic or spawning species if we let m count the life stages, such as larva, pupa, adult, and set $s_m=0$ for the breeding phase. Here, the prodigious output of many spawners is typically balanced by the low probability of survival of the pre-adult phases. This is consistent with the prior single-habitat, non-overlapping generations models, where the survival of the non-breeding phases is incorporated into the effective birth rate (sometimes denoted as the ”biotic potential”) in the adult phase. Splitting out the pre-adult phases, as done here, allows us to apportion the carrying capacities and mechanisms to the appropriate phases. In Figure 6, we show a hypothetical two-stage metamorphic species where the high death rate of the adults is sufficient to eliminate breeding survivors, but enough juveniles survive to allow overall growth until a saturation is reached near the non-breeder carrying capacity. We note that the differential model, Equation $\left(2\right)$, with a higher death than birth rate, just leads to extinction, so is not valid for metamorphic species. The life stage problem has been addressed previously using matrix formulations (Leslie 1945, DeAngelis et al. 1980, Caswell 2006), but this is not necessary if all births transition to the next stage in the cycle, as shown here for metamorphic species.

6. Threshold Effects

Some species may require a critical density, $x_c$, to survive or thrive, as discussed by Allee et al. (1949). This critical density could arise from the necessity to find a mate, establish a breeding colony, hunt cooperatively, school or flock, etc. Such threshold effects can be incorporated into the multi-habitat models above if the appropriate birth or death rates, $b_m$ or $d_m$, are respectively multiplied or divided by another density-dependent factor, $\Phi \left(x\right)$, that increases from a small residual value, $ϵ<1$, to 1 as x transitions through $x_c$. $ϵ$ should be chosen so that the death rate exceeds the birth rate, for example, $ϵ=0$ for breeding thresholds or $ϵ=d_m/(b_m+d_m)$ for survival thresholds. Some standard functions with this property are given by

$$\Phi \left(x\le x_c\right)=ϵ,\Phi \left(x>x_c\right)=1,$$

(8a)

$$\Phi \left(x\right)=1-\left(1-ϵ\right)exp\left[-{\left(x/x_c\right)}^{\eta }\right],$$

(8b)

$$\Phi \left(x\right)={\displaystyle \frac{1+ϵ}{2}}+{\displaystyle \frac{1-ϵ}{\pi }}arctan\left[{\displaystyle \frac{x-x_c}{\gamma }}\right],$$

(8c)

$$\Phi \left(x\right)={\displaystyle \frac{1+ϵ}{2}}+{\displaystyle \frac{1-ϵ}{2}}erf\left[{\displaystyle \frac{x-x_c}{\sigma }}\right].$$

(8d)

Equation $\left(8a\right)$ is a step function at $x=x_c$, representing a sharp, definite threshold. This modification does not change the dynamics discussed above except when the populations dip below $x_c$, whereupon they eventually go extinct. Equation $\left(8b\right)$ represents a more gradual transition, characterized by the $\eta$ parameter. Equations $\left(8c,d\right)$ are more probabilistic descriptions of this viability, representing a spread of effective critical densities given by a Lorentzian distribution, $\left[1/\left[1+{\left(x-x_c\right)}^2/{\gamma }^2\right]\right]$, or a Gaussian, $\left[exp\left[-{\left(x-x_c\right)}^2/{\sigma }^2\right]\right]$, respectively. Because the modifying functions in Equations $\left(8b-d\right)$ extend beyond the critical density, they may have a quantitative effect on the population dynamics, but the qualitative features discussed above remain. Such threshold effects may select against species whose population fluctuates wildly due to overbreeding (May and Oster 1976, Thomas et al. 1980).

7. Conclusions

We have shown, by the simple expediency of separating the density dependencies of the birth and death rates, that we can avoid the anomalous, unphysical behavior of the customary logistic model, Equations $\left(1a,b\right)$, for populations in excess of the carrying capacity. This procedure is essential for modeling the effects of sequential habitat variation encountered by seasonal breeders. The biological interpretation of the controlling density-dependent birth and death mechanisms is thereby more transparent, even in a single habitat. Simultaneously operative limiting mechanisms, including threshold effects, can be treated by multiplicative density-dependent functions. The disadvantages include numerical integration in place of analytical solutions and a great increase in the number of model parameters, some of which may be indeterminate in fitting population dynamics measurements. Some parameters can be specified independently by methods such as sampling, tagging and tracking. We have also demonstrated the deviations between continuous differential and discrete difference models, leading to time lags and possible inductive instabilities in the latter case.

These models may be useful to evaluate the factors limiting the abundance and distribution of migrating and resident seasonal species. In particular, they indicate how populations of migratory birds may be limited by non-breeding season habitat availability. More realistic behavior may require the inclusion of several other phenomena neglected here, such as immigration, stochastic processes and multi-species interactions. Application of the models for age-structured breeders may also be restricted, as shown for periodic breeders in an accompanying study (Pine and Rappole 2008).