- Email: [email protected]
Ecosystems as evolutionary complex systems: Network analysis of fitness models
Environmental Modelling & Software 22 (2007) 693e700 www.elsevier.com/locate/envsoft
Ecosystems as evolutionary complex systems: Network analysis of fitness models Brian D. Fath a,*, W.E. Grant b b
a Biology Department, Towson University, Towson, MD 21252, USA Department of Wildlife and Fisheries Sciences, Texas A&M University, College Station, TX 77843, USA
Received 25 September 2005; received in revised form 24 October 2005; accepted 15 December 2005 Available online 22 March 2006
Abstract Understanding and managing ecosystems as biocomplex wholes is the compelling scientific challenge of our times. Several different systemtheoretic approaches have been proposed to study biocomplexity and two in particular, Kauffman’s NK networks and Patten’s ecological network analysis, have shown promising results. This research investigates the similarities between these two approaches, which to date have developed separately and independently. Kauffman (1993) has demonstrated that networks of non-equilibrium, open thermodynamic systems can exhibit profound order (subcritical complexity) or profound chaos (fundamental complexity). He uses Boolean NK networks to describe system behavior, where N is the number of nodes in the network and K the number of connections at each node. Ecological network analysis uses a different Boolean network approach in that the pair-wise node interactions in an ecosystem food web are scaled by the throughflow (or storage) to determine the probability of flow along each pathway in the web. These flow probabilities are used to determine systemwide properties of ecosystems such as cycling index, indirect-to-direct effects ratio, and synergism. Here we modify the NK model slightly to develop a fitness landscape of interacting species and calculate how the network analysis properties change as the model’s species coevolve. We find that, of the parameters considered, network synergism increases modestly during the simulation whereas the other properties generally decrease. Furthermore, we calculate several ecosystem level goal functions and compare their progression during increasing fitness and determine that at least at this stage there is not a good correspondence between the reductionistic and holistic drivers for the system. This research is largely a proof of concept test and will lay the foundation for future integration and model scenario analysis between two important network techniques. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Boolean networks; Coevolution; Ecological modeling; Fitness landscapes; Network analysis
1. Introduction One goal of theoretical ecosystem ecology is to identify and quantify system-level concepts and find general patterns of ecosystem organization. One promising method has been to conceptualize ecosystems as networks connected by their transfer and exchange of energy and matter within and across system boundaries. Several different developments of this
* Corresponding author. Fax: þ1 410 704 2405. E-mail address: [email protected] (B.D. Fath). 1364-8152/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2005.12.023
conceptualization have been realized. Independently, they have added significantly to our understanding of ecosystems yet there has been a lack of integration with these methods because of the different terminology, notation, history, disciplinary genesis, emphasis, and application. The main goal of this project is to find linkages between two commonly used Boolean representations of ecological networks. In particular, we link ecosystem theory based on network analysis to Kauffman’s theory of self-organized systems in order to test the hypothesis that network properties of homogenization, amplification, indirect effects, and synergism increase as an ecosystem coevolves to higher fitness levels.
694
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
2. Background 2.1. Ecological network analysis Ecological network analysis is a mathematical methodology to represent material, energy and information within an ecosystem as a network of nodes (compartments, components, storages, etc.) and connections (links, flows, etc.). Network analysis provides a system-oriented perspective to identify and quantify direct and indirect relations between all the objects in the system. After the first application of inputeoutput analysis to investigate flow distribution in ecosystems (Hannon, 1973), several formulizations of ecological network analysis have arisen including Ascendency Analysis (Ulanowicz, 1980, 1986, 1997) and Network Environ Analysis (Patten, 1978, 1981, 1982, 1985). Bernard Patten used mathematical systems theory as a foundation for studying ecosystems (Patten et al., 1976; Patten, 1978, 1981, 1982). He stressed the utility of the inclusione exclusion principle of set theory as a way to formalize the transactions that naturally occur in food webs. A binary interaction exists in ecological networks, simplified often as a question of ‘‘who eats whom,’’ but more broadly as the transfer of conservative energyematter between any two entities in the system. Much of the subsequent work in network ecology builds on this basic premise of direct energyematter transactions between coupled binary pairs. These transactions form the basis of both direct and indirect ecological relations, such as predation (direct), neutralism (direct), altruism (direct), mutualism (indirect) and competition (indirect) that are of importance to community ecology. Some of the primary findings of this research include the importance of indirect effects as they propagate through the myriad of network connections (Higashi and Patten, 1989) and synergism, individual compartments in an ecosystem gaining positive value from being embedded in a larger network (Patten, 1992; Fath and Patten, 1998). 2.2. Ecological network properties Several network properties have been developed with four in particular: amplification, indirect effects, homogenization, and synergism, used most regularly to investigate ecosystem behavior. Since they have been described elsewhere, only a brief description is provided here (see Fath and Patten, 1999a for details). The four properties relate the distribution and contribution of conservative energyematter flow through the network’s many direct and indirect pathways. One measure of resource distribution is given in the direct flow intensity, or transfer efficiency, matrix G, whose values, gij ¼ fij/Tj, represent the likelihood of flow along a given path, where fij corresponds to the flow from compartment j to compartment i, and Tj ¼ Sj(si) ¼ 0,n fij is the total sum of flow through compartment j including input and output boundary flows ( fi0 and out f0j, respectively). Tin at steady state. In the direct flow j ¼ Tj intensity matrix, G, all elements have a non-negative value less than one (0 gij < 1) and can be interpreted as a probability of flow along each pathway. Using standard inputeoutput
analysis techniques, an integral flow intensity matrix, N, is computed from the convergent power series: N ¼ G0 þ G1 þ G2 þ G3 þ . ¼ ðI GÞ
1
ð1Þ
where I is the multiplicative identity matrix. Elements nij in the matrix N include the contribution of direct (m ¼ 1) and indirect (m > 1) pathways, and therefore are always greater than or equal than the values of G. The G and N matrices are used to define the amplification, homogenization, and indirect effects properties. A specific quantitative test exists to determine each property (Table 1). Amplification occurs whenever an off-diagonal element of the integral flow matrix is greater than one (nij > 1). The integral flow from j to i, can exceed one when cycling drives more than the equivalent of one unit of input flow over the pathway. This property was observed in several of the small-scale models but is rare in large-scale models (Fath, 2004). The homogenization property compares the resource distribution between the direct and integral flow intensity matrices. It was observed that, due to the contribution of indirect pathways, flow in the integral matrix was more evenly distributed or more homogenized than that in the direct matrix, meaning that flow is comprised of contributions from many parts of the network. Network homogenization occurs when the coefficient of variation N, CV(N), is less than the coefficient of variation of G, CV(G), because this indicates that the network flow is more evenly distributed in the integral matrix (Fath and Patten 1999b). Indirect effects are calculated as the integral contributions minus the direct and initial boundary input (Indirect ¼ N I G). The indirect to direct effects ratio is a measure of the relative strength of these two factors. When the ratio is greater than one, then indirect effects are greater than direct effects. The fourth property, network synergism is based on a net flow intensity matrix, D, where dij ¼ ( fij fji)/Ti. Unlike the other series in which the elements are non-negative, entries Table 1 Four network properties
1. 2.
Property
Test
Amplification Homogenization
nij > 1 for isj CVðGÞ >1 CVðNÞ
3.
Ratio of indirect to direct effects
n P n P nij iij gij i¼1 j¼1 n P n P
>1 gij
i¼1 j¼1
4.
Synergism P þutility P >1 j utilityj
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
in D can be positive or negative (1 dij 1). The elements of D represent the relative utility between that (i,j ) pairing. An integral utility matrix U, is obtained from the power series as: U ¼ D0 þ D1 þ D2 þ D3 þ .ðI DÞ
1
ð2Þ
This methodology is used to determine qualitative relations between any two components in the network such as predation, mutualism, competition, etc. Quantitatively, synergism arises when the integral positive utility exceeds negative utility because of mutualistic relations in the system and is calculated as the ratio of the magnitude of the positive and negative utilities. 2.3. Ecological goal functions Objective functions, goal functions, optimization criteria, extremal principles, and orientors are examples of criteria that have been used to track ecosystem growth and development at a system level (e.g., Mu¨ller and Leupelt, 1998). Holistic ecological goal functions typically track certain thermodynamic or informational aspects of the ecosystem. This holistic approach to monitoring, assessing, even orienting at the ecosystem level contrasts sharply with the reductionistic approach which measures the fitness of an organism in its environment. In fact this schism between evolutionary ecology and ecosystem ecology is a major challenge to modern ecological theory (Odling-Smee et al., 2003). Odum (1969) proposed several trends to be expected in the growth and development of ecosystems from early to mature stages (Table 2). Five such goal functions, minimize specific entropy production, maximize energy throughflow, maximize exergy storage, maximize retention time, and maximize cycling, are considered here (Table 3). Table 2 Trends to be expected in ecosystem development (after Odum, 1969) Ecosystem attribute
Developmental stage
Mature stage
<1>
w1
High
Low
695
Table 3 Ecological goal functions tracked during evolution of model 1. Minimize specific entropy production (SPD) (Prigogine and Wiame, 1946). Decrease in the respiration to biomass ratio 2. Maximize energy throughflow (TST) (Odum, 1983). Increase in the internal energy flow 3. Maximize exergy storage (TSS) (Jørgensen and Mejer, 1977). Exergy storage (biomass) and information increase due to shift to more complex species composition. 4. Maximize retention time (RET) (Cheslak and Lamarra, 1981). Biological components develop mechanisms to increase time lags to maintain the energy stores longer. 5. Maximize cycling (TSC) (Morowitz, 1968). Energy and matter are retained in the system through network feedbacks.
The least specific entropy production principle, which originates from Prigogine and Wiame (1946) and Prigogine (1980, 1947), is the ratio of entropy production and the mass of system producing the entropy. In ecological terms, minimizing specific entropy means minimizing the entropy production (maintenance cost) for a given biomass (structure). In other words, an ecosystem tends to become more energetically efficient and the respiration to biomass ratio decreases over time. Maximum power principle is taken as a second thermodynamic goal function. Lotka and Odum have hypothesized that an ecosystem develops towards maximizing power (Lotka, 1922; Odum and Pinkerton, 1955), which is interpreted here as the highest possible throughflow of energy. A third perspective states that exergy storage increases during ecosystem development (Jørgensen et al., 2000; Jørgensen, 2002; Fath et al., 2001). This property, maximizing exergy storage becomes the third hypothesis considered herein. Exergy is the maximum amount of work the system can perform relative to an environmental reference state; therefore, a system with greater exergy is moved further from its reference state and further from thermodynamic equilibrium. Jørgensen and Mejer (1977, 1979) have introduced the application of exergy to account for ecosystem development by using as a reference the same system at thermodynamic equilibrium at the same temperature and pressure. In ecosystems an increase in exergy corresponds to an increase in the overall biomass and an increase in the species’ internal organizationdthis informational component is currently measured as genetic complexity (Marques et al., 1997; Jørgensen et al., 2005). Throughflow and exergy storage are related by the retention time of the system, so a fourth hypothesis, maximum retention time, is also considered here. Finally, the system cycling (Finn, 1976) described above has also been considered as an ecological goal function.(Morowitz, 1968) The goal is to find a pattern between how these various thermodynamic-based ecological goal functions (Fath et al., 2004) respond to changes in fitness in the coevolutionary fitness models described below.
Community energetics Gross production/community respiration (P/R ratio) Gross Production/standing crop biomass (P/B ratio) Biomass supported/unit energy flow (B/E ratio) Food chains
Low
High
Linear
Weblike
Nutrient cycling Mineral cycles Nutrient exchange rate
Open Rapid
Closed Slow
Undeveloped Poor Poor
Developed Good Good
2.4. Kauffman’s NK model
High Low
Low High
Stuart Kauffman uses binary Boolean networks to find general laws of system self-organization (Kauffman, 1993, 1996, 2000). His main thesis is that biological systems are composed
Overall homeostasis Internal symbiosis Nutrient conservation Stability (resistance to external perturbations) Entropy Information
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
696
of autonomous agents, or self-replicating systems that perform work, which are ‘‘co-constructing and propagating organization’’ (Kauffman, 2000, p. 5). An emphasis is placed on co-construction and coevolution because of the cybernetic feedback that makes agents adapt to other agents at the same time they modify their own environment. There recently has been renewed interest in the impact species have on each other and on their environment (e.g., Jones et al., 1997; Odling-Smee et al., 2003). Coevolution and indirect effects are both manifestations of interacting networks. In his NK model, Kauffman (2000) addresses species coevolution by coupling the influence from genes of one species to genes of another species. The basic module of the NK model represents an organism with N genes, each having two alleles, 0 and 1. The contribution of each gene to the fitness of the organism depends on the allele of that gene and the alleles of K other genes in its genome, called ‘‘epistatic’’ inputs. In this simple model there are 2N combinations of alleles that influence fitness. Each allele is randomly assigned a fitness contribution value for each of the 2N combinations in which it occurs. The overall fitness of each combination of alleles (i.e., each genotype) is taken as the mean of the fitness values of the constituent alleles. The result is a fitness landscape, such that every allele combination has a specific fitness value (Table 4 shows an example for N ¼ 3). When there is a flip in one allele from 0 to 1 or vice versa, the fitness contribution of the gene changes. If the result is higher fitness, then the allele shift is accepted, if not, then it is not accepted. Kauffman found that when the number of connections to other genes, K, is low the system quickly evolves to a global fitness maximum. As the number of connections increases there are more local peaks until the point when the system is completely interconnected (K ¼ N 1) and the resulting fitness landscape is fully random. The more local peaks that occur, then the more improbable it is to ‘‘climb’’ to the global peak, resulting, on average, in an overall lower fitness. However, Kauffman maintains that fitness landscapes are not random but instead are generated by the coevolutionary interactions of the various species. Therefore, the next step is to link NK models of various species. Table 4 Example of fitness landscape 1
2
3
Fitness value w1
Fitness value w2
Fitness value w3
Average fitness w
0 0 0 0 1 1 1 1
0 0 1 1 0 0 1 1
0 1 0 1 0 1 0 1
0.2 0.7 0.5 0.3 0.5 0.1 0.9 0.6
0.5 0.1 0.9 0.3 0.4 0.5 0.2 0.8
0.4 0.2 0.8 0.1 0.4 0.3 0.8 0.4
0.37 0.33 0.73 0.23 0.43 0.30 0.63 0.60
There are eight possible binary combinations of three genes. Each is assigned a random fitness value between 0 and 1, and the fitness for each allele combination is the mean of the three values. This procedure is used to construct a fitness landscape. For example, starting with each gene expressing a 0, the fitness is 0.37. If the allele on the first gene flips to ‘‘1’’ then fitness increases to 0.43. This simple model has only one fitness peak at (0,1,0).
In the multi-species version of Kauffman’s NK model, the fitness value of each allele depends not only on the allele of that gene and on the alleles of K epistatic genes, but also on the alleles of C other genes in each of S other species. If there are two species coupled together, then each gene has K þ C inputs, and a table of random fitness contributions is generated that has 2(KþC ) combinations. A model in which each species is connected with S other species has 2(KþCS ) possible states, so the number of possible states grows combinatorically. The fitness of the species is calculated as the mean of the fitness values of the alleles in its current genotype; each species is assumed to be isogenic. Now, when one species evolves (a flipping of an allele on a gene) this likely has ramifications for the other species by deforming the overall fitness landscape. Kauffman found that in general coevolving systems coupled in this manner behave either in an ordered or chaotic regime, separated by a phase transition depending on the number of couplings. We have recreated Kauffman’s multi-species NK model here to investigate the fitness of coevolving species with a particular interest in understanding how coevolutionary processes affect ecosystem properties. A few modifications to the original model as presented above are noted. Each time step during the simulation, any one of four events, randomly chosen, may happen. (1) A randomly chosen species may evolve to a new genotype via recombination, if the randomly chosen new genotype has a higher fitness value than the current genotype. A randomly chosen species may be replaced by a new species that (2) may have a different K than the current species, but has the same C and S, (3) may have a different C, but has the same K and S, or (4) may have a different S, but has the same K and C, if the new species has a higher fitness value than the current species. Thus, as species evolve or are replaced by invading species, they change their own fitness landscape (Kauffman, 1996, 2000) as well as the fitness landscape of the other species. The above restrictions could be relaxed in future research to study more general cases, but for now the model was used to generate a time series of connectance matrices. We apply ecological network analysis to each matrix. Eventually, it would be useful to look at models that have more realistic ecosystem structures by using the methodologies developed in Fath (2004) or perhaps to see if over time species in the models naturally evolve into a configuration similar to a trophic structure. However, that is beyond the scope of this paper. Here we present the initial results from this research, which uses a five-species model coevolving for 10 time steps under two different species coupling regimes. 3. Integrated model In the first simulation, all species were initialized with S ¼ 1 (i.e., each species is connected with one other species), and in the second simulation all species were initialized with S ¼ 4 (connected to all other species). Every time interval during the simulation, we generate a connectance matrix based on the current fitness and S values of the set of species. Elements of the connectance matrix are equal to 0 if the fitness values of
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
Sp0
0.607
0.520
Sp4
Sp1
0.657
0.415 Sp3
Sp0
0.607
0.656
Sp4 0.657
Sp2
0.227 Sp3
Sp4 0.657
0.227 Sp3
Sp2
Sp2
0.417
0.417
(a)
(b)
(c)
Sp0
Sp1
0.227
0.454
0.227 Sp3
Sp2
Sp0
0.225
Sp4
Sp1
0.225
0.454 Sp3
Sp1
0.227
0.417
Sp0
Sp4
Sp0
0.427
Sp1
0.227
697
0.225
0.138
Sp4 0.454
Sp2
0.138
0.232 Sp1
0.232
Sp3
0.138 Sp2
0.417
0.417
0.232
(d)
(e)
(f)
Fig. 1. Model evolution during ten time steps. New figures are given when adaptations are selected: (a) initial configuration and weighting factors. Configuration at (b) t ¼ 3, (c) t ¼ 5, (d) t ¼ 6, (e) t ¼ 9, (f) t ¼ 10.
the genes of the ‘‘to’’ species are not affected by the genes of the ‘‘from’’ species, and diagonal elements are equal to 0, that is, species are not connected to themselves. Values of the other elements of the connectance matrix are calculated as the fitness value of the ‘‘to’’ species divided by its S value, that is, the sum of all elements ‘‘to’’ a given species is equal to its fitness value. The elements of the connectance matrix represent the fitness contribution among connected species. In order to apply ecological network analysis to these matrices, we assume that elements of the connectance matrix represent relative rates of energy flow among the set of species. Obviously, fitness is not flow, but in a more general sense the fitness represents a measure of influence between species. The flow probability between two compartments is the proportion of flow to total throughflow ( gij ¼ fij/Tj) where Tj is the total throughflow into compartment j. This could also be interpreted as the probability of influence between two compartments (Patten et al., 1976). Here we assume that the fitness contribution (from 0 to 1) can be used as a measure of the weighted influence. This allows us to apply ecological network analysis to each matrix and calculate the four ecological network properties described above. We then examined the temporal dynamics of these properties as the set of species co-evolve through different fitness landscapes to test the hypothesis that cycling index, homogenization, amplification, indirect effects, and synergism increase as the ecosystem coevolves. Note that ecological network analysis is a steady-state analysis; however, we treat the model generated from each time step as a snapshot in time. As the system changes over time, we can determine the network properties of the system in that particular state. One other assumption is needed to run the analysis, which is that the model ecosystems, as open systems, receive external input. Energy enters the system largely through primary
producer and lower trophic level species. Usually, for a model this size (5 compartments) external input into one compartment is enough, but in some of these simulations the first compartment is eliminated after which time there would be no further input available to higher trophic levels. Therefore, a unit of input is given to each of the five compartments.
4. Numerical simulation results In the first simulation, each species is connected to one other species. The connectance values can change at each time step given the occurrence of a randomly chosen event, as described above. Fig. 1 shows the digraph changes for the first simulation in which all species are originally connected Table 5 Example of two connectance matrices generated by the first simulation model Sp 1
Sp 2
Sp 3
Sp 4
Sp 5
T¼2 Sp 1 Sp 2 Sp 3 Sp 4 Sp 5
0 0 0 0 0.61
0.52 0 0 0 0
0 0.42 0 0 0
0 0 0.42 0 0
0 0 0 0.66 0
T¼3 Sp 1 Sp 2 Sp 3 Sp 4 Sp 5
0 0 0 0 0.61
0.66 0 0 0 0
0 0.23 0 0 0
0 0.23 0.42 0 0
0 0 0 0.66 0
Reading from columns to rows, at time 2, Sp 2 affects Sp 1 (0.52), Sp 3 affects Sp 2 (0.42), Sp 4 affects Sp 3 (0.42), Sp 5 affects Sp 4 (0.66), and Sp 1 affects Sp 5 (0.61). At T ¼ 3, Sp 2 is replaced by a new Sp 2 that is affected by Sp 3 and Sp 4 (overall fitness is higher (0.46 versus 0.42). The new Sp 2 also has caused a change in the fitness value of Sp 1.
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
698
20 18 16 14 12 10 8 6 4 2 0
2.5
5
parameter value
parameter value
6
4 3 2 1
2 1.5 1 0.5
0 0
2
4
6
8
10
0
time step Cycling index
homogenization
0 synergism
2
4
6
8
10
time step
i/d ratio Cycling index
Fig. 2. Behavior of network properties for first simulation.
homogenization
i/d ratio
synergism
Fig. 3. Behavior of network properties for second simulation. Synergism is plotted on the alternate y-axis.
Table 6 Connections in Simulation A (species initially connected to one species) and Simulation B (species initially connected to four species) T
A: # links
B: # links
0 1 2 3 4 5 6 7 8 9 10
5 5 5 6 6 5 4 4 4 5 7
20 19 19 19 19 19 19 18 18 18 14
simulation and dropping again near the end (Fig. 3; note in the figure that synergism is plotted on the alternate y-axis). Next we observed the progression of certain goal functions to see how they also responded to changes in the fitness model (Figs. 4 and 5). In both simulations, the total system throughflow and total system storage increase. Retention time increases initially and then decreases and the overall trend for specific dissipation is to increase. Cycling, which was plotted on Figs. 2 and 3 as one of the network properties, decreases during the simulations. Based on these results only the maximize exergy storage and maximize power goal functions behaved as expected; however, the models used in this proof of concept stage were not ecologically realistic. Future research will involve further refinement of the model by considering other initial structural patterns, using different sized models, and allowing the model to run longer. 5. Discussion One theme of the International Environmental Modelling and Software Society (iEMSs) Conference, which the papers in this issue represent, is how to acquire through science and modeling a better understanding of complexity in environmental systems, primarily the sustainability of coupled ecological and social systems. One of the main motivations for Kauffman’s Boolean networks was to identify fundamental 35
property value
to one other species. Table 5 shows two matrices generated by the model at time steps 2 and 3. For example, we see that in the third time step a new species 2 appears which is also dependent on species 4 and the overall connectance or fitness from species 2 to species 1 increases. Changes such as these continue through to the end of the simulation after 10 time steps. When the ecological network properties of these connectance matrices from each time step are calculated we find the following: amplification does not occur at any time step; the cycling index, homogenization, and ratio of indirect to direct effects all decrease over time; and the synergism parameter rises until T ¼ 9 at which time it starts to drop by the end of the simulation (Fig. 2). In the second simulation, each species is initially linked to four other species (complete graph minus self-loops). Several changes occur immediately, most notably, the connection between Sp 4 and Sp 5 is lost. During the 10-step simulation the system becomes more articulated, meaning there are fewer connections between species, but these changes would only be accepted if the overall fitness of the species increases. One simple measure to consider is the total number of connections in the system during each time step (Table 6). We see a similar pattern in the network parameters in the second simulation as well. Amplification does not occur at any time step. Cycling index and indirect effects ratio decrease, while in this simulation homogenization bounces around but is fairly flat. Synergism also oscillates reaching a peak in the middle of the
0.1
30 25 20
0.05
15 10 5 0
0 0
2
4
6
8
10
time step TSS
RET
TST
SDP
Fig. 4. Behavior of ecological goal functions for first simulation. TSS and RET are plotted on the alternate y-axis.
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
property value
35
0.1
30 25 20
0.05
15 10 5 0
0 0
2
4
6
8
RET
TST
of influence between the compartments. This allows the application of network analysis techniques to determine the values of specific network properties. In particular, we found that network synergism appears to respond positively as fitness increases, and the other properties respond negatively. This paper represents the first attempt to integrate the two Boolean techniques; further research is needed to more deeply understand the interrelation between them.
10
time step TSS
699
References SDP
Fig. 5. Behavior of ecological goal functions for second simulation. TSS and RET are plotted on the alternate y-axis.
principles relating network structure and behavior. He found that system behavior ranged from very stable (‘‘deeply ordered’’ in Kauffman’s terminology) to chaotic depending on the number of nodes and connections, and that systems tended to evolve toward an intermediate level of connectivity (‘‘the edge of chaos’’). These results indicate that there is an optimal complexity for system design, which balances the opposing evolutionary tendencies toward increased stability (and hence less adaptability) and toward increased flexibility (with the concomitant risk of losing system integrity). In the ecological examples given above, we can see this tendency for systems to evolve toward an intermediate level of connectivity (Table 6). The network models described herein, although theoretical in nature, provide insight into system organization and due to their reliance on first principles of Boolean analysis are applicable to many different system types. Therefore, they provide a good platform for modeling integrated systems. This paper presents an initial foray into two different important network approaches with the aim of gaining better confidence in how they may be linked or they complement each other. The ecological network properties inform about structure, function, and relations, but they have traditionally been applied to static networks. Kauffman’s NK model gives a way to let the models evolve over time, and thus see how the properties change over time. Since the ecological goal functions are assumed to track ecosystem development, they have already been used for structural dynamic models (Nielsen, 1992) and model calibration (Jørgensen and Fath, 2004). It is also important to find other useful management applications for these methodologies. In the case of managing biological and societal systems interacting together and coevolving on correlated landscapes, a hierarchical management approach dealing with the subsystems shows potential to find a compromise solution between complexity levels of the different subsystems (Alexander, 1964).
6. Conclusions In conclusion, we have recreated Kauffman’s multi-species NK model and used it to investigate the coevolution of a simple model ecosystem. Furthermore, we have used the fitness values generated by the model as surrogates for the probability
Alexander, C., 1964. Notes on the Synthesis of Form. Harvard University Press, Cambridge. Cheslak, E.F., Lamarra, V.A., 1981. The residence time of energy as a measure of ecological organization. In: Mitsch, W.J., Bossermann, R.W., Klopatek, J.M. (Eds.), Energy and Ecological Modelling. Elsevier, Amsterdam, pp. 591e600. Fath, B.D., 2004. Network analysis applied to large-scale cyber-ecosystems. Ecol. Modell 171, 329e337. Fath, B.D., Patten, B.C., 1998. Network synergism: emergence of positive relations in ecological systems. Ecol. Model 107, 127e143. Fath, B.D., Patten, B.C., 1999a. Review of the foundations of network environ analysis. Ecosystems 2, 167e179. Fath, B.D., Patten, B.C., 1999b. Quantification of resource homogenization using network flow analysis. Ecol. Model 123, 193e205. Fath, B.D., Patten, B.C., Choi, J.S., 2001. Complementarity of ecological goal functions. J. Theor. Biol. 208 (4), 493e506. Fath, B.D., Jørgensen, S.E., Patten, B.C., Strasˇkraba, M., 2004. Ecosystem growth and development. Biosystems 77, 213e228. Finn, J.T., 1976. Measures of ecosystem structure and function derived from flow analysis. J. Theor. Biol. 56, 363e380. Hannon, B., 1973. The structure of ecosystems. J Theor. Biol. 41, 535e546. Higashi, M., Patten, B.C., 1989. Dominance of indirect causality in ecosystems. Amer. Nat 133, 288e302. Jones, C.G., Lawton, J.H., Shachak, M., 1997. Positive and negative effects of organisms as physical ecosystem engineers. Ecology 78, 1946e1957. Jørgensen, S.E., Fath, B.D., 2004. Modelling the selective adaptation of Darwin’s Finches. Ecol. Model. 176, 409e418. Jørgensen, S.E., Mejer, H.F., 1977. Ecological buffer capacity. Ecol. Model 3, 39e61. Jørgensen, S.E., Mejer, H.F., 1979. A holistic approach to ecological modelling. Ecol. Model 7, 169e189. Jørgensen, S.E., Mejer, S.E., 2002. Integration of Ecosystem Theories: A Pattern, third ed. Kluwer Academic Publishers , Dordrecht. Jørgensen, S.E., Patten, B.C., Strasˇkraba, M., 2000. Ecosystems emerging: 4. Growth. Ecol. Model 126, 249e284. Jørgensen, S.J., Ladegaard, N., Debeljak, M., Marques, J.C., 2005. Calculations of exergy for organisms. Ecol. Model 185, 165e175. Kauffman, S.A., 1993. The Origins of Order: Self-Organization and Selection in Evolution. Oxford University Press, New York. Kauffman, S.A., 1996. At Home in the Universe: The Search for Laws of SelfOrganization and Complexity. Oxford University Press, New York. Kauffman, S.A., 2000. Investigations. Oxford University Press, New York. Lotka, A.J., 1922. Contribution to the energetics of evolution. Proc. Nat. Acad. Science 8, 147e151. Marques, J.C., Pardal, B.C., Nielsen, S.N., Jørgensen, S.E., 1997. Analysis of the properties of exergy and biodiversity along an estuarine gradient of eutrophication. Ecol. Model 102, 155e167. Morowitz, H.J., 1968. Energy Flow in Biology: Biological Organization as a Problem in Thermal Physics. Academic Press, New York. Mu¨ller, F., Leupelt, M. (Eds.), 1998. Eco Targets, Goal Functions, and Orientors. Springer-Verlag, New York, p. 619. Nielsen, S.N., 1992. Strategies for structural-dynamical modelling. Ecol. Model 63, 91e102. Odling-Smee, F.J., Laland, K.N., Feldman, M.W., 2003. Niche Construction: The Neglected Process in Evolution. Princeton University Press, Princeton, 472 pp.
700
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
Odum, E.P., 1969. The strategy of ecosystem development. Science 164, 262e 270. Odum, H.T., 1983. Systems Ecology: An Introduction. John Wiley, New York. Odum, H.T., Pinkerton, R.C., 1955. Time’s speed regulator: the optimum efficiency for maximum power output in physical and biological systems. Am. Sci. 43, 331e343. Patten, B.C., 1978. Systems approach to the concept of environment. Ohio J. Sci. 78, 206e222. Patten, B.C., 1981. Environs: the superniches of ecosystems. Amer. Zoo. 21, 845e852. Patten, B.C., 1982. Environs: relativistic elementary particles or ecology. Amer. Nat. 119, 179e219. Patten, B.C., 1985. Energy cycling in the ecosystem. Ecol. Model. 28, 1e71.
Patten, B.C., 1992. Energy, emergy and environs. Ecol. Modell 62, 29e69. Patten, B.C., Bosserman, R.W., Finn, J.T., Cale, W.G., 1976. Propagation of cause in ecosystems. In: Patten, B.C. (Ed.), Systems Analysis and Simulation in Ecology, vol. IV. Academic Press, New York, pp. 457e579. 593. Prigogine, I., 1980. From Being to Becoming: Time and Complexity in the Physical Sciences. WH Freeman and Co, San Francisco. Prigogine, I., Wiame, J.M., 1946. Biologie et Thermodynamique des Phenomenes Irreversible. Experientia II, 451e453. Ulanowicz, R.E., 1980. An hypothesis on the development of natural communities. J. Theor. Biol. 85, 223e245. Ulanowicz, R.E., 1986. Growth and Development: Ecosystem Phenomenology. Springer-Verlag, New York. Ulanowicz, R.E., 1997. Ecology, The Ascendent Perspective. Columbia University Press, New York.
Ecosystems as evolutionary complex systems: Network analysis of fitness models Brian D. Fath a,*, W.E. Grant b b
a Biology Department, Towson University, Towson, MD 21252, USA Department of Wildlife and Fisheries Sciences, Texas A&M University, College Station, TX 77843, USA
Received 25 September 2005; received in revised form 24 October 2005; accepted 15 December 2005 Available online 22 March 2006
Abstract Understanding and managing ecosystems as biocomplex wholes is the compelling scientific challenge of our times. Several different systemtheoretic approaches have been proposed to study biocomplexity and two in particular, Kauffman’s NK networks and Patten’s ecological network analysis, have shown promising results. This research investigates the similarities between these two approaches, which to date have developed separately and independently. Kauffman (1993) has demonstrated that networks of non-equilibrium, open thermodynamic systems can exhibit profound order (subcritical complexity) or profound chaos (fundamental complexity). He uses Boolean NK networks to describe system behavior, where N is the number of nodes in the network and K the number of connections at each node. Ecological network analysis uses a different Boolean network approach in that the pair-wise node interactions in an ecosystem food web are scaled by the throughflow (or storage) to determine the probability of flow along each pathway in the web. These flow probabilities are used to determine systemwide properties of ecosystems such as cycling index, indirect-to-direct effects ratio, and synergism. Here we modify the NK model slightly to develop a fitness landscape of interacting species and calculate how the network analysis properties change as the model’s species coevolve. We find that, of the parameters considered, network synergism increases modestly during the simulation whereas the other properties generally decrease. Furthermore, we calculate several ecosystem level goal functions and compare their progression during increasing fitness and determine that at least at this stage there is not a good correspondence between the reductionistic and holistic drivers for the system. This research is largely a proof of concept test and will lay the foundation for future integration and model scenario analysis between two important network techniques. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Boolean networks; Coevolution; Ecological modeling; Fitness landscapes; Network analysis
1. Introduction One goal of theoretical ecosystem ecology is to identify and quantify system-level concepts and find general patterns of ecosystem organization. One promising method has been to conceptualize ecosystems as networks connected by their transfer and exchange of energy and matter within and across system boundaries. Several different developments of this
* Corresponding author. Fax: þ1 410 704 2405. E-mail address: [email protected] (B.D. Fath). 1364-8152/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2005.12.023
conceptualization have been realized. Independently, they have added significantly to our understanding of ecosystems yet there has been a lack of integration with these methods because of the different terminology, notation, history, disciplinary genesis, emphasis, and application. The main goal of this project is to find linkages between two commonly used Boolean representations of ecological networks. In particular, we link ecosystem theory based on network analysis to Kauffman’s theory of self-organized systems in order to test the hypothesis that network properties of homogenization, amplification, indirect effects, and synergism increase as an ecosystem coevolves to higher fitness levels.
694
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
2. Background 2.1. Ecological network analysis Ecological network analysis is a mathematical methodology to represent material, energy and information within an ecosystem as a network of nodes (compartments, components, storages, etc.) and connections (links, flows, etc.). Network analysis provides a system-oriented perspective to identify and quantify direct and indirect relations between all the objects in the system. After the first application of inputeoutput analysis to investigate flow distribution in ecosystems (Hannon, 1973), several formulizations of ecological network analysis have arisen including Ascendency Analysis (Ulanowicz, 1980, 1986, 1997) and Network Environ Analysis (Patten, 1978, 1981, 1982, 1985). Bernard Patten used mathematical systems theory as a foundation for studying ecosystems (Patten et al., 1976; Patten, 1978, 1981, 1982). He stressed the utility of the inclusione exclusion principle of set theory as a way to formalize the transactions that naturally occur in food webs. A binary interaction exists in ecological networks, simplified often as a question of ‘‘who eats whom,’’ but more broadly as the transfer of conservative energyematter between any two entities in the system. Much of the subsequent work in network ecology builds on this basic premise of direct energyematter transactions between coupled binary pairs. These transactions form the basis of both direct and indirect ecological relations, such as predation (direct), neutralism (direct), altruism (direct), mutualism (indirect) and competition (indirect) that are of importance to community ecology. Some of the primary findings of this research include the importance of indirect effects as they propagate through the myriad of network connections (Higashi and Patten, 1989) and synergism, individual compartments in an ecosystem gaining positive value from being embedded in a larger network (Patten, 1992; Fath and Patten, 1998). 2.2. Ecological network properties Several network properties have been developed with four in particular: amplification, indirect effects, homogenization, and synergism, used most regularly to investigate ecosystem behavior. Since they have been described elsewhere, only a brief description is provided here (see Fath and Patten, 1999a for details). The four properties relate the distribution and contribution of conservative energyematter flow through the network’s many direct and indirect pathways. One measure of resource distribution is given in the direct flow intensity, or transfer efficiency, matrix G, whose values, gij ¼ fij/Tj, represent the likelihood of flow along a given path, where fij corresponds to the flow from compartment j to compartment i, and Tj ¼ Sj(si) ¼ 0,n fij is the total sum of flow through compartment j including input and output boundary flows ( fi0 and out f0j, respectively). Tin at steady state. In the direct flow j ¼ Tj intensity matrix, G, all elements have a non-negative value less than one (0 gij < 1) and can be interpreted as a probability of flow along each pathway. Using standard inputeoutput
analysis techniques, an integral flow intensity matrix, N, is computed from the convergent power series: N ¼ G0 þ G1 þ G2 þ G3 þ . ¼ ðI GÞ
1
ð1Þ
where I is the multiplicative identity matrix. Elements nij in the matrix N include the contribution of direct (m ¼ 1) and indirect (m > 1) pathways, and therefore are always greater than or equal than the values of G. The G and N matrices are used to define the amplification, homogenization, and indirect effects properties. A specific quantitative test exists to determine each property (Table 1). Amplification occurs whenever an off-diagonal element of the integral flow matrix is greater than one (nij > 1). The integral flow from j to i, can exceed one when cycling drives more than the equivalent of one unit of input flow over the pathway. This property was observed in several of the small-scale models but is rare in large-scale models (Fath, 2004). The homogenization property compares the resource distribution between the direct and integral flow intensity matrices. It was observed that, due to the contribution of indirect pathways, flow in the integral matrix was more evenly distributed or more homogenized than that in the direct matrix, meaning that flow is comprised of contributions from many parts of the network. Network homogenization occurs when the coefficient of variation N, CV(N), is less than the coefficient of variation of G, CV(G), because this indicates that the network flow is more evenly distributed in the integral matrix (Fath and Patten 1999b). Indirect effects are calculated as the integral contributions minus the direct and initial boundary input (Indirect ¼ N I G). The indirect to direct effects ratio is a measure of the relative strength of these two factors. When the ratio is greater than one, then indirect effects are greater than direct effects. The fourth property, network synergism is based on a net flow intensity matrix, D, where dij ¼ ( fij fji)/Ti. Unlike the other series in which the elements are non-negative, entries Table 1 Four network properties
1. 2.
Property
Test
Amplification Homogenization
nij > 1 for isj CVðGÞ >1 CVðNÞ
3.
Ratio of indirect to direct effects
n P n P nij iij gij i¼1 j¼1 n P n P
>1 gij
i¼1 j¼1
4.
Synergism P þutility P >1 j utilityj
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
in D can be positive or negative (1 dij 1). The elements of D represent the relative utility between that (i,j ) pairing. An integral utility matrix U, is obtained from the power series as: U ¼ D0 þ D1 þ D2 þ D3 þ .ðI DÞ
1
ð2Þ
This methodology is used to determine qualitative relations between any two components in the network such as predation, mutualism, competition, etc. Quantitatively, synergism arises when the integral positive utility exceeds negative utility because of mutualistic relations in the system and is calculated as the ratio of the magnitude of the positive and negative utilities. 2.3. Ecological goal functions Objective functions, goal functions, optimization criteria, extremal principles, and orientors are examples of criteria that have been used to track ecosystem growth and development at a system level (e.g., Mu¨ller and Leupelt, 1998). Holistic ecological goal functions typically track certain thermodynamic or informational aspects of the ecosystem. This holistic approach to monitoring, assessing, even orienting at the ecosystem level contrasts sharply with the reductionistic approach which measures the fitness of an organism in its environment. In fact this schism between evolutionary ecology and ecosystem ecology is a major challenge to modern ecological theory (Odling-Smee et al., 2003). Odum (1969) proposed several trends to be expected in the growth and development of ecosystems from early to mature stages (Table 2). Five such goal functions, minimize specific entropy production, maximize energy throughflow, maximize exergy storage, maximize retention time, and maximize cycling, are considered here (Table 3). Table 2 Trends to be expected in ecosystem development (after Odum, 1969) Ecosystem attribute
Developmental stage
Mature stage
<1>
w1
High
Low
695
Table 3 Ecological goal functions tracked during evolution of model 1. Minimize specific entropy production (SPD) (Prigogine and Wiame, 1946). Decrease in the respiration to biomass ratio 2. Maximize energy throughflow (TST) (Odum, 1983). Increase in the internal energy flow 3. Maximize exergy storage (TSS) (Jørgensen and Mejer, 1977). Exergy storage (biomass) and information increase due to shift to more complex species composition. 4. Maximize retention time (RET) (Cheslak and Lamarra, 1981). Biological components develop mechanisms to increase time lags to maintain the energy stores longer. 5. Maximize cycling (TSC) (Morowitz, 1968). Energy and matter are retained in the system through network feedbacks.
The least specific entropy production principle, which originates from Prigogine and Wiame (1946) and Prigogine (1980, 1947), is the ratio of entropy production and the mass of system producing the entropy. In ecological terms, minimizing specific entropy means minimizing the entropy production (maintenance cost) for a given biomass (structure). In other words, an ecosystem tends to become more energetically efficient and the respiration to biomass ratio decreases over time. Maximum power principle is taken as a second thermodynamic goal function. Lotka and Odum have hypothesized that an ecosystem develops towards maximizing power (Lotka, 1922; Odum and Pinkerton, 1955), which is interpreted here as the highest possible throughflow of energy. A third perspective states that exergy storage increases during ecosystem development (Jørgensen et al., 2000; Jørgensen, 2002; Fath et al., 2001). This property, maximizing exergy storage becomes the third hypothesis considered herein. Exergy is the maximum amount of work the system can perform relative to an environmental reference state; therefore, a system with greater exergy is moved further from its reference state and further from thermodynamic equilibrium. Jørgensen and Mejer (1977, 1979) have introduced the application of exergy to account for ecosystem development by using as a reference the same system at thermodynamic equilibrium at the same temperature and pressure. In ecosystems an increase in exergy corresponds to an increase in the overall biomass and an increase in the species’ internal organizationdthis informational component is currently measured as genetic complexity (Marques et al., 1997; Jørgensen et al., 2005). Throughflow and exergy storage are related by the retention time of the system, so a fourth hypothesis, maximum retention time, is also considered here. Finally, the system cycling (Finn, 1976) described above has also been considered as an ecological goal function.(Morowitz, 1968) The goal is to find a pattern between how these various thermodynamic-based ecological goal functions (Fath et al., 2004) respond to changes in fitness in the coevolutionary fitness models described below.
Community energetics Gross production/community respiration (P/R ratio) Gross Production/standing crop biomass (P/B ratio) Biomass supported/unit energy flow (B/E ratio) Food chains
Low
High
Linear
Weblike
Nutrient cycling Mineral cycles Nutrient exchange rate
Open Rapid
Closed Slow
Undeveloped Poor Poor
Developed Good Good
2.4. Kauffman’s NK model
High Low
Low High
Stuart Kauffman uses binary Boolean networks to find general laws of system self-organization (Kauffman, 1993, 1996, 2000). His main thesis is that biological systems are composed
Overall homeostasis Internal symbiosis Nutrient conservation Stability (resistance to external perturbations) Entropy Information
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
696
of autonomous agents, or self-replicating systems that perform work, which are ‘‘co-constructing and propagating organization’’ (Kauffman, 2000, p. 5). An emphasis is placed on co-construction and coevolution because of the cybernetic feedback that makes agents adapt to other agents at the same time they modify their own environment. There recently has been renewed interest in the impact species have on each other and on their environment (e.g., Jones et al., 1997; Odling-Smee et al., 2003). Coevolution and indirect effects are both manifestations of interacting networks. In his NK model, Kauffman (2000) addresses species coevolution by coupling the influence from genes of one species to genes of another species. The basic module of the NK model represents an organism with N genes, each having two alleles, 0 and 1. The contribution of each gene to the fitness of the organism depends on the allele of that gene and the alleles of K other genes in its genome, called ‘‘epistatic’’ inputs. In this simple model there are 2N combinations of alleles that influence fitness. Each allele is randomly assigned a fitness contribution value for each of the 2N combinations in which it occurs. The overall fitness of each combination of alleles (i.e., each genotype) is taken as the mean of the fitness values of the constituent alleles. The result is a fitness landscape, such that every allele combination has a specific fitness value (Table 4 shows an example for N ¼ 3). When there is a flip in one allele from 0 to 1 or vice versa, the fitness contribution of the gene changes. If the result is higher fitness, then the allele shift is accepted, if not, then it is not accepted. Kauffman found that when the number of connections to other genes, K, is low the system quickly evolves to a global fitness maximum. As the number of connections increases there are more local peaks until the point when the system is completely interconnected (K ¼ N 1) and the resulting fitness landscape is fully random. The more local peaks that occur, then the more improbable it is to ‘‘climb’’ to the global peak, resulting, on average, in an overall lower fitness. However, Kauffman maintains that fitness landscapes are not random but instead are generated by the coevolutionary interactions of the various species. Therefore, the next step is to link NK models of various species. Table 4 Example of fitness landscape 1
2
3
Fitness value w1
Fitness value w2
Fitness value w3
Average fitness w
0 0 0 0 1 1 1 1
0 0 1 1 0 0 1 1
0 1 0 1 0 1 0 1
0.2 0.7 0.5 0.3 0.5 0.1 0.9 0.6
0.5 0.1 0.9 0.3 0.4 0.5 0.2 0.8
0.4 0.2 0.8 0.1 0.4 0.3 0.8 0.4
0.37 0.33 0.73 0.23 0.43 0.30 0.63 0.60
There are eight possible binary combinations of three genes. Each is assigned a random fitness value between 0 and 1, and the fitness for each allele combination is the mean of the three values. This procedure is used to construct a fitness landscape. For example, starting with each gene expressing a 0, the fitness is 0.37. If the allele on the first gene flips to ‘‘1’’ then fitness increases to 0.43. This simple model has only one fitness peak at (0,1,0).
In the multi-species version of Kauffman’s NK model, the fitness value of each allele depends not only on the allele of that gene and on the alleles of K epistatic genes, but also on the alleles of C other genes in each of S other species. If there are two species coupled together, then each gene has K þ C inputs, and a table of random fitness contributions is generated that has 2(KþC ) combinations. A model in which each species is connected with S other species has 2(KþCS ) possible states, so the number of possible states grows combinatorically. The fitness of the species is calculated as the mean of the fitness values of the alleles in its current genotype; each species is assumed to be isogenic. Now, when one species evolves (a flipping of an allele on a gene) this likely has ramifications for the other species by deforming the overall fitness landscape. Kauffman found that in general coevolving systems coupled in this manner behave either in an ordered or chaotic regime, separated by a phase transition depending on the number of couplings. We have recreated Kauffman’s multi-species NK model here to investigate the fitness of coevolving species with a particular interest in understanding how coevolutionary processes affect ecosystem properties. A few modifications to the original model as presented above are noted. Each time step during the simulation, any one of four events, randomly chosen, may happen. (1) A randomly chosen species may evolve to a new genotype via recombination, if the randomly chosen new genotype has a higher fitness value than the current genotype. A randomly chosen species may be replaced by a new species that (2) may have a different K than the current species, but has the same C and S, (3) may have a different C, but has the same K and S, or (4) may have a different S, but has the same K and C, if the new species has a higher fitness value than the current species. Thus, as species evolve or are replaced by invading species, they change their own fitness landscape (Kauffman, 1996, 2000) as well as the fitness landscape of the other species. The above restrictions could be relaxed in future research to study more general cases, but for now the model was used to generate a time series of connectance matrices. We apply ecological network analysis to each matrix. Eventually, it would be useful to look at models that have more realistic ecosystem structures by using the methodologies developed in Fath (2004) or perhaps to see if over time species in the models naturally evolve into a configuration similar to a trophic structure. However, that is beyond the scope of this paper. Here we present the initial results from this research, which uses a five-species model coevolving for 10 time steps under two different species coupling regimes. 3. Integrated model In the first simulation, all species were initialized with S ¼ 1 (i.e., each species is connected with one other species), and in the second simulation all species were initialized with S ¼ 4 (connected to all other species). Every time interval during the simulation, we generate a connectance matrix based on the current fitness and S values of the set of species. Elements of the connectance matrix are equal to 0 if the fitness values of
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
Sp0
0.607
0.520
Sp4
Sp1
0.657
0.415 Sp3
Sp0
0.607
0.656
Sp4 0.657
Sp2
0.227 Sp3
Sp4 0.657
0.227 Sp3
Sp2
Sp2
0.417
0.417
(a)
(b)
(c)
Sp0
Sp1
0.227
0.454
0.227 Sp3
Sp2
Sp0
0.225
Sp4
Sp1
0.225
0.454 Sp3
Sp1
0.227
0.417
Sp0
Sp4
Sp0
0.427
Sp1
0.227
697
0.225
0.138
Sp4 0.454
Sp2
0.138
0.232 Sp1
0.232
Sp3
0.138 Sp2
0.417
0.417
0.232
(d)
(e)
(f)
Fig. 1. Model evolution during ten time steps. New figures are given when adaptations are selected: (a) initial configuration and weighting factors. Configuration at (b) t ¼ 3, (c) t ¼ 5, (d) t ¼ 6, (e) t ¼ 9, (f) t ¼ 10.
the genes of the ‘‘to’’ species are not affected by the genes of the ‘‘from’’ species, and diagonal elements are equal to 0, that is, species are not connected to themselves. Values of the other elements of the connectance matrix are calculated as the fitness value of the ‘‘to’’ species divided by its S value, that is, the sum of all elements ‘‘to’’ a given species is equal to its fitness value. The elements of the connectance matrix represent the fitness contribution among connected species. In order to apply ecological network analysis to these matrices, we assume that elements of the connectance matrix represent relative rates of energy flow among the set of species. Obviously, fitness is not flow, but in a more general sense the fitness represents a measure of influence between species. The flow probability between two compartments is the proportion of flow to total throughflow ( gij ¼ fij/Tj) where Tj is the total throughflow into compartment j. This could also be interpreted as the probability of influence between two compartments (Patten et al., 1976). Here we assume that the fitness contribution (from 0 to 1) can be used as a measure of the weighted influence. This allows us to apply ecological network analysis to each matrix and calculate the four ecological network properties described above. We then examined the temporal dynamics of these properties as the set of species co-evolve through different fitness landscapes to test the hypothesis that cycling index, homogenization, amplification, indirect effects, and synergism increase as the ecosystem coevolves. Note that ecological network analysis is a steady-state analysis; however, we treat the model generated from each time step as a snapshot in time. As the system changes over time, we can determine the network properties of the system in that particular state. One other assumption is needed to run the analysis, which is that the model ecosystems, as open systems, receive external input. Energy enters the system largely through primary
producer and lower trophic level species. Usually, for a model this size (5 compartments) external input into one compartment is enough, but in some of these simulations the first compartment is eliminated after which time there would be no further input available to higher trophic levels. Therefore, a unit of input is given to each of the five compartments.
4. Numerical simulation results In the first simulation, each species is connected to one other species. The connectance values can change at each time step given the occurrence of a randomly chosen event, as described above. Fig. 1 shows the digraph changes for the first simulation in which all species are originally connected Table 5 Example of two connectance matrices generated by the first simulation model Sp 1
Sp 2
Sp 3
Sp 4
Sp 5
T¼2 Sp 1 Sp 2 Sp 3 Sp 4 Sp 5
0 0 0 0 0.61
0.52 0 0 0 0
0 0.42 0 0 0
0 0 0.42 0 0
0 0 0 0.66 0
T¼3 Sp 1 Sp 2 Sp 3 Sp 4 Sp 5
0 0 0 0 0.61
0.66 0 0 0 0
0 0.23 0 0 0
0 0.23 0.42 0 0
0 0 0 0.66 0
Reading from columns to rows, at time 2, Sp 2 affects Sp 1 (0.52), Sp 3 affects Sp 2 (0.42), Sp 4 affects Sp 3 (0.42), Sp 5 affects Sp 4 (0.66), and Sp 1 affects Sp 5 (0.61). At T ¼ 3, Sp 2 is replaced by a new Sp 2 that is affected by Sp 3 and Sp 4 (overall fitness is higher (0.46 versus 0.42). The new Sp 2 also has caused a change in the fitness value of Sp 1.
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
698
20 18 16 14 12 10 8 6 4 2 0
2.5
5
parameter value
parameter value
6
4 3 2 1
2 1.5 1 0.5
0 0
2
4
6
8
10
0
time step Cycling index
homogenization
0 synergism
2
4
6
8
10
time step
i/d ratio Cycling index
Fig. 2. Behavior of network properties for first simulation.
homogenization
i/d ratio
synergism
Fig. 3. Behavior of network properties for second simulation. Synergism is plotted on the alternate y-axis.
Table 6 Connections in Simulation A (species initially connected to one species) and Simulation B (species initially connected to four species) T
A: # links
B: # links
0 1 2 3 4 5 6 7 8 9 10
5 5 5 6 6 5 4 4 4 5 7
20 19 19 19 19 19 19 18 18 18 14
simulation and dropping again near the end (Fig. 3; note in the figure that synergism is plotted on the alternate y-axis). Next we observed the progression of certain goal functions to see how they also responded to changes in the fitness model (Figs. 4 and 5). In both simulations, the total system throughflow and total system storage increase. Retention time increases initially and then decreases and the overall trend for specific dissipation is to increase. Cycling, which was plotted on Figs. 2 and 3 as one of the network properties, decreases during the simulations. Based on these results only the maximize exergy storage and maximize power goal functions behaved as expected; however, the models used in this proof of concept stage were not ecologically realistic. Future research will involve further refinement of the model by considering other initial structural patterns, using different sized models, and allowing the model to run longer. 5. Discussion One theme of the International Environmental Modelling and Software Society (iEMSs) Conference, which the papers in this issue represent, is how to acquire through science and modeling a better understanding of complexity in environmental systems, primarily the sustainability of coupled ecological and social systems. One of the main motivations for Kauffman’s Boolean networks was to identify fundamental 35
property value
to one other species. Table 5 shows two matrices generated by the model at time steps 2 and 3. For example, we see that in the third time step a new species 2 appears which is also dependent on species 4 and the overall connectance or fitness from species 2 to species 1 increases. Changes such as these continue through to the end of the simulation after 10 time steps. When the ecological network properties of these connectance matrices from each time step are calculated we find the following: amplification does not occur at any time step; the cycling index, homogenization, and ratio of indirect to direct effects all decrease over time; and the synergism parameter rises until T ¼ 9 at which time it starts to drop by the end of the simulation (Fig. 2). In the second simulation, each species is initially linked to four other species (complete graph minus self-loops). Several changes occur immediately, most notably, the connection between Sp 4 and Sp 5 is lost. During the 10-step simulation the system becomes more articulated, meaning there are fewer connections between species, but these changes would only be accepted if the overall fitness of the species increases. One simple measure to consider is the total number of connections in the system during each time step (Table 6). We see a similar pattern in the network parameters in the second simulation as well. Amplification does not occur at any time step. Cycling index and indirect effects ratio decrease, while in this simulation homogenization bounces around but is fairly flat. Synergism also oscillates reaching a peak in the middle of the
0.1
30 25 20
0.05
15 10 5 0
0 0
2
4
6
8
10
time step TSS
RET
TST
SDP
Fig. 4. Behavior of ecological goal functions for first simulation. TSS and RET are plotted on the alternate y-axis.
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
property value
35
0.1
30 25 20
0.05
15 10 5 0
0 0
2
4
6
8
RET
TST
of influence between the compartments. This allows the application of network analysis techniques to determine the values of specific network properties. In particular, we found that network synergism appears to respond positively as fitness increases, and the other properties respond negatively. This paper represents the first attempt to integrate the two Boolean techniques; further research is needed to more deeply understand the interrelation between them.
10
time step TSS
699
References SDP
Fig. 5. Behavior of ecological goal functions for second simulation. TSS and RET are plotted on the alternate y-axis.
principles relating network structure and behavior. He found that system behavior ranged from very stable (‘‘deeply ordered’’ in Kauffman’s terminology) to chaotic depending on the number of nodes and connections, and that systems tended to evolve toward an intermediate level of connectivity (‘‘the edge of chaos’’). These results indicate that there is an optimal complexity for system design, which balances the opposing evolutionary tendencies toward increased stability (and hence less adaptability) and toward increased flexibility (with the concomitant risk of losing system integrity). In the ecological examples given above, we can see this tendency for systems to evolve toward an intermediate level of connectivity (Table 6). The network models described herein, although theoretical in nature, provide insight into system organization and due to their reliance on first principles of Boolean analysis are applicable to many different system types. Therefore, they provide a good platform for modeling integrated systems. This paper presents an initial foray into two different important network approaches with the aim of gaining better confidence in how they may be linked or they complement each other. The ecological network properties inform about structure, function, and relations, but they have traditionally been applied to static networks. Kauffman’s NK model gives a way to let the models evolve over time, and thus see how the properties change over time. Since the ecological goal functions are assumed to track ecosystem development, they have already been used for structural dynamic models (Nielsen, 1992) and model calibration (Jørgensen and Fath, 2004). It is also important to find other useful management applications for these methodologies. In the case of managing biological and societal systems interacting together and coevolving on correlated landscapes, a hierarchical management approach dealing with the subsystems shows potential to find a compromise solution between complexity levels of the different subsystems (Alexander, 1964).
6. Conclusions In conclusion, we have recreated Kauffman’s multi-species NK model and used it to investigate the coevolution of a simple model ecosystem. Furthermore, we have used the fitness values generated by the model as surrogates for the probability
Alexander, C., 1964. Notes on the Synthesis of Form. Harvard University Press, Cambridge. Cheslak, E.F., Lamarra, V.A., 1981. The residence time of energy as a measure of ecological organization. In: Mitsch, W.J., Bossermann, R.W., Klopatek, J.M. (Eds.), Energy and Ecological Modelling. Elsevier, Amsterdam, pp. 591e600. Fath, B.D., 2004. Network analysis applied to large-scale cyber-ecosystems. Ecol. Modell 171, 329e337. Fath, B.D., Patten, B.C., 1998. Network synergism: emergence of positive relations in ecological systems. Ecol. Model 107, 127e143. Fath, B.D., Patten, B.C., 1999a. Review of the foundations of network environ analysis. Ecosystems 2, 167e179. Fath, B.D., Patten, B.C., 1999b. Quantification of resource homogenization using network flow analysis. Ecol. Model 123, 193e205. Fath, B.D., Patten, B.C., Choi, J.S., 2001. Complementarity of ecological goal functions. J. Theor. Biol. 208 (4), 493e506. Fath, B.D., Jørgensen, S.E., Patten, B.C., Strasˇkraba, M., 2004. Ecosystem growth and development. Biosystems 77, 213e228. Finn, J.T., 1976. Measures of ecosystem structure and function derived from flow analysis. J. Theor. Biol. 56, 363e380. Hannon, B., 1973. The structure of ecosystems. J Theor. Biol. 41, 535e546. Higashi, M., Patten, B.C., 1989. Dominance of indirect causality in ecosystems. Amer. Nat 133, 288e302. Jones, C.G., Lawton, J.H., Shachak, M., 1997. Positive and negative effects of organisms as physical ecosystem engineers. Ecology 78, 1946e1957. Jørgensen, S.E., Fath, B.D., 2004. Modelling the selective adaptation of Darwin’s Finches. Ecol. Model. 176, 409e418. Jørgensen, S.E., Mejer, H.F., 1977. Ecological buffer capacity. Ecol. Model 3, 39e61. Jørgensen, S.E., Mejer, H.F., 1979. A holistic approach to ecological modelling. Ecol. Model 7, 169e189. Jørgensen, S.E., Mejer, S.E., 2002. Integration of Ecosystem Theories: A Pattern, third ed. Kluwer Academic Publishers , Dordrecht. Jørgensen, S.E., Patten, B.C., Strasˇkraba, M., 2000. Ecosystems emerging: 4. Growth. Ecol. Model 126, 249e284. Jørgensen, S.J., Ladegaard, N., Debeljak, M., Marques, J.C., 2005. Calculations of exergy for organisms. Ecol. Model 185, 165e175. Kauffman, S.A., 1993. The Origins of Order: Self-Organization and Selection in Evolution. Oxford University Press, New York. Kauffman, S.A., 1996. At Home in the Universe: The Search for Laws of SelfOrganization and Complexity. Oxford University Press, New York. Kauffman, S.A., 2000. Investigations. Oxford University Press, New York. Lotka, A.J., 1922. Contribution to the energetics of evolution. Proc. Nat. Acad. Science 8, 147e151. Marques, J.C., Pardal, B.C., Nielsen, S.N., Jørgensen, S.E., 1997. Analysis of the properties of exergy and biodiversity along an estuarine gradient of eutrophication. Ecol. Model 102, 155e167. Morowitz, H.J., 1968. Energy Flow in Biology: Biological Organization as a Problem in Thermal Physics. Academic Press, New York. Mu¨ller, F., Leupelt, M. (Eds.), 1998. Eco Targets, Goal Functions, and Orientors. Springer-Verlag, New York, p. 619. Nielsen, S.N., 1992. Strategies for structural-dynamical modelling. Ecol. Model 63, 91e102. Odling-Smee, F.J., Laland, K.N., Feldman, M.W., 2003. Niche Construction: The Neglected Process in Evolution. Princeton University Press, Princeton, 472 pp.
700
B.D. Fath, W.E. Grant / Environmental Modelling & Software 22 (2007) 693e700
Odum, E.P., 1969. The strategy of ecosystem development. Science 164, 262e 270. Odum, H.T., 1983. Systems Ecology: An Introduction. John Wiley, New York. Odum, H.T., Pinkerton, R.C., 1955. Time’s speed regulator: the optimum efficiency for maximum power output in physical and biological systems. Am. Sci. 43, 331e343. Patten, B.C., 1978. Systems approach to the concept of environment. Ohio J. Sci. 78, 206e222. Patten, B.C., 1981. Environs: the superniches of ecosystems. Amer. Zoo. 21, 845e852. Patten, B.C., 1982. Environs: relativistic elementary particles or ecology. Amer. Nat. 119, 179e219. Patten, B.C., 1985. Energy cycling in the ecosystem. Ecol. Model. 28, 1e71.
Patten, B.C., 1992. Energy, emergy and environs. Ecol. Modell 62, 29e69. Patten, B.C., Bosserman, R.W., Finn, J.T., Cale, W.G., 1976. Propagation of cause in ecosystems. In: Patten, B.C. (Ed.), Systems Analysis and Simulation in Ecology, vol. IV. Academic Press, New York, pp. 457e579. 593. Prigogine, I., 1980. From Being to Becoming: Time and Complexity in the Physical Sciences. WH Freeman and Co, San Francisco. Prigogine, I., Wiame, J.M., 1946. Biologie et Thermodynamique des Phenomenes Irreversible. Experientia II, 451e453. Ulanowicz, R.E., 1980. An hypothesis on the development of natural communities. J. Theor. Biol. 85, 223e245. Ulanowicz, R.E., 1986. Growth and Development: Ecosystem Phenomenology. Springer-Verlag, New York. Ulanowicz, R.E., 1997. Ecology, The Ascendent Perspective. Columbia University Press, New York.