Download Latest Version Here.
I suppose this is a sort of science blog. It shows the thought process behind some simulations. This page should be considered a notebook of work in progress, where models of spatial replicators are hewn. People are welcome to contribute their own work to be linked, or can send corrections and comments. The page can be understood as a novel method for writing scientific papers. Any contributors will be authors. I intent to develop a wiki based system where anyone in the world can collaborate to write papers together, efficiently, openly, and honistly. One page on the wiki will be for discussion, results posting etc... another page will be for editing the actual paper. When all authors are satisfied, the paper may be submitted to a journal. Maybe this approach will massively increase productivity, and quality of publications. This is becuase a large number of people with complementary skills will be free to contribute, knowing that their contributon will be properly recognized.
Mistakes are always present, but hopefully decrease in frequency down the page. All the intermediate models and versions of the software are available online using the 'here' links. The above version is the latest version. However, versions downloadable from within the text are specific or each example.
Download a Spatial Network Simulator Backbone (compiles with XCode).
You can download this very simple skeleton program that simulates using Eular integration, a network of agents on a grid, with diffusion between grid squares. The initial versions produce unstable ecosystems due to self-self interaction coefficients being initialized to negative values. Later versions downloadable in the text, incrementally solve these problems.The code is in C++ and works with Cocoa on Mac OS X. It is intended to form a sort of template that you can use to write more complex (and more realistic) simulations. Below I describe some investigations. The version you downloaded above contains randomly assigned non-density dependent growth terms, and so generally results in unstable population explosions, see below. Xi is the concentration of species i, and kji is the effect of species j on species i. A matrix is kij values is randomly assigned from -1 to 1, and avector of xi values is assigned from 0 to 100. The number of species (m) = 5.
(Should be signma j = 0 -> m, not sigma i = 0 to m). The above growth law results in explosions of the population, see below...
Figure 1. Population sizes for a 5 species 'cell' (one square in the grid) with no diffusion. [y-axis = population size, x-axis = time]. The system is initialized randomly 6 times (see discontinuities) On the 5th initialization, a stable network results, however, on all other occasions there is instability, with species exploding. Species go extinct if their concentration drops below 0.001.
Effectively, the above model represents a simulation of interacting Malthusian replicators, where individual replicators can grow indefinately by posessing a positive k value for interaction with themselves, and/or form cycles of coupled species that are self-reinforcing for each others growth. Very unrealistic. A simple but powerful extension of the above model is to introduce density-dependent growth.
Single-Species Lotka-Voltera Equations
The LV equations are similar to those produced above except for the introduction of an extra term that reduces the growth rate of a species xi by a multiple (1-xi/K), where K is a positive constant. Thus, the growth of each species is inhibited by itself. When xi = K, there can be no more growth of xi, since (1-xi/K) becomes zero. The equation for a single species is quite simple.
Above we see that kii is the constant determining how much xi influences the rate of change of itself. K is the concentration of xi at which xi cannot increase any more. It is the carrying capacity of species xi. You can explore the dynamics of these equations with mathematica (download the file here).
The above simulation with mathematica shows that as K (the carrying capacity) is increased, so the resting state of y increases.
The above simulation with mathematica shows that as R is increased, so the rate at which the same carrying capacity (K) is reached is increased.
Multi-Species Lotka-Voltera Equations
LV equations can be used to describe interactions between species. Keeping the same structure as the above LV equation we subsume all effects due to other species on xi, in one term and multiply it by (1-xi/K) which represents the density dependence of xi imposed by itself, see below.
Above we see an equation for describing the density dependent growth of species, influenced by other species. kji is randomly initialized independent of kij. The growth depends on two terms. The first in the bracket is the multiplicative influence of other species (and for kii of xi itself) on xi. The second term in the bracket is the inhibitory influence of xi on its own growth. Only where the two terms cancel out is the growth of xi zero. You can download the second version of the spatial network simulator that generates and models randomly initialized networks of replicators described by the above equations here. Behaviour is still unstable because the sigma term can always be greater than the 1/K term given real parameters of kij and K. To ensure stability, i.e. that at some value of xi, there is no further increase of xi, it would be necessary to ensure that the sigma term could become less than the 1/K term with increasing xi.
The third version of the program downloadable here, contains an implementation of the general kind of interaction relation shown below. This resembles more closely an extension of the preditor-pray equations for two populations.
(The sum on the right should not apply for the case where j == i). Each species xi is described by the logistic equation (the first term), as well as being influnced by interactions with other species. The interactions may be mutually enhancing, preditor-pray, or mutually destructive, if kji can take values from -1 to 1. The equation is different from the previous one because this time there is a linear growth term for xi, which was not included before. Mutualism for example can arise if kji and kij are both positive for the same values of i and j. In this case, species i enhances the growth rate of j, and vice versa. Networks produced above are unstable. The dynamics tend to be unstable, since it is possible for the term on the right to override any density-dependent repression due to xi since it is not bultiplied by (1-xi/K).
More complex competition dynamics can be included if we allow xi values to interact within the 1/K term, e.g.
(Both sums should exclude j = i). This is similar to the previous equation except that there are other competition coefficients cji that affect how other species alter the carrying capacity of this species. There are many methods for constraining the ecosystem interactions to make them more realistic. Alex Penn discusses these in her thesis (Alex Penn, DPhil Thesis 2006). For example, D.S. Wilson has used a metacommunity model consisting also of patches in which within-patch dynamics are modeled as multispecies Lotka-Voltera models. (The prevention of diffusion results in increased differences between patches. You can try this with the above model also by switching off diffusion.)
Ni is equivalent to xi. S is equivalent to m. Ki is the carrying capacity of i as in previous equations. Ri is equivalent to kii. alphaij is equivalent to cji. This equation is similar to the above competition type model, without the preditor-pray component that we had considered. The implementation of the more general model, with mutualism, preditor-pray, and competition dyanmics included can be downloaded here. cij is a random number drawn from the range 0 - 1, whereas kij is a random number drawn from the range -1 to 1. The inclusion of non-zero kij's allows instability due to the fact that the carrying capacity term does not limit it. There is no diffusion between patches, each patch starts with a random population density vector. Thus, each patch can be thought of as a testtube holding the species, but with noise in the initial conditions. The dynamics are complex. Considering competition alone, there can be competitive exclusion or co-existence. There are potentially multiple stable equilibria. Possible extensions to the model are vast, e.g. make K dependent on the patch, K oscillates depending on environmental rhythms, kij and cij terms can be made variable dependent on species plasticity, the competitive effect can be made non-additive (e.g. A + B competes with C, but A and B alone do not succesfully compete with C).
Can Between-Ecosystem Selection Occur in the Wild?
Alex Penn has been doing important work on artificial ecosystem selection in her DPhil thesis. Here we look to see whether similar selection dynamics can occur 'in the wild', i.e. on a surface without explicit selection. The main point of this work is to generate time-series data that can be analysed statistically in order to automatically detect what might be considered 'higher-level' organizations present in ecosystem group structure. Although we use Lotka-Voltera equations, the code can easily be altered to model chemical reaction networks as well. However, LV systems are easier to deal with (guaranteed stability) so we use those for now. Spatial replicator models are not new, neither is the search for multilevel selection in spatial systems. For example, Hogeweg and Takeuchi have implemented a Hypercycle on a surface, which is however unstable to parasites when noise to replication rates is introduced. Stephan-Otto Attolini and Stadler have modeled hypercycle evolution in space. Spatial models by Scheuriing and Czaran et al of the co-evolution of replicators and metabolism show that co-existence of replicators is possible in space when there is diffusion limitation, due to the advantage of rarity. See also here for Szabo et al's Nature paper on this. Spatial models have been applied at the ecosystem level, e.g. Savill and Hogewag model preditor-prey waves on a surface. Spirals are found to be a higher order structural feature. Much more biologically realistic spatial ecosyste models have been implemented of-course, see Ault, Luo and Wang's model of spotted seatrout. A review of some work using the ecosystem modeling tool STELLA is found here. The spatial dynamics of logistic type growth equations that combine diffusion have been analysed for small networks by Lee et al. Effectively we are considering a specific class of reaction diffusion systems in which the diffusing entities are replicators described by some dynamical system. The distinguishing feature of this approach is a concern with the attractor structure of networks on each patch, and an analysis of the extent to which such networks, e.g. metabolic or species ecosystem networks can implement a holistic inheritance system without the need for genetic change. This type of evolution is likely to have been very important for the evolution of metabolism, prior to the existence of template replicators. It may also be important in ecosystems of species. Basically, anything you can do on a patch, you can do in space, and then you can see what happens. The results are not always intuitively obvious.
Assume that any one species consists of clones, i.e. all samples from the species will give us individuals with the same kii, cii and Ki values. Also, the kij, kji, cij and cji values will be the same, e.g. influence of other species on this species, and influence of this species on other species. Later we will allow this type of genetic variation and so allow microevolution of species interactions. Species are randomly distributed at first over a surface consisting of patches. The species diffuse over the surface. Interactions between species only occur within a patch. The values for cij, kij, K and kii are taken from Alex Penn's DPhil thesis. No skew of cij is used so far, rather the distribution is uniform random. Where kdiffusion = 0.01, in comparison to kii = 1.0, i.e. diffusion between patches is 100 times slower than replication within a patch, the dynamics observed in the program downloadable here are observed. Multiplicative noise has been introduced at each timestep of +- 50% the change of concentration.
The lightness of a patch relates to its total population density. The patch on the left above shows the resting state after all oscillations have died down. Kdif = 0.1, K = {100, 1000}, Num species (m) = 20, Initial Xi = rand()*100, kii = {1.5 - 2.5}, kij = 0, cij = {0,2.0}, T = 0.1. The patch on the right above shows another run without noise. There is still multiple values for resting population densities given the same set of interaction coefficients, carrying capacities and intrinsic growth rates. Note that this is not observed in all the runs. Some runs settle down into a completely plain grid, with homogenious total population density.
The surface reaches an equilibrium determined by the rate of species production on a patch (which obviously must differ from patch to patch), and diffusion between patches. It is interesting to note that the stable total population density of each patch is different. Is this because the species are in different attractors, or is there the same ranking of concentrations but with just a difference in the total concentration? To see whether this is the case, it would be interesting to visualize in some way the rank of population densities within each patch, to see how much the ranks of concentrations differed between patches. Alternatively one could do cross-correlation of rankings between each patch in the population! This would loose any idea of spatial structure however.
Note that in contrast to Alex's thesis, there is no explicit selection function due to pipetting, rather, species diffuse in proportion to their concentration. Also, there is no explicit fitness function that determines the extent of this proliferation of species, rather, species that are more effectively produced within the patch can obviously contribute to more diffusion than species that are not produced in high concentration within the patch. The question remains, is there just one attractor for all patches or are there different attractors? If there are different attractors then presumably it is possible for attractors to compete with each other for occupancy of the surface. An attractor takes over the surface occupied by another attractor by sending its species into the patch of the other attractor, and converting the attractor.
We want to know what the different patches are doing, e.g. what their compositions are, what the diversity of species is within that patch etc... First we visualize the entropy (H) of each patch using a different color coding system (pink). The modified program that does this is downloadable here. The entropy of species starts high, but often decreases, i.e. the diversity of the population tends to decrease. Sometimes at steady state there is variation in entropy between patches (left), with a few patches being able to maintain high diversity, see below. HIgh entropy is shown in bright pink, and low entropy is shown as black. Most often however, there is no variation between patches , of entropy at steady state (right).
A simple extension is to include 'higher-order' fitness contributions due to some feature of the ecosystem structure within a patch. Just as in Alex's thesis, we introduce a fitness benefit due to maximization of within patch entropy, i.e. species diversity. This is scaled such that maximum entropy = 1 and minimum entropy = 0. This is converted into a fitness effect by using the entropy measure to scale the carrying capacity of all species within the patch. The new carrying capacity becomes Knew = Kold(1+scaledEntropy). This means that patches with a maximum entropy have a carrying capacity double that of patches with minimum entropy. The code for this extension can be downloaded here. Pressing the start button over and over, one gets qualitatively different behaviour each time.
The first two pictures at the top show the stable states without diffusion, but with entropy dependent carrying capacities. You can see that the surface is now much more heterogeneous. Strangely, on the left, patches with high entropy have high total population density, but on the right, patches with high entropy, have low total population density.
The next two pictures show the system with diffusion. You you run the program, you see that with diffusion there are waves. It is not clear that the resting entropy is any different in the two cases.
Several questions need to be addressed, and this will require some quantitative analysis.
1. Does between-patch mixing (selection) dynamics improve the diversity measure over the surface? To test this, we must examine the diversity attained with and without diffusion. This is because diffusion is the means by which 'innoculation' between patches occurs. Without diffusion we are effectively doing random search from many different starting concentration vectors. We thus need to measure whether there is a statistically significant improvement in surface fitness as a result of diffusion. This would indicate that somehow innoculation is allowing conversion of some patches in low fitness attractors to a higher fitness attractor. [I need to plot a graph of diffusion rate verses resting total entropy in the system.]
2. Does including a skewed interaction coefficient distribution (i.e. more low strength interactions, few high strength interactions) alter the capacity for diffusion to improve total surface diversity? Alex found in her thesis that using skewed interaction coefficients did increase the capability of artificial selection to select for initial concentration vectors that maximized a fitness measure (a random linear sum of the concentration vector).
Another method of visualization is to calculate the 'central' (mean) ecosystem composition over all patches. Then to plot the distance from this mean position for each patch. An important question is, "How does one distinguish a rain forrest from grassland?". It is largely on the basis of species that are most hetergeneously distributed between the two. I suppose a rainforrest has trees but a grassland does not, although it is possible that other species concentration (or biomass) differences are greater between the two ecosystems. What is the equivalent in a vector of species concentrations? The easiest thing to do is to plot the frequency of the main species on each patch. Each species is assigned its own color (there are 20 species). The lightness of the color is the frequency of that species on the patch. Now this is truely amazing. Download the new code here. The static view dows not do justice to the beauty of the dynamics of this system!!!
See the 4 runs below. On the top right, at stable state there are two distinct ecosystem types, the first has blue species as max species, and the other has green species as max species. Entropy is low. On the top right, there are stable regions with high population density, high entropy and a purple species as max species. However (and this can only be appreciated with an animation), the darker regions are oscillating between several different top species. On the bottom left, the system has reached steady state, with a few cells having a different community composition. On the bottom right, all cells quickly attain the same community composition. Random initializations of the network result in these qualitatively very different behaviours. Carrying out statistics on randomly initialized networks would loose these distinctions. This is why visualization is useful. Please run the program if you can and see the beautiful patterns that result. (Click to download).
The most striking picture is shown below for a diffusion rate of 0.001 which seems to produce the most heterogeneous systems. Below you see that there are clearly distinct attractors in which the network exists, one with yellow species being high, the other with red species being high. The patches slowly change shape.
Whatever is happening, it is clear that re-organization does not continue indefinately. Eventually some stable surface configuration is reached. The question is, is this the optimum possible stable configuration, that maximizes the fitness function?
An important manipulation that was found to be important (Penn, 2006) was the production of skewed interaction matrices according to the function kij = K*x^D, where K is the maximum value of the interaction matrix, x is a random number between 0 and 1, and D is a constant from 1.5 -> 3. Some networks were simulated with these types of interaction matrices also. You can change the line of code that generates the random distribution quite easily in world.cpp.
Analysis of Fitness Maximization As a Function of Diffusion Rate
We try to answer the question whether diffusion rate results in a maximization of the total surface fitness. In order to do this, we generate 500 different interaction matrices, and undertake the following experiments on each one.
1. Initialized with a uniform random distribution of species, what is the final total entropy after a run where diffusion rate = {0,0.0001,0.001,0.01,0.1,1,10}? Plot the final entropy distribution for each value of diffusion rate. [y-axis = final total entropy per trial, x-axis = initial total concentration of species on the surface, or x-axis = total sum over the interaction matrix, or x-axis = total sum over carrying capacities].
2. Using tha above data, plot the mean entropy over all trials, for each diffusion rate value. [y-axis = mean total entropy, x-axis = diffusion rate.]
3. Add lines onto the above graph for different values of skew, {D = 1, 1.5, 2. 2.5, 3}. Include error bars.
4. How much variability is there between interaction matrices, and can their behaviour be clustered in any way?
Download the code to run the batch trials here. Batch results are written to batchData.data. Some data generated with the model is made available below. You can analyse this data yourself if you like. The mathematica code I used to analyse the data is available here.
1. Dataset 1. Trials conducted with this the program above, with D = 1.0. The format of the data is space seperated with the following fields.
NetworkNumber DiffusionRate TrialNumber TotalSurfaceEntropy
Results
Below you see the total entropies over the surface (H) after 3000 time steps of size 0.1, for 500 randomly generated networks. Each network was tested with 5 different diffusion rates from (0.000001 to 0.01). The average replication rate of a species was 1.0. The results below are for D = 1.0, i.e. interaction matrices with no skew. You can see that there are 2 qualitatively different types of result. The first kine of network shows no response to diffusion rate w.r.t. increasing resting total entropy. The second kind of network does show increased entropy with higher diffusion rates. These two classes of network may better be visualized if we plot the maximum difference between total entropies per network vs. the difference in diffusion rates that resulted in that maximum difference. From the graph below it even looks like there may be an optimum diffusion rate of 0.001, at which resting entropy of the surface is maximized, but this needs to be confirmed.
We will visualize all the other datasets with this kind of graph and then try other visualization techniques. This method of visualization is obviously not very useful due to the large number of boring networks that obscure the few interesting ones.
There is no obvious effect of D visable from these graphs. The most obvious thing to make this data easier to see is to take the means over all randomly generated networks, for each diffusion value. This graph is shown below.
Above you can see that the mean final H values plotted over all 500 randomly generated networks for each different D (skew) value, at different diffusion rates. Higher D values result in higher entropies, irrespective of diffusion rate. The increase in final entropy as a function of diffusion rate is greatest for the low D values, but there is an increase in mean H with increasing diffusion rate, for all D values.
The graph below shows (on the x-axis) the maximum difference in final entropy obtained in each network, when that network was run with different diffusion rates. On the y-axis is shown the difference in the diffusion values at which the maximum and minimum entropies were obtained. The graph is for D = 3.0. (i.e. high skew). Similar graphs are observed for the other skew values.
We see that for the networks that had the greatest difference in entropy due to changing the diffusion rate, the diffusion rate that produced the greatest entropy tended to be 4 orders of magnitude higher than the diffusion rate that produced the lowest entropy. We would like to have a graph of mean D(max H)- D(min H) value for bins of Max H - Min H. We would expect that as the Max H - Min H increased, there would be an increase in mean D(max H) - D(min H). The graphs are shown below for each value of D> The red line is a mean value of D(maxH) - D (Min H) over 10 unit x-axis ranges, the blue line is a linear regression fit, and the green line is a quadratic regression fit .
You can see that the largest differences in entropy occured between runs with the greatest difference in diffusion rates. Only a few randomly initialized networks had large differences in entropy due to diffusion. It would be most interesting to analyse these particular networks. In the program downloadable here, we add noise to diffusion as well as replication, and we save each network to a file in a format that can be 'cut and pasted' into the function world::world(1). This is a rather tacky way of doing this, but it allows some flexibility, and it produces a readable file. If you want to adjust parameters, this is done at the top of "world.h".
Is Ecosystem Selection Really Happening?
Unfortunately, the above experiments do not show that 'selection' has been important. It is necessary to show that the selection coefficient promoting high entropy on each patch actually alters the resting total surface entropy. This means re-running the above experiments with different values of SEL. If selection is affecting the optimization of surface entropy then we would expect higher SEL values should give higher resting total surface entropies. This extra control is required since the increased surface entropy could have been merely a result of increased diffusion rate.
The above graph shows the surface entropy without higher order fitness effects (red) and with higher order fitness effects doubled (blue). We see that with higher order fitness effects there is increased final surface entropy for all diffusion rates except for very high diffusion rates of 0.1. With selection coefficient = 2(blue) the maximum final surface entropy is greatest. The above results are for D = 1.5. The results arenot striking suggesting that this kind of diffusion alone is insufficient for 'naturalizing' selection dynamics.
A fitness function that might be less well favoured by diffusion alone is the randomly assigned linear-weight fitness function.
It would be useful to analyse the dynamics of the networks that displayed interesting behaviour. There are several ways to visualize the attractor structure of ecological networks. One method is to exhaustively explore the attractor structure of the network, i.e which attractor is reached for an evenly spaced sample of initial conditions?
Allowing Genetic Variation within Species
1. A variation of the ecosystem selection model is to assume that within a species there is genetic variation. This is implemented as slight noise in the interaction coefficients of species due to diffusion. Of-course, there should be genetic variation within patches as well, so we should observe noise in the interaction coefficients within the patches. Diffusion itself should sample from these within patch interaction coefficients. As an approximation we introduce noise in the within patch interaction coefficients, and allow diffusion to 'copy' this noise from one patch to another in proportion to the extent of diffusion.
Genetic Variation Between Species (Mutation) Verses the Introduction of Novel Species (Immigrants and Tourists and Peacekeepers).
2. Another variation of the ecosystem model is to assume that new species are introduced into the population, at some point or points on the surface of the ecosystem (invaders). These new species are assigned a random set of interaction coefficients and entered into the appropriate data structures. The simulation is run for some time, and then another species is introduced. A simple version of such a program is downloadable here. A species is introduced every 1000 time steps. No attempt is made to check to equilibrium before the introduction of a new species.
This is a very simplistic community interaction model. Other options are cascade or niche approaches. Some runs from the community iteration model downloadable here are shown below.
Above we see the number of species (diversity) of the ecosystem over time, for a surface of 20x20 patches with diffusion = 0.0001, D = 1.5, and a new species being introduced every 10 time units (in the scale of the above graph), equivalent to every 1000 time steps. a selection coefficient of 1.0 is used, i.e. species with double the entropy on a patch have double the carrying capacity. The individual species concentrations are shown below for the same run.
You can see that there appear to be periods of stability, in which invaders are not able to destabilize the existing ecosystem, however, inevitably, the invaders win after some time. The run ended when the number of species increased to 26, and the eular integration failed. We looked at the principle components of the entire system over time and found that there were quite clear clusters, and periods of steady movement through phase space. This is visualized below in a scatter plot in the dimensions of the first two and three principle components. The colour of the circle represents the number of species in the ecosystem at that time point. Blue circles are ecosystem with few speices, and red circles are ecosystems with many species.
Unfortunately, visualizing the ecosystem in this way looses information about all the structure of the surface, and really one must look by eye to the surface to see the interesting topologies that arise. For example, sometimes two ecosystem types co-exist on the surface, not being able to invade each others territories. Each ecosystem type is distinguished by containing a different distribution of species. The fact that new species are only injected into some of the patches on the surface contributes to this diversity.
Webworld: A method of constraining ecological interactions to make them much more realistic!
Caldarelli, Drossel, Higgs, McKane and Quince developed a model that combines evolutionary dynamics with ecological dynamics. Mikel Maron did his MSc thesis (University of Sussex) on an extension to this fascinating system and we explore this system below. Species are defined as having a set of characturistics, e.g. a binary string. Interactions between species (strength, direction) depends on the interaction between binary or integer strings, of K possible features and length L. Typically, L = 10, and K = 500. Characturistics are for example 'sharp claws' or 'red hair' or 'long neck' (or 'nice arse'?). Speciation, i.e. the generation of new species from old ones, is an operation on this characturistic list, therefore unlike above where a new species has a randomly generated interaction matrix, species here have interaction matrices constrained by the species that already exist. N.B. Such a generative story might be told for a chemical reaction network in which characturistics are active sites or groups on the chemical in question. The system is initialized with an antisymmetric KxK feature matrix of normally distributed random numbers, i.e. Mij = -Mji, where Mij tells us how feature i performs against feature j. The following equation describes the Score Sij of species i against species j...
Thus, the species with the highest score is the preditor, and the species with the lowest score has its score shifted to 0, it is the pray. This score is used to constrain the interaction matrix. External resources enter the system at rate R, and distributed to species according to Sij scores. The abiotic environment is species 0, and has a fixed set of characturistics. Speices with Sie > Sei are primary produces, where e = 0, the environment. Higher R results in larger and more diverse communities. The ecological dynamics are described by the following equation.
There is a constant death rate, uniform and equal. The second term is the number of resources given to species i from species j, and depends on gij, the 'functional response'. Lambda is the ecological efficiency (0-1), generally Lambda = 0.1. The final term is predation by species j on species i. gij is given by the equation...
extended to several species. b describes the level at which saturation of predation occurs. With multiple preditors this becomes...
alphaki represents competition between and within species. It is 1 for intra-species competition, and < 1 for inter-species competition. Its value is determined by the degree of similatiry between species, qki, which is the proportion of features the species have in common, see equation...
alpha_ki = c + (1-c).q_ki,
where c is the minimal level of competition. c has range [0-1]. Increasing the minimum competition decreases food web diversity. The response function gij is further modified by considering the fact that predictors cannot apply unlimited effort to predate multiple pray. Rather, it is an evolutionary stable stratergy for them to divide efforts according to which pray afford the greatest transfer of resources. If maximum effort (fij) must sum to 1, maximization of resources can be obtained by...
resulting in the final equation for gij(t).
It is necessary to iterate the calculation of fij and gij for all species until convergence, whereupon gij can be used to calculate the change in populations. This is all rather complicated and it will be helpful to have a simple example. The advantage of the webworld constraints is that mutations cannot cause independent changes to the interaction matrix Sij, i.e. one species cannot specifically increase its predation of other species and reduce is predation by other speices. A species can only alter its feature, and this results in a normally distirbuted set of mij values which is used to calculate Sij. This is meant to produce a realistic constraint on the capacity for adapations to produce 'magical' species, uber-preditors. [Try this in the previous community iteration model and see what happens.] Click here for the paper we have submitted on an extension of webworld.
The 'Evolution' of Metabolism
How could molecules as complex as nucleotides arise? Can we concieve of a self-organizing process occuring in chemistry that might somehow promote increasingly complex molecules? Maximization of Entropy Production Principle suggests there might be such a principle, however, this is a thermodynamic theory, and so it does not address the kinetic constraints that might be crucial in determining how fast any given system may be able to maximize entropy dissipation. Why might increasingly specific catalysts arise? Might chemical 'individuals' have similar organizational properties to 'ecosystem individuals',e.g. multiple stable states. Of-course, the ecosystem work above has defined two stages of evolution. The first experiments showed that selection and heredity could occur between ecosystem attractors of the same chemical network. Later we investigated the processes that could actually change the range of possible attractors that an ecosystem unit could have, these included genetic micro-mutation resulting in heterogeneity in the interaction coefficients of individuals within a species, however, we also considered the possibility of 'invasion' by external species. Although in chemistry there exist autocatalytic individuals, it has been assumed that there is little scope for 'micromutation' in chemistry, for example, a chemical X that is capable of autocatalytic growth is generally unlikely to be capable of autocatalytic growth given some chemical modification to its structure. There may be some chemical systems that exhibit this property to a small extent however, e.g. autocatalytic chemicals that produce other slightly modified autocatalytic chemicals. In order to investigate these chemical systems, I developed a simple model of chemistry that conserves mass and energy and implemented it on a surface with diffusion.
Click here to view the Evolution of Metabolism Paper in Progress.
Automated Detection of 'Autonomous' Entities.
Is it possible to write an algorithm that detects 'ecosystem individuals', i.e. groups of species that grow and die together in some sort of unit that has higher order fitness benefits? In the graph below I plot the total concentration of species over the whole surface during a run.
This corresponds to the run shown below in which one type of patch takes over another. The results of some other runs are shown below.
In the above run, one type of ecosystem patch takes over the surface, whilst another type dies out. D = 3.5. Diff = 0.001
The most complicated dynamics exhibited by the surface seems to be the eventual takeover by one ecosystem attractor of all the other types of patch on the surface. More interesting dynamics arise when within species getting change is permitted.