New paper: Ecological modeling of paleocommunity food webs

Roopnarine, P. D. 2009. Ecological modeling of paleocommunity food webs. in G. Dietl and K. Flessa, eds., Conservation Paleobiology, The Paleontological Society Papers, 15: 195-220.

Find the paper here:
http://zeus.calacademy.org/roopnarine/Selected_Publications/Roopnarine_09.pdf
or here
http://zeus.calacademy.org/publications/

Prey dynamics

The figure here is very similar to the one in the previous post, but these results are for the guild of shallow infaunal suspension feeders (primarily clams). The main difference is the more regular increase in the number of species that become extinct as the perturbation magnitude () increases. Another interesting note is that this guild is not the only driver, or any driver at all, of the behaviours exhibited by the guild of predators. Those predators may or may not prey on members of this guild, and also have an array of prey in other guilds. So the oscillatory behaviour seen at higher perturbation levels is probably system-wide. And it is system-wide because of indirect effects via network links. One wonders what a summary of the results would look like, and what the implications are for individual species population dynamics.

1. For example, even at a very low perturbation level, maximum sustainable population sizes oscillate wildly before settling down to a new stable state (which can in fact be the initial one, or zero, indicating extinction). One would assume that population sizes would follow this trend, if the timescales of the perturbation and population growth were sufficiently close. What if they are not? How does this affect what one would actually observe for a given species?
2. What is the distribution of stable states over the perturbation range? Are the oscillations observed at high perturbation level convergent, i.e. if run long enough they would also settle to a new stable state? Or are they asymptotic, but never settle down, or settle to two alternative states? One way to find out would be to simply run the series for many additional steps. Another would be model the oscillations themselves, and see if the convergence is linear or asymptotic. And what is the perturbation range of the bifurcations? At what point do we begin to observe oscillation/bifurcation, and is it synchronous throughout the community? Only one way to find out, but I’ll probably have to write some Sed/Awk or Perl scripts to handle these large datafiles.

Lighting up an ecosystem

One of the final pieces needed to explain the critical/threshold point in a bottom-up CEG perturbation is an understanding of which species become extinct, and what the species dynamics look like during the cascade. Therefore, what I’ve done is to modify the basic simulation to capture the demographic properties (technically, the carrying capacities) of each species; results in HUGE output files. Shown here in this figure are the species dynamics for the guild of shallow infaunal carnivores (e.g. naticid snails) at three different perturbation levels. Note that the levels correspond to a low secondary extinction response, and the two critical points identified earlier. The top row of figures plot the dynamics of surviving species, and lower show those of the species which become extinct. The community dynamics were recorded for 250 steps beyond the initial perturbation. The first thing to note is that species become extinct very quickly. Beginning K for each species is standardized at 1, and the species that become extinct have, on average, lower in-degrees, i.e. lower numbers of prey, than do surviving species (statistical tests to follow later). That result matches expectations of the CEG combinatoric model.

The other thing to note is that at the low and mid-perturbation levels ( 0.2 and 0.55), species’ K respond immediately to a perturbation of the producer guilds, oscillate for several steps, but eventually settle down to a new stable K. This is a transition to new stable states for the species populations. At the perturbation level which coincides with the major critical point of secondary extinction, however, there is no indication that the species ever settle to a new state. Instead, there seems to be bifurcation and subsequent alternation between two alternative stable states; the species are lit up by the disturbance (each species is given a different colour in the corresponding plot for easy distinction). The series is quasi-periodic though, in that the system never returns to quite the same point on alternating steps. It is possible that the series eventually converge to a single, or two, stable points, but the current data cannot address that. Therefore, I’ll next repeat this simulation, but extend the data capture to 1,000 steps (2Gb file). Hopefully that will give some indication of whether the series are converging, diverging, are stable, or perhaps chaotic. Also, the results for the other guilds need to be examined.

Recipe for building a scientific model

The following is the second part of the excerpt (see previous post) on model-building. Apologies for the paleontological and evolutionary biology jargon to any casual readers who might stumble across this. I’ll try to embellish and explain another time, or just leave me a question. Same for the citations.

“…the supreme goal of all theory is to make the irreducible basic elements as simple and as few as possible without having to surrender the adequate representation of a single datum of experience.” (Einstein, 1934). There are numerous views and opinions on the purpose for constructing a scientific model, but they all fall roughly between the end-points on a scale of general scientific goals: problem-solving (Grimm and Railsback, 2005) and constructing an underlying theory based on observations (Randall, 2005). For example, Darwin developed the qualitative model of natural selection to explain a theory of evolution that he deduced underlay many of his observations of the natural world. In the process he assumed the existence of a mechanism for the inheritance of traits. Evolutionary biologists of the early 20th century were subsequently faced with reconciling the problem between Mendelian genetics and concepts of natural selection at that time. Models of quantitative and population genetics were developed that went a long way toward reconciliation (Fisher, Wright, Haldane). It is not difficult to appreciate differences between the models and presentations of Darwin, and for example, Fisher, and that brings us to the first guideline.

A model should be appropriate to the problem or question being addressed. Darwin set out to explain patterns which he had observed, and his work is remarkable in its weaving together of so many major themes into the theory of evolution, including inheritance of traits, population growth and the struggle for existence, competition within and among species, and the history of life as recorded in the fossil record. Formulating his theory did not require a detailed understanding of any of these themes comparable to what we have available to us today. Moving beyond Darwin, however, in order to explain observations, test the theory of natural selection, and search for underlying mechanisms, requires the construction of far more specific models, and that is precisely what Fisher did. Darwin presented a convincing case of modification with descent, but it was up to Fisher and others to provide adequate models of natural selection capable of explaining the data. Fisher provided a specific and quantitative model explaining how natural selection can be based on the Mendelian “particulate character of the hereditary elements” (Fisher, 1958). This and related models, plus the later revolution in molecular biology, were crucial to our current understanding of natural selection as an evolutionary mechanism. Nevertheless, the details of these models are not at all essential when incorporating evolutionary theory into many biological models. For example, Sepkoski (1984) did not need them for his kinetic models of Phanerozoic diversity, even though all arguments of taxon origination are based on evolution. How then does one decide on the proper ingredients for a model?

Define the domain of the model. Presumably, all things in nature are connected, but most of these connections are never relevant to understanding the problem immediately at hand. One should define, at the outset, those connections or relationships which are causal, and those that are contingent (Bohm, 1957). To paraphrase an example from Bohm, one can quite accurately deduce a law of gravitation by dropping balls of various sizes and composition. Now drop a sheet of paper. Does the difference in motion suggest a new law, or inadequacy of the old one? The paper does eventually meet the ground, and precisely because of the same law of gravitation to which the balls are subject. Its different motion is a result, of course, of air resistance. If your goal is to deduce a law of gravity, then one realizes that gravity here is causal in the net journey of the paper, but the specific steps taken along the journey are contingent upon the sheet’s air resistance. Air resistance is irrelevant to a law of gravity. If your goal, on the other hand, is to model the motion of sheets of paper, then accounting for gravity alone is not sufficient and air resistance is indeed relevant.

Constrain model specificity. Getting back to the balls for a moment, is air resistance relevant to them? If one’s measurements were precise enough, then you would likely discover differences based on different surface textures of the balls, perhaps variation over time corresponding to air temperature, and so on. But those details should be unimportant when describing the motion of a dropped ball, because mass and sphericity dominate minor differences in ball composition. Therefore, the addition of air temperature as a causal parameter is unnecessary. There is often a temptation and tendency to strive for realism in biological models to the point at which the model becomes too complex. This point can be recognized either when the addition or variation of parameters fail to alter model output in a predictable manner, or variation of model output cannot be explained directly as a result of parameter variation. An accurate simulation of reality is not a model unless one can point to exactly why the simulation reproduces all the aspects that make it realistic. Simulations are not necessarily models.

Having stated that, a major caveat must be given immediately. Ecological systems are mostly complex, that is, they are the result of the interactions of numerous independently acting or semi-independent entities. Simple models of those entities can sometimes produce complicated results, particularly when causal relationships within the model are nonlinear (as we will see later on). The behavior of a system of entities can be modeled as a collection of simple models, or ignore that lower level altogether and treat the entire system as an entity. These approaches are not likely to yield the same results. The former can quickly become a hopelessly realistic simulation, while the latter simply cannot be interpreted at a level relevant to individual entities. A proper model will search for a middle path, and though the search can be as much art as it is science, that’s where the fun is.

Scientific models

The following is a little excerpt from a paper that I’m writing on the role of food web models (focusing on mine) in the developing field of conservation paleobiology.

A scientific model is a set of formal statements which are hypothesized to represent or explain some aspect of the natural world. The natural world is complex, and models are necessarily simplifications to some degree. They often reflect educated guesses of the relationship, or nature of a relationship between known processes or observed phenomena. Models do, however, allow us to make predictions about the natural world. Furthermore, models allow us to explore hypothetical possibilities, even when our data are insufficient or incomplete. Conservation paleobiology is concerned with the application of paleontological ideas and data to conservation biology, that is, understanding the circumstances under which species are threatened with extinction, become extinct, or survive those circumstances. These are the concerns that our models must address, and to be of use to conservation biology, they must make predictions of future biodiversity events on the basis of the fossil record.

Scientific models need not be mathematical, but in this paper I will focus on mathematical models for two reasons. First, the paper is intended to be instructive about the creation and use of models in conservation paleobiology. Almost any such model must be expected to contain statements of ecology and paleoecology, and mathematical modeling in ecology has matured into a rich and useful discipline. Therefore, “model” in this paper will always refer to a mathematical formulation. Second, when formulated mathematically, a model gives us much freedom to question our hypotheses and explore implications and consequences outside the bounds of empirical observations. This is particularly useful when data are unknown or incomplete, or the hypothesis itself is incomplete. Mathematical modeling should not be taken, however, as an exclusive alternative to empiricism. The two approaches to science are complementary: a Platonic approach that presumes the existence of fundamental truths, and the Aristotelian approach rooted in empirical observation. If there is any order in Nature, and scientific experience tells us that there is, then our models can be based on our understanding of fundamental principles, and should be expected to be congruent with empirical data to the extent that the model reflects the entirety of the processes producing the data. Incongruence will often require a reformulation of the model, or even revision of the principles themselves, a task made easier and more precise when the model has been stated and communicated mathematically. Alternatively, one should also always bear in mind the division between data and interpretation. While data perhaps seldom lie, scientific theories are based on histories of important misinterpretation, and there are numerous examples of models that predicted as yet unobserved data or phenomena, while being at odds with interpretations of existing data.

Tipping point II

Red ellipses represent the two threshold points.

The story so far: We have a food web of a shallow coastal marine community from the Late Miocene of the Dominican republic. The metanetwork comprises 29 guilds, 139 guild-level links, and 130 species. A perturbation of the system, where all three primary producer guilds plus detritus were systematically and incrementally removed from derived species-level networks, results in the typical CEG result: that is, a relatively flat and low level of secondary extinction () over a broad range of perturbation magnitude (), succeeded by a rapid transition to a state of high secondary extinction. In fact, for this community, there are two transitions. The first occurs at , and represents a very minor but secular increase in . The second transition occurs at and represents a catastrophic increase in . Topological-only perturbation of the system makes it very clear that these transitions correspond exactly to two stages of the perturbation: First, the complete extinction or removal of benthic autotrophs and complete disruption of the particulate detritus supply. The second and greater transition occurs at the complete extinction of the benthic macroalgae and macrophytes. Accompanying the second transition is the complete extinction of the benthic herbivore guilds which specialize on the macroalgae and macrophytes (and derived detritus), comprising families such as the Phasianellidae, Cerithiidae, Vitrinellidae, Haminoeidae and Retusidae. This is accompanied by extinction of species in other more generalist guilds that include macroalgae in their diet.

Extinction of those heterotrophic taxa is not itself the cause of the major tipping point though. Simulations where the perturbation is specifically removal of these herbivores result in very low levels of secondary extinction, with no tipping point or threshold. The obvious question then is, why does extinction of the macroalgae drive the system to a new state? The qualitative answer is that the complete loss of this resource, and the bottom-up propagation to the herbivores, in turn cause intense top-down cascades of compensatory responses from higher level consumers. These cascades propagate throughout the network, even to the remaining source of production, the phytoplankton. The result is a tremendous loss of species. A very curious thing, however, is that phytoplankton productivity in the network is almost 3 times greater than macroalgal productivity, reflecting the much greater diversity of planktivores. So why does the collapse coincide with loss of the macroalgae?

Perturbation of macroalgae (top) and phytoplankton (bottom), and resulting loss of autotrophic productivity because of top-down effects (left), and secondary extinction of heterotrophic species (right).

I performed two separate perturbations to answer this question. First, I perturbed the system by removing macroalgae only, and second by removing phytoplankton only. The top row of the second figure shows the results of the first experiment. Secondary loss of autotrophic resources (left column) as a result of top-down effects is effectively zero. Secondary extinction of heterotrophs (right column) is significant but not dramatic. There is a mild increase in the region of , which represents the loss of the specialized herbivore guilds. Removing phytoplankton had a more dramatic impact, reflecting the greater overall dependence of the community on phyloplankton resources. There is a clear threshold, occurring at approximately . At this point, resource loss to the community is great enough to trigger the catastrophic top-down cascades and feedback within the network. Therefore, it seems that in the previous experiment, where all resource guilds were perturbed, the complete loss of macroalgae triggers the top-down cascades and compensatory feedback that in turn deplete phytoplankton resources to the point where the system transitions to a higher state of secondary extinction. This conclusion is supported by the fact that when all producer guilds are perturbed, the contribution or perturbation of phytoplankton at the tipping point is 38%, whereas when only phytoplankton are perturbedm the tipping point occurs at 50%.

Some closing observations:

1. Topological analyses of network vulnerabilities are likely to underestimate the severity of link losses when those links have variable interaction strengths, and the nodes have varying properties. In the case of a biological community, species could and are likely to alter interaction strengths to compensate for lost resources (i.e. links). Topological vulnerability analyses should be well suited for networks with static properties, perhaps such as power grids and the internet (though I’m no expert here!), but are ill-suited for dynamic networks, such as those describing transportation, metabolic/physiologic and ecologic systems.
2. An hierarchically structured, directed network such as an ecological community should be resistant to a broad array of random perturbations. This is a function of both the underlying link distributions (as already understood in the case of static networks or graphs), as well as the compensatory abilities of consumer species, and the variance of dietary breadth. The network is, however, vulnerable to the loss of highly linked nodes. Here I am referring specifically to basal, autotrophic nodes, and not necessarily keystone consumer species. Not all autotrophic nodes are equal, however, as shown in the above results. Nevertheless, because of the complexity of the species interactions and the hierarchical divisions of ecological functions, there should be strong nonlinearities in the network responses. This is borne out by the differences between the topological-only and fully dynamic simulation results. The nonlinearities are expressed as two or more alternative states of secondary extinction, separated by rather sharply defined thresholds of perturbation. I can think of no way in which to analytically predict the threshold points, but heuristically I would argue that they should exist in every ecological community.
3. Perturbation of top-level consumers are observed in nature to often result in top-down cascading effects, compatible with such notions as keystone predators. I will show in later results that the CEG model captures all this. The results will also show, however, that while top-down effects can be locally catastrophic, i.e. for individual species or groups of closely linked species, they are never globally catastrophic in the manner in which bottom-up perturbations are. This conclusion has implications for understanding the role of ecological collapse in large scale extinctions observed in the fossil record. It also has implications for the ongoing biodiversity crisis, where species far removed from the “tops” of food webs are increasingly threatened by climate change and habitat destruction.
4. An close examination of many of the results presented in this blog will show apparent “bifurcation” of the results, e.g. beyond the threshold point in the lower right graph above. These observations suggest that there is more than one type of species-level network that can be derived from the same metanetwork. So, while the higher-level organization of the community is the same, networks are being generated that vary enough in their interspecific link topologies to yield very different responses to the same level of perturbation. I believe that this is a statistical property of the underlying trophic link distributions and the resulting multinomial probabilities from which the species-level networks are drawn stochastically. In the case of the above results, where one set of networks is significantly more resistant than the other (i.e. they have a much higher tipping point), this mathematical feature of the model is not likely to be of great relevance ecologically. That is because the lower threshold is already so high, in this case, 50% shutdown of primary productivity. Those are catastrophic environmental conditions and would occur with very low frequency in nature on a large scale. There are cases, however, such as the Early Triassic Lystrosaurus zone community, where there seem to be multiple alternative states at very low perturbation levels. Those communities would very likely have experienced frequent low-level perturbations, and then one has to consider whether: (1) this feature of the model is a mathematical artifact, in which case one wonders about the constraints necessary to prohibit it in nature, or (2) the feature is real, and then one wonders how species within a community cope with such a situation.

Progress

Ah, finally back to work. The lungs seem to cooperating again, and while not 100%, boy, it felt good to be back to the Academy. Also, being back meant that I was able to finally begin testing some of these couch-bound speculations from the past couple of weeks!

Shown in this first figure is a metanetwork representation of a Late Miocene, shallow marine community from the Dominica Republic. These data were compiled by one of my graduate students, Rachel Hertog. Recall that each sphere represents a set of one or more species that potentially share the same predators and prey. The guilds are colour-coded, but we’ll ignore that for now. The links between guilds represent sets of trophic interactions. This paleocommunity has 29 guilds and 139 guild-level links. The guilds range from phytoplankton to epifaunal benthic carnivores to pelagic carnivorous fish. There are 130 species in the community.

I ran topological-only and fully dynamic simulations of bottom-up perturbation on species-level networks derived from this metanetwork. The perturbation is a progressive reduction of primary productivity, implemented as a progressive reduction in the size of all four primary production guilds. The second figure shows the results. The fully dynamic results (in yellow) exhibit the typical CEG result. The topological-only results are shown in aqua. As expected, topological extinction underestimates the scale of extinction possible, and follows a predicted “exponential-type” of increase. Notice, however, that the pattern takes a little “hop” at an approximate perturbation magnitude of 0.67. Note, also, that it is at precisely this point that the dynamic results show the typically rapid increase in secondary extinction level.

Topological-only and dynamic results of bottom-up perturbation.

An examination of the topological results reveal that the hop is due to the complete extinction, at that point, of two guilds: epifaunal herbivores, and shallow infaunal herbivores. These two guilds consume the macroalgae/seagrass guild exclusively. That this is the predicted point of extinction can be checked analytically. Remarkably, it explains the discontinuous/catastrophic increase seen in the dynamic results. There are five guilds that include these herbivore guilds as prey, and they are all carnivorous or omnivorous macroinvertebrates and fish, shallow infaunal, epifaunal and pelagic. All these consumers have a wide array of prey guilds, but the loss of the two herbivore guilds represent a significant enough loss of resources that the compensation of the predators, represented in the model as increases of interaction strength with remaining prey, causes top-down cascades strong enough to in turn cause the rise in secondary extinctions. The dramatic increase is, of course, also a function of the fact that the other producer guilds and their primary consumers are being perturbed. One of the next steps will be to repeat these simulations, but to perturb the herbivores only. And the really next big and tedious step is to work out the analytical predictions of the topological scenario.

Three paths to network collapse

There are three responses to bottom-up perturbation that have been observed in both the models and simulations of real communities. First, there is the “typical” CEG response (described in an earlier post). The low levels of secondary extinction at low levels of perturbation are basically as expected and predicted theoretically by the topological collapse of the network as the perturbation propagates through the network.

The second response is also part of the typical CEG response, and is the nearly discontinuous and significant increase in secondary extinction, implicating a threshold response. It isn’t clear at this point exactly what precipitates the threshold change, that is, is it the loss of highly linked species? Or is it threshold responses in compensatory increases of link strength? Or something else? But, at the very least I know that it does not involve topological effects, because the threshold is neither predicted mathematically, nor is it present in the results of topological extinction-only simulations. It is therefore definitely a feature of the demographic (Lotka-Volterra-type) interspecific interactions.

Finally, there is the occasionally observed high level of secondary extinction at low levels of perturbation. Surprisingly, I think that this is not a demographic effect, as we’ve long hypothesized. It is not likely to be, since it occurs in the topological extinction simulations of the 3-guild model! Those results point to the presence of at least two distinctly different topological classes of species-level networks. One class responds according to the mathematical predictions, and the other exhibits the anomalous response. I believe that this response is essentially a result of constructing networks that have very low linkages in the critical places. The probability of this occurring must lie in the multinomial probabilities of species linkages based on guild trophic link distributions, and the uncertainty of these probabilities, or the entropy of the metanetwork, increases under certain combinations of guild species richnesses.

Nonlinear cascades

All the simulations described so far are of bottom-up perturbations to the basal level producer guilds. Network theory predicts that given networks such as food webs, with power law-like link distributions, the networks should be robust against random removal (extinction) of nodes, while being highly vulnerable to the removal (perhaps targeted) of highly linked (hub) nodes. This of course is a topological prediction, since it in no way incorporates dynamics of link strengths, compensatory modification of link strengths, or extinction thresholds (e.g. those Allee effects). Then, why do the CEG predictions and simulations of topological effects have a gradually, mildly exponential, rate of increase of secondary extinctions as the number of nodes removed is increased? The answer is two-part:

1. The probability of link loss increases with the in-degree of the consumer. Therefore most of the links being lost at any given level of perturbation are lost by highly linked species. But those species are also the most resistant to extinction.
2. The probability of secondary extinction increases almost linearly for consumers of very low degree, but almost not at all for the most highly linked species, until levels of perturbation are very high. Therefore most of the extinctions that occur at low to mid- levels of perturbation are of poorly connected species.

The nonlinear increase seen in the CEG simulations is therefore likely a response to a threshold being reached where highly linked consumers, though still robust to topological extinction, initiate significantly devastating top-down cascades because of compensatory increases of link or interaction strengths.

It is also therefore reasonable to hypothesize that very high levels of secondary extinction at low perturbation levels is the result of having a few highly connected, upper-level consumers. This could explain the great difference, at those perturbation levels, between the topological expectations and the simulations. Should low diversity communities, or communities with low diversities of high trophic level consumers, then be limited to those consumers being very specialized?

Allee effects

This is a very nice and detailed explanation of Allee Effects, a summary term for (mostly) stochastic events that heighten the probability of population extinction when population levels become very low. These effects are included in the CEG simulations as the MVP parameter, or Minimum Viable Population. MVP is taken as a proportion of the initial K (carrying capacity), and if realizable K falls to or below this level, the population is considered to be functionally extinction.

Powered by ScribeFire.