Population Dynamics Simulations: Predator-Prey Systems

The Lotka-Volterra equations, independently derived by Alfred Lotka in 1925 and Vito Volterra in 1926, are among the most studied and celebrated models in mathematical biology. They describe the interaction between a prey population (N) and a predator population (P): dN/dt = αN − βNP; dP/dt = δNP − γP. In these equations, α is the prey growth rate in the absence of predation, β is the per-capita predation rate, δ is the predator's conversion efficiency (the fraction of consumed prey converted into new predators), and γ is the predator death rate in the absence of prey.

The qualitative prediction of the model is periodic oscillation in both populations, with predator cycles lagging behind prey cycles by approximately one quarter-period. This prediction appears to be borne out by the Hudson's Bay Company's fur trade records, which show striking 9-11 year cycles in lynx and snowshoe hare pelt counts from 1845 to 1935, one of the most frequently cited empirical validations of a theoretical biology model.

However, the classical Lotka-Volterra model has well-known limitations that motivate extensions. The model assumes a linear functional response: predation rate is exactly proportional to prey density, so a predator in an environment with twice as many prey will kill exactly twice as many per unit time. This is unrealistic at high prey densities, where predators become satiated and the marginal kill rate diminishes. C.S. Holling's 1959 experiments established a more realistic Type II functional response: predation rate rises with prey density but asymptotes toward a maximum handling-rate-limited value. At very high prey densities, the Type III response applies, where there is also a sigmoidal "switching" behavior as predators focus on abundant prey types. These extensions produce qualitatively different dynamics: while the classical model's trajectories are closed orbits (conservative), the Type II system can exhibit true limit cycle attractors or, under appropriate parameter values, deterministic chaos.

We numerically integrate both the classical Lotka-Volterra system and its Holling Type II extension (Rosenzweig-MacArthur model) using a 4th-order Runge-Kutta algorithm with adaptive step size control (Dormand-Prince method) to ensure numerical accuracy near the system's fixed points. Initial conditions and parameters are estimated by fitting the model to the Hudson's Bay Company pelt records (MacLulich 1937, reanalyzed by Elton and Nicholson 1942), which serve as the empirical benchmark.

Parameter estimation minimizes the sum of squared distances in phase space between the model trajectory and the empirical (hare, lynx) time series, using a genetic algorithm for global search followed by local refinement with L-BFGS-B. We then map the parameter space systematically to characterize the qualitative behavior of the Rosenzweig-MacArthur model across a grid of carrying capacity K values and Type II saturation constants h, computing Lyapunov exponents to distinguish stable equilibria (all exponents negative), limit cycles (one zero exponent), and chaos (positive maximal Lyapunov exponent). Phase portraits are computed for representative parameter values in each regime.. Phase portraits are computed for representative parameter values in each regime.