Individual-based predator-prey model for biological coevolution: Fluctuations, stability, and community structure

We study an individual-based predator-prey model of biological coevolution, using linear stability analysis and large-scale kinetic Monte Carlo simulations. The model exhibits approximate $1∕f$ noise in diversity and population-size fluctuations, and it generates a sequence of quasisteady communities in the form of simple food webs. These communities are quite resilient toward the loss of one or a few species, which is reflected in different power-law exponents for the durations of communities and the lifetimes of species. The exponent for the former is near $-1$, while the latter is close to $-2$. Statistical characteristics of the evolving communities, including degree (predator and prey) distributions and proportions of basal, intermediate, and top species, compare reasonably with data for real food webs.

Figure 1

(Color online) Time series for the Shannon-Wiener diversity index (a) and the population sizes (b). The jagged curves in the background show data sampled every 2048 generations to show the rapid fluctuations, while the curves in contrasting color (brightness) in the foreground are running averages over 524 288 generations, emphasizing the slower fluctuations. The three sets of curves represent all species (black/light gray; black/yellow online), producer species (light gray/dark gray; green/violet online), and consumer species (medium gray/light gray; red/cyan online).

Figure 2

(Color online) (a) Species abundance distributions (SADs), normalized to the number of species. Solid curves represent “full connected communities,” and dashed curves represent “core communities,” both extracted as described in the text. Both were sampled every 256 generations. The data were averaged over 12 independent simulation runs, and the error bars represent standard errors, based on the differences between runs. The dot-dashed curve is a fit to the curve describing all core species, using Eq. (14) with parameters $C=13.569$, $\mu =0.00210967$, and $\beta =1.89338$, interpolating between log-series and log-normal forms. (b) Histograms of the full-community species richness, the species richness of the core communities, and the Shannon-Wiener diversity for the sampled communities. The latter is seen to be an excellent approximation for the species richness of the core communities.

Figure 3

(Color online) Time series for the same simulation run shown in Fig. 1, displaying quantities that measure the evolutionary activity: the magnitude of the logarithmic derivative of the Shannon-Wiener diversity index $\mid dS∕dt\mid$ (averaged over 16 generations) (a), and the size of extinctions per generation (b). The horizontal lines in each part (solid, long-dashed, and short-dashed) indicate different cutoff levels used to define quiet and active periods as discussed in Sec. 4D.

Figure 4

(Color online) PSDs for the Shannon-Wiener diversities (a), population sizes (b), and the extinction sizes and number of species going extinct per generation (c), all averaged over 12 independent runs of $2^{25}$ generations. The error bars are standard errors, estimated from the variations between the individual runs. The dashed straight line with slope $-1.29$ in (a) is a weighted fit to the PSD for the overall diversity over the whole frequency range. The dashed straight lines in (b) and (c) are guides to the eye with the same slope.

Figure 5

(Color online) Histograms for the lifetimes of individual species. The dashed, straight line is a guide to the eye, corresponding to a $t^{-2}$ power law. Data averaged over 12 independent simulation runs.

Figure 6

(Color online) (a) Histogram of the logarithmic derivative of the overall diversity $dS∕dt$. The data were averaged over 16 generations for each run and then over 12 independent runs. The dashed curve is a Gaussian fit. The time series of the magnitude of this quantity for one particular simulation run is shown in Fig. 3a. (b) Histograms for the durations of active and quiet periods for various values of the cutoff $y_c$ for $\mid dS∕dt\mid$, averaged over 12 independent runs. The dashed, straight line with slope $-1.07$ is a weighted fit to the histogram for the quiet periods for $y_c=0.010$ between 10 and ${10}^6$ generations. The steep histogram curve with only three data points between 10 and 100 generations corresponds to the active periods, which are always very short.

Figure 7

(Color online) Histograms of the durations of quiet periods, obtained from the time series of extinction sizes, an example of which is shown in Fig. 3b. Data averaged over 12 independent simulation runs. The straight, dashed line with slope $-1.07$ is a guide to the eye, based on the fit to the QSS duration distribution based on $\mid dS∕dt\mid$, shown in Fig. 6b.

Figure 8

(Color online) Major populated species shown vs time for the same simulation run shown in Figs. 1, 3. The horizontal lines correspond to the species label, and the gray scale (color online) to the population size. Black: $n_I>1000$. Gray (red online): $n_I∊[101,1000]$. (a) Producers. (b) Consumers.

Figure 9

(Color online) Food webs representing QSS core communities for the simulation shown in Figs. 1, 3, 8 at times near $22\times {10}^6$ generations (a), near $27\times {10}^6$ generations (b), and at the end of the simulation near $34\times {10}^6$ generations (c). The core communities were identified as described in Sec. 4A. The thickness and head size of the arrows correspond to the magnitude of $M_{IJ}$, and the area of the circles to the stationary population size, as calculated analytically from Eq. (5). Light gray (red online) lines without arrowheads connect nearest neighbors in genotype space. Species No. 609477 (label medium gray, red online) persists from the first to the last snapshot, (a) to (c). Species No. 682497, No. 681473, and No. 943617 (labels dark gray, blue online) persist from (a) to (b), while species No. 85703 (label light gray, green online) persists from (b) to (c).

Figure 10

(Color online) Histograms of the multiplication ratio $n_i(t+1)∕n_i\left(t\right)$ [exponential of the invasion fitness, Eq. (12)] for low-density invaders. Each curve is averaged over 14 different core communities. Black: nearest neighbors in genotype space against stable QSS communities. Dark gray (red online): All species not in the community against stable QSS communities. Light gray (orange online): All species not in the community against random feasible communities. Inset: The part of the distributions for potentially successful invaders $n_i(t+1)∕n_i\left(t\right)>1$ with the $y$ axis on logarithmic scale. The tails appear nearly exponential.

Figure 11

(Color online) Log-linear plot of histograms of the birth cost $b_I$ (peaked to the left) and producer coupling to the external resource ${\eta }_I$ (peaked to the right), for members of QSS core communities (shaded areas) and full communities (lines only). Both parameters are selected away from the uniform distributions of the full species pool (horizontal dashed line).

Figure 12

(Color online) Histograms of community-level parameters. (a) The interspecies interaction strengths $\mid M_{IJ}\mid$ for QSS core (shaded) and full (line only) communities. (b) The effective interaction strength $\Theta$ (left) and resource coupling $\mathcal{E}$ (right) for QSS core communities, defined in Eq. (7).

Figure 13

(Color online) Histograms for the linkage density $z=L∕S$ (main part) and the directed connectance $C=L∕S^2$ (inset), for the model communities. The histograms for the core communities are shaded.

Figure 14

(Color online) (a) Unscaled degree distributions for the QSS core communities, showing the number of prey (solid), the number of predators (dashed), and the total degree (dot-dashed). The distributions were individually normalized for each community, and then averaged over all communities. The error bars are standard errors, based on the spread between the 12 simulation runs. Average linkage density is $\overline{z}=1.17$. The inset on log-linear scale shows that these distributions decay at least exponentially with increasing argument. (b) Same as (a), for the full communities. Average linkage density is $\overline{z}=2.58$. (c) Histograms over individual sampled communities of the correlation coefficient between a species’ number of predators and prey. Black filled with light gray (yellow online): QSS core communities. Medium gray (red online): full, connected communities. Dark gray (blue online) with drop lines: empirical data for seventeen natural food webs.

Figure 15

(Color online) Scaled, cumulative degree distributions for the model, compared with results for the 17 empirical food webs. Data for individual, real food webs are shown in the background as isolated data points. Average data for the simulated and real webs are shown as data points with error bars indicating standard error, connected by heavy lines. Dark gray with circles (blue online): simulated core communities. Medium gray with squares (red online): simulated full communities. Black with diamonds: empirical communities. Details of the averaging are discussed in the text. (a) Number of predators (vulnerability, or outdegree). (b) Number of prey (generality, or indegree). (c) Number of predators plus prey (total degree).

Figure 16

(Color online) Histograms for the proportions of basal, intermediate, and top species. Basal: thin black line filled with light gray (yellow online). Intermediate: line of intermediate thickness and gray scale (red online) without filling. Top: thick, dark gray (blue online) line without filling. The vertical, dashed lines in matching gray scale (color) and thickness represent the mean proportion for each class. (a) Core communities. Mean proportions: 0.33 for basal, 0.31 for intermediate, and 0.38 for top species. (b) Full communities. Mean proportions: 0.14 for basal, 0.76 for intermediate, and 0.13 for top species.

Figure 17

(Color online) Normalized correlation functions for the matrix elements $M_{IJ}$ for the two schemes described in the text, each based on a single realization of $\mathbf{M}$ for $L=6$, 8, and 9. (a) The scheme described in Ref. 17. (b) The modified scheme introduced here.

Figure 18

(Color online) Histogram for the pseudorandom $M_{IJ}$ produced by the scheme introduced here for $L=9$, corresponding to 262 144 individual matrix elements. The straight, black lines represent the triangular shape of the theoretical distribution.