Under what ecological conditions do species coexist, even though some of them replicate more slowly than others? An analogous problem exists in the evolution of metabolism. Below I review some standard results from previous ecological models. In particular I seek an understanding of how to demonste the co-existence of replicators.
Demonstrating Co-existence in Interspecific Competition Lotka-Volterra Models.
The classic Lotka-Volterra model of competition between N species is described here. The equations contain an exponential growth term multiplied by r1,the Malthusian growth-rate constant. If this were the only term, then the solution to the differential equation would be N(t) = N(0)e^kt. Instead, a reverse reaction is typically introduced, at rate proportional to N^2. This reverse reaction represents intra-specific competition. The intra-specific competition coefficient is set here to unity. In addition, the term on the far right represents inter-species competition.
For two species the equations are...
If alpha12 < 1, interspecific competition by species 2 on species 1 is less than intraspecific competition by species 1 on itself. The same applies to alpha21. Consider the phase space of this system for various values of K and alpha.
There are 4 possible fixed points (steady-states). The nullclines can have 4 qualitatively different configurations. In the first two, only one species survives, irrespective of the initial conditions. In the third, only one species survives, depending on initial conditions relative to the stable manifold of the saddle-point node in the center. In the fourth condition, both species coexist. This occurs when K1 > K2 x alpha12 AND K2 > K1 x alpha21, i.e. for equal carrying capacities, this means that the co-existence only occurs when the interspecies competition term is less than the intra-species competition term (1). Note that the means of demonstrating co-existence here has been to identify the stable states using the phase plots. This is an important technique.
Demonstrating Co-existence in Cross-Catalytic Networks of Autocatalysts with Unlimited Precursors.
The LV competition model is formally equivalent [DEMONSTRATE THIS] to 2nd-order replicator equations (Ref. 7) of the form...
in which replicators linearly influence the replication rates of other replicators, i.e. where fk(x) = the sum of alpha(kj) x Xj over all j. The subtracting sum on the right is a selection constraint in order to keep the total concentration of replicators equal in the reactor.
It is useful here to consider a simple example of the above dynamical equation, with two species, and show how it relates to the above Lotka-Volterra model in which we showed the conditions for co-existence above.
Catalytic networks (Ref 8). actually cross-catalytic networks) are subsets of 2nd order replicator networks in which alpha(kj) >= alpha(jj), i.e where cross-catalysis is always greater than auto-catalysis. Co-existence has been studied in such networks. For example, linear chains in catalytic networks do not show co-existence. The hypercycle is 'dynamically stable' but 'structurally unstable' to topological changes such as short cuts and parasites, the introduction of which leads to dynamic instability. These are two important and distinct varieties of stability.
The notion of co-existence has been addressed formally in different ways in replicator equation systems (Ref 9). From reference (4) "A replicator equation is said to be permanent if there is a compact subset, C, of the interior of the state space Sn, such that any trajectory starting in the interior of Sn will eventually end up in C. In other words, there is a (possibly very thin but nite) repulsive "skin" on the boundary of the state space." There are then conclusions that can be made about permenance from the topology of the catalytic network. For example, 2 disjoint hypercycles cannot coexist! (Ref 10).
Demonstrating Co-Existence due to Product Inhibition
Eors Szathmary (ref 12.) and Guenter Von Kiedrowski (Ref 13). have described how product inhibition leads to parabolic growth and can result in co-existence of replicators. Parabolic growth is described by dx/dt = k x^p, where 0<p<1. Often it is found that p = 1/2. Szathmary and Gladkih (1989) examined a set of replicators described by the equation below...
....(1)
Note that total population size is kept constant due to the sum on the right. Consider first the standard conditions of exponential growth where p = 1. METHOD 1 (VISUALIZE PHASE SPACE):Imagine two species x1 and x2, both with p = 1, but where k1 > k2, i.e. where x1 has higher malthusian growth rate than x2. At equilibrium we require that dx1/dt = dx2/dt = 0. This occurs at the fixed points {x1 = C, x2 = 0} or {x1 = 0, x2 = C}, i.e. there is never any co-existence. If k1 > k2, then x1 is stable. This is clear from the phase space of this system. The phase plot is drawn using CurveGraphics.nb for Mathematica.
On the Y axis is [X2], and on the X-axis is [X1]. We see that the trajectory moves linearly towards the positive [X2] axis, and [X1] = 0. This is consistant with k2 > k1. Note that the sum of X1 + X2 remains the same in the phase space, i.e. a trajectory can only move in a straight line in this system, since X1 + X2 always equals C. METHOD 2 (CHECK STABILITY OF ZERO FIXED POINT): It is also useful to be able to show that a species can invade when rare. In this case this was done by finding a function for x1'(t) in terms of constants and x1(t). Setting x1(t) small means we can ignore higher-order terms in x1(t), and so analyse the coefficient only in the 1st order term. This turns out to be (k1 -k2) which is always positive when k1 > k2, implying the x1 can increase from very small values. The opposite coefficient was found for x2. This approach is effectively linearizing at x1(t) =0 to show that 0 is an unstable attractor. Note that competitive exclusion in the above model depends on the mass conservation constraint. If both species could invade when rare, then dynamic co-existence would be possible.
Next consider the parabolic growth of just one species. First it would be helpful to understand under what conditions parabolic growth occurs, and why! Von Kiedrowski described the following "minimal" reaction steps; i.e. A + B + C <---> ABC ---> C2 <---> 2C. You can download the mathematica simulation of such a system here. The simulation code is reproduced below for the rates shown, a > b, a2 > b2, k << a,b. with an initially low concentration of C.
A and B decrease sigmoidally, and C increases sigmoidally. ABC shows a peak, associated with maximum rate of production of C. In fact it is clear that dC/dt = k[ABC] since ABC is the reactant immediately producing new C (in the double stranded C2 form). However, the new C2 cannot directly catalyse its own production but must dissociate to the single stranded C state. It is due to this inhibitory effect of C on its own production due to its tendency to bind to itself, and prevent incorperation of A + B, (which can bind only to single stranded C), that growth is sub-exponential. Compare the growth rate of C when the dissociation of C2 -> 2C is very fast, i.e. a2 = 1, and b2 = 0. See below...
In the above graphs, a2 = 1, and b2 = 0, i.e. dissociation of C2 into 2C is virtually instantaneous, and there is no product inhibition, i.e. no binding of C together to form C2. You can see that the concentrations of C2 follow the concentration of ABC. The rise of [C] is also much sharper in the second case. Below you see the concentration [ABC] + 2[2C] + [C] over time for various values of a2, from 1(virtually instantaneous in red), to 0.000001, i.e. very slow dissociation (black).
[ABC] + 2[2C] + [C] 
Time
The rise of C templates is slowed down if there is product inhibition. The point of inflection at which the second derivative of growth rate = 0, occurs at lower concentrations as product inhibition is increased, thus the above curves do not represent a "mere shifting to the right". It is useful to plot the initial growth rate dc/dt against the initial concentration [C] of template. Von Kiedrowski plotted the initial dc/dt against initial c, and found the square root relation c' = a + bc^0.5, rather than c' = ac, to be expected from exponential growth.
Below we derive the equations that demonstrate parabolic and sub-exponential growth at the begining of the run. These equations are important since by finding isomorphisms with other systems and these equations we will show that the other systems also exhibit sub-exponential growth.
[DO THIS!]
However, we do not need the above equations to show what happens GIVEN parabolic or sub-exponential growth. For a start we can integrate the growth laws, and obtain the following two solutions to the Malthusian and Sub-exponential growth laws.
The above solutions show that malthusian growht is exponential, and parabolic growth is, parabolic, i.e. proportional to kt^2.
So, given the existence of parabolic and sub-exponential growth, what are the consequences for selection and co-existence? The proof for co-existence devised by Szathmary and Gladkih depends on METHOD 2 above. They showed that "any varient of a subexponentially growing entity could invade when rare". Here I apply method 1 and method 2 to the parabolic case. Consider the phase space of two replicators x1 and x2, that grow parabolically. This gives the following picture. Both nullclines are overlapping, such that whatever concentration combinations one begins with, the system always moves towards co-existence. You can see that even when the species concentrations are very near 0, they move inexorably towards the nullcline.
Now I apply method 2. Note that due to the selection constraint, effectively we just have a one dimensional dynamical system, since x1 can be expressed purely in terms of x1! Note that all the trajectories lie along a line! See the following equations.
The top two equations are the application of (1) for two parabolic replicators, in our examples we have chosen k1 = 0.8 and k2 = 0.5. The third equations shows that the total replicator number is constant, and so replacing x2 with C-x1, in the equation for x1', we can obtain an expression for x1' purely in terms of x1. This is a one dimensional dynamical system which has the phase line shown below.
The same phase line can be obtained for X2, and looks the same except that the stable state is in a different place. There is a stable attractor at a positive value of X1 between 0 and 1, and unstable attractors at x1 = 0, and x1 = 1. This means that X1 can invade when stable. The same is true for X2, (phase line not drawn here), therefore both species can invade when rare, and hence dynamic coexistence is guarenteed, in stark contrast to the exponential case whose phase space was shown previously. Szathmary and Gladkih demonstrated the same result for any exponent p between 0 and 1. They call this property "Survival of Everybody".
Is there a relationship between the reason for co-existence in the LV equations where intra-species competition is greater than inter-species competition, and the above example of co-existence in parabolic replicators? Yes, in the second case, x1 inhibits its own growth always more than x2 can inhibit its growth. [DEMONSTRATE THIS].
Demonstrating Co-existence due to Spatial Structure and Diffusion Limitation
(Ref 14.)
Conditions for Co-existence in Recycling Autocatalytic Ecologies.
1. Competing Autocatalytic Cycles.
Consider the following chemical system in which two autocataytic particles Y1 and Y2 compete for the same food source X1, and produce the SAME waste material W1 which is recycled to X1 again. Whereas co-existence is always possible for independent non-interacting autocatalytic cycles that exist in their own seperate biospheric re-cycling loops, these two autocatalysts may not always be able to co-exist. The mathematica files for these systems can be downloaded here. The following reactions describe the system (Cellerator for Mathematica format).
The properties of an autocatalyst Y1 are described completely in terms of k1, k1' and d1. Therefore, we want to know how the autocatalyst 'phenotype' {k1,k1', d1} determines its competiveness with another autocatalyst 'phenotype', in this type of interaction. Some example runs are shown in the attached mathemetica file downloadable here.
Let us find first the nullcines for Y1 and Y2. These are given by setting the differential equation dY1/dt = 0, and dY2/dt = 0, respectively. We obtain [Y1]* = k1 [X1*] - d1 & 0, and [Y2]* = k2 [X1*] - d2 & 0, respectively. i.e. we see that line of values for which Y1 = 0, is linearly dependent upon the resting state of X1. A steady state with [Y1] > 0 and [Y2] > 0, (stable or unstable) can only exist if the resting value of X1 is greater than the smaller of the two ratios d1/k1 and d2/k2. This concentration (d/e) we will call the food threshold of that autocatalytic particle. So, next we need to know, for what values of k1,d1, k2 and d2, and the INITIAL CONDITIONS, [X1] is able to be maintained above min(d/k). We plot [X*], [Y1]* and [Y2]* against various parameters.
Plot 1 : {k1,k2} with both autocatalysts having the same decay rates e1 = e2. We used the following parameters, {A = 0.001; p = 1; R = 0.0001 ; k1d = 0.0001; d = 0.001, d2 = 0.001};{ k2d = 0.0001};, and allowed k1 and k2 to vary from 0 to 0.0015. The graph below LEFT shows the product of template concentrations, against various values of k1 and k2. We found that co-existence only occured when growth rates were equal. Even a small difference in growth rate is enough to send the slower growing species to very low steady state concentration. On the RIGHT is a graph of Log[concentration] for two species, k1 = 0.01, and k2 = 0.0105, a very small difference in growth rates. Neverthaless, k2 cannot persist at high concentration.
Plot 2 : {e1,e2} with both autocatalysts having the same growth rates k1 = k2. Let us set k1 = k2 = 0.01. All parameters are set asfollows : {A = 0.001; p = 1; R = 0.0001 ; k1 = 0.01; k1d = 0.0001; }; {k2 = 0.01; k2d = 0.0001}.
Above we see three ways of visualizing co-existence. On the left, the height = 1, if Y1 and Y2 are both present in concentration above 10^-10. We see that co-existence only occurs when the decay rates of both autocatalysts is equal! In the middle, height indicates the proportion on existing autocatalyst that is Y1. When the decay rate of Y1 is less than the decay rate of Y2, then Y1 is present at 100%. Note that when both Y1 and Y2 have high decay rates i.e. > 0.01, there is noise in this proportion since both species are present at very low concentration. On the right is shown the product of Y1 and Y2. Again, you see that high template concentration is only maintained for low decay rates, and co-existence is restricted to autocatalysts with equal decay rates.
Above is shown the steady state values of Y1, Y2, X1, and W1 respectively, for various values of autocatalyst decay rates, where both autocatalysts had the same growth rates. So in conclusion, it seems that even a small decay rate difference is enough to send the faster decay autocatalyst to a very low concentration, FOR THESE PARTICULAR INITIAL VALUES OF X!!!
Above you see that even for a very small difference in decay rates between Y1 and Y2, the faster decaying species is sent to a very low concentration. Thus we conclude that there can be no co-existence of (exponentially growing) autocatalysts with different decay rates, and equal growth rates, in the above re-cycling system.
But, we have not explored the parameter space fully, for example, we need to consider the effect of increasing initial X, [Xo]. Remarkably, as pointed out to me by Dov Stekel, by increasing [Xo] we alter the stability conditions for the above system. Even for "exponential replicators", if there is a reversible reaction Y + Y -> Y + X, then we can obtain co-existence of Y1 and Y2, even for different rates e1/k2 and e2/k2 as long as initial X1 is high enough. See below for a numerical demonstration and analytical proof of this.
This co-existence is not possible without a reversible reaction Y +Y -> Y + X. Of-course, such a reversible reaction cannot occur in templates since the activated monomers cannot be reproduced by hydrolysis of the newly formed strand back into two short single strands. The equivalent of this reaction for templates that contain an irreversible step in the autocatalytic cycle is decay of the inactive constituent, back to X! We show the analysis of these dynamics below.
A Proof of the Co-existence of Reversible Autocatalytic Cycles in Re-Cycling Reactors.
The equations describing the competing autocatalytic cycles in a recycling reactor are.
The following calculations are in a mathematica file downloadable here.
STEP 1: Simplify the System of Equations in Some Meaningful Way.
I have simplified the equations by removing W and allowing Y1 and Y2 to decay to X1 directly. This makes the system 3D rather than 4D. If there is no autocatalyst, then [X] does not change, as opposed to the full model in which [X] and [W] would take some time to come to a simple equalibrium determined by R/pA. I am now assuming that this equilibrium is instantaneous and results in [X]/[W] = pA/R, i.e. [W] = R[X]/pA. The results may be different from the above simplification if the initial conditions involve a long transient, for example, if the reactor is initialized with only W, and pA is very small. In this case the autocatalysts could die out before suffiient [X] could be created. Above, I have made the reasonable assumption that the biosphere has reached a re-cycling equilibrium before autocatalysts are introduced.
STEP 2: Determine the Fixed Points.
Next I determine the fixed points of the system by setting dX/dt = dY1/dt = dY2/dt = 0. I get 4 fixed points. At FP1, Y1 and Y2 = 0, At FP2 and FP3 Y1 and Y2 exist in isolation, and at FP3 there is "co-existence", as yet with indeterminate stability. Eventually, we seek parameter settings for which FP4 is stable.
STEP 3: Determine the Jacobian Matrix.
Next I determine the (3x3) Jacobian matrix from which I can obtain the Eigenvalues of FP4. This is shown in line 1 in the figure below. For FP4 I substitute the Y1 and Y2 values in terms of X[t] into the Jacobian obtaining line 2 in the figure below. When mathematica solves for the Eigenvalues we obtain some impenetrable nonsense. (3) for the Eigenvalues, of which there are 3, 1 being zero. This means that the determinent of the Eigenvalues is zero, which means that the fixed points are not isolated, but a line.
Strogatz has a very useful diagram in his book that classifies the stability of fixed points according to the trace and determinent of the linearized matrix at the fixed point, i.e. the Jacobian matrix calculated for values at the fixed point
STEP 4: Determine the Parameters and Initial Conditions for which the Eigenvalues are Less than Zero at Each Fixed Point.
Next I want to know for what values of paramaters and [X] these two non-zero Eigenvalues are less than zero, i.e. the conditions for stability of FP4, i.e. when the trace (the sum of eigenvalues) is less than zero. Solving the inequality (Eigenvalue1 + Eigenvalue2 + Eigenvalue3 < 0), for each of the Eigenvalues seperately we can obtain the conditions below.
[X] > d1/k1, and [X] > d2/k2.
for both non-zero Eigenvalues. Therefore, there IS a condition in which co-existence is stable, and that is (as we expected) when X > d1/d2 and X > d2/k2, i.e. greater than the autocatalytic threshold of both particles. However, if I now solve for the sum of the eigenvalues being less than zero, I get the expression, X must be greater than...
STEP 5: Try to Eliminate Variables from the System by Using Additional Reasonable Constraints.
But now we need to know an expression for [X] in terms of parameters and initial conditions. Since the system is conservative we have a means of doing this. We assume initially very low concentrations of Y1 and Y2. If the initial [Y1] and [Y2] << [Xinitial], then [Xinitial] = X[t] + Y1[t] + Y2[t]. Since we know the final values of Y1[t] and Y2[t] in terms of X[t], we can write the following equation by substituting expressions for Y1[t] and Y2[t] in terms of X[t] into the conservation equation above. Solving this equation for X[t] in terms of the parameters and Xi, gives us the equation below for X[t]!. This is precisely what we need in order to test the conditions for the inequalities that determine stability of co-existence.
Solving the inequality, we obtain an initial condition for Xi, above which co-existence is possible, assuming Y1 is the autocatalyst with the lowest threshold e/k. Xi must be higher as d1 decreases, and less as also k1 increases. In the case where d1 = d2, and k1 = k2, and k1d = k2d, i.e. when the replicators are identical, then [X] > d/k, as expected.
Solving the above inequality using the other value of X we obtained (i.e. the value of X at which the sum of eigenvalues is less than zero), we obtain the expression, Xi must be greater than,...
for stability.Presumably, we now no longer have to use a max term to calculate the largest ratio of e/k first, but can just plug both values. Both equations make the correct predictions for Xi for which coexistence is possible.
A more elegent proof would have been to prove the instability of the other three fixed points, from which the stability of the 4th fixed point follows. This would have involved a simpler set of Eigenvalue calculations, compared to the mess above. In conclusion co-existence of exponential replicators, with unequal growth thresholds, is possible ABOVE A MINIMAL INITIAL FOOD CONCENTRATION, in a recycling system if the replicators decay to FOOD supra-exponentially, i.e. in proportion to their concentration squared, i.e. if there is self limitation, as in the LV equations. Without this second-order decay term, co-existence is not possible. This is a potentially structurally stable co-existence since small changes to food threshold do not destroy co-existence.
Imagine that species decay to different waste products. The only effect of this is to alter the efficiency of recycling. Effectively, e can be considered the net rate with which Y is re-cycled back into food. Therefore, we need no further analysis of that case.
1. Parabolicly Growing Competing Autocatalytic Cycles, with Single and Double Strand Decay.
Above we have considered only an autocatalyst that grows according to the biomolecular rearrangement, X + Y <----> 2Y. In our later artificial chemistry this might be represented by aba + aab <------> 2aab for example. However, we can also envisage a more complex autocatalytic cycle that undergoes parabolic growth, rather than exponential growth, e.g.
aba + aab <-----> aababa
aababa ------> aabaab
aabaab <-----> 2 aab [Biased towards aabaab, i.e. product inhibition between aab molecules.
The diagram below shows an equivalent system. All reactions of Y1 and Y2 should be considered reversible, except for step XY -> YY. Effectively, this is the case that Lifson and Lifson (1999) (ref. 15). examined of an autocatalyst experiencing 'single-strand decay', i.e. Y1 decays, but XY1 and YY1 do not decay. Lifson and Lifson showed that under these conditions, when there were two parabolic replicators with single strands experiencing exponential decay, that the replicators survived together only when e1/k1 = e2/k2. Do we confirm Lifson and Lifson's findings in this recycling system? Yes. For the reaction network shown below, there is no co-existence except when the above condition is fulfilled. The mathematica file of the system below is available here.
For the parameter settings {A = 1; p = 1; R = 0.1; k1 = 1; k1d = 0.2; k12 = 1; k13 = 0.001; k13d = 1000.0; k2 = 1; k2d = 0.2; k22 = 1; k23 = 0.001; k23d = 1000.0; d = 0.1; d2 = 0.1}, I varied k1 and k2 systematically (bottom left), and found that co-existence occured only when they e1/k1 = e2/k2, for fixed e1 and e2, confirming Lifson and Lifson's findings.
Also, (Top right) I found that keeping one autocatalyst's parameters constant, and systematically varying the others e1 and k1 values, it was evident that for only a fixed ratio of e1/k1 was there co-existence, and this ratio was equal to e2/k2 = (1/0.1) = 10. Szathmary and Scheuring later demonstrated that when double strands were allowed to decay that co-existence could again be reestablished under certain parameter conditions. In particular, there must be a high X1 steady-state, or most of the strands must be in the double stranded state (see Szathmary and Scheuring (2001) for a proof of this). We implement decay of XY and YY strands, at rate dB = 0.002. This very effectively results in co-existence once again! See below...
As predicted, above we see co-existence of parabolic replicators when there is ds (0.002) and ss (0.1) decay. Thus, we expect that if two autocatalysts compete for the same resource within the generative chemistry,there should be dynamic stability only if there is decay of the active and inactive constituents of the autocatalytic cycle. Otherwise, if there is decay of the active autocatalyst only, we expect survival of the fittest, or "extinction of the unfittest", as Lifson and Lifson put it.
It is far superior to have an analytical demonstration of co-existence than a numerical one since one can prove a result for well defined parameter ranges. Szathmary and Scheuring (2001) managed to do this even for re-cycling systems of the sort that I have simulated above. It will be extremely useful to understand how they did this, so I will reproduce their reasoning below... I hope to be able to use similar analytical methods to deduce conditions for co-existence in the following cases also!
[SHOW ANALYTICAL METHODS HERE]
2. Scavenger Cycles.
A scavenger cycle is shown below. It has two replicators Y1 and Y2, Y2 feeds on the waste of Y1, and produces another waste material, W2. This is recycled back into W1. What are the conditions for co-existence in this case? The mathematica file for this analysis can be downloaded here.
The above network can be modelled using the full equations below.
STEP 1: Simplify the System in a Meaningful Way.
As before, a simplification is possible. Assume that re-cycling of Y2 to X occurs directly at rate e2. The system simplifies to a set of 4 differential equations below. For further simplificaiton remove W1, assuming that a proportion E of W1 goes to X1, and a proportion (1-E) goes to Y2. The only difference between this simplification and the predator prey case later is that the later growth of Y2 is not limited by the rate of decay of Y1 into W1, nor by the rate of re-cycling of W1 into X. Thus a completely new analysis will not be necessary, merely a re-casting of constants.
STEP 2: Calculate the Fixed Points.
The system has some rather complicated looking fixed points for Y1 and Y2.
There are 4 fixed points.
1. In FP1, both autocatalysts are extinct. This happens only when [X] = 0.
2. At FP2 the scavenger is extinct, but the primary species exists at positive concentration only if [X] > (d1xEf)/k1. It would be interesting to know what the basin of attraction of this fixed point was, and whether or not it was stable.
3. The next two fixed points FP3, and FP4 are the two roots of a quadratic characturistic equation. I can calculate the condition for the fixed point of Y2 being positive for FP3 and FP4, i.e. Y2 > 0.I find this occurs when X is greater than...
in both cases. This makes sense, since when Ef = 0, i.e. all of the flux from Y1 goes to Y2, and no flux goes directly to X, then this fixed point can exist only when...
[X] must be greater than the above term, for this co-existent fixed state to have positive Y1. Conversely, if Ef = 1, i.e. flux from Y1 is completely recycled back to X, without being able to feed Y2, then X must be infinity (division by zero), for positive Y2, this means, that there is no positive value of [X] for which this fixed state can have positive Y2.
4. However, there is no positive real solution for the value of X* for which Y1* = 0, for either of the expressions for the steady-states of Y1 in terms of X*. Consider first FP3, here we see that the Y1 is always positive for positive parameters and initial conditions. It only has a solution when d2 = 0. Conversely, the Y2 concentration at FP4 is always negative, tending to zero but never reaching it. This means that FP4 cannot exist for any positive concentrations of X, and so we can ignore it.
In conclusion, there are 3 possible fixed points in the positive orthant. FP1 has no species, Fp2 has only Y1, and FP3 has both. Next we determine the stability, and the basins of attraction of these fixed points.
STEP 3: Determine the Jacobian Matrix.
STEP 4: Determine the Parameters and Initial Conditions for which the Eigenvalues are Less than Zero at Each Fixed Point, i.e. the conditions for stability of each fixed point.
1. At FP1 there are three Eigenvalues. The first two are 0 and -d2. The third Eigenvalue is -d Ef + k1 X[t]. The fixed point is stable when the trace of eigenvalues is less than zero, (E1 + E2 + E3 < 0), i.e. when [X] < (d2 + d x Ef)/k1. Total extinction is not possible unless resting [X] is below this value.
2. At FP2 (Y1 is positive, but Y2 is extinct) there are again three Eigenvalues. They are show below.
The first Eigenvalue is zero.
The second eigenvalue is negative when [X] > (d x Ef)/ k1, i.e. when [X] is greater than the autocatalytic threshold of Y1. This is quite understandable, otherwise Y1 could not exist.
ii.The third Eigenvalue is negative when X is less than
. Presumably this must be a value of [X] above which Y2 could actually exist, but below which it cannot, hence making this a stable fixed point in the later case when [X] is below this value. So, calculating the trace E1 + E2 + E3 < 0, we find that this fixed point is stable when X is LESS THAN...
3. At FP3, the eigenvalues are extremely complicated and may be complex. IF [X] > the above value at which FP2 is stable,i.e. when FP2 is unstable, AND when FP1 is unstable, i.e. when [X] > (d2 + d x Ef)/k1, then FP3 must be stable, since no other fixed point is possible. The mathematica file can be downloaded here.
In conclusion, the stability of the fixed points depends on the resting value of [X]. All three fixed points can be stable given the correct resting values of [X]. We now attempt to an expression for the stability of FP3 in terms of parameters and initial conditions.
STEP 5: Try to Eliminate Variables from the System by Using Additional Reasonable Constraints.
Again, assume that initially Xo = X + Y1 + Y2, giving the same expression for X in terms of resitng X, Y1 and Y2, and the initial condition Xo. However, since the resting values of Y1 and Y2 are very complicated for FP3, we cannot easily get an expression for Xo under these conditions so instead we look at the stability of FP2.
In[5] shows the equality approximination relating the initial condition of X to the steady states of X at the attractor FP2. We see that at this attractor, resting X should relate to initial Xo as shown in Out[5]. Since we have an expression previously for the value of resting [X] above which FP2 is not stable, we can test for what value of Xo this occurs, In[8]. Out[6] shows that the initial X must be greater than the term on the right in order for there to be co-existence of Y1 and Y2. In[7] simplies this to produce Out[7]. However, we must also show that FP1 is unstable for this value of Xo. FP1 is unstable when [Xo] > (d2 + d x Ef )/k1. Therefore we conclude what when the two expressions for instability of FP1 and FP2 are satisfied, FP3 is stable.
Above you will see that co-existee i.e. stability of FP3 only occurs for small e2, and small e1. Note however that e1 can be much larger than e2, for co-existence to be possible, in fact two orders of magnitude larger in this case, where Xo is set to 1.
We have derived an expression based on two inequalities for the co-existence of scavenger-host populations in a re-cycling reactor.
We now test the above predictions empirically. We would like to plot numerically derived terms for co-existence as a function of e1 and e2.
Above we see that there is some agreement with the real system and the predictions. However, there is expected to be some deviation since the numerical integration is not carried out for infinity, and Y1 and Y2 initial conditions are not equal to zero, as assumed in the predictions. You can see the range of values for which Y2 persits takes a notch out of the concentration of Y1.
Below I simulate this system for another set of parameters.
The above graphs show the resting states of various concentrations as decay rates of each autocatalyst are varied. The parameter settings are as previously. You can see that co-existence occurs when e2 is low and when e1 is low.
The above graphs show that as initial [Xo] is increased, there is a greater extent of co-existence as predicted. The graphs show the proportion of strands that are Y1. As Xo is increased, the scavenger obtains a much higher proportion in the system.
Above we see the steady state values of Y1 and Y2 as initial [Xo] is increased. The graph on the top shows an inset for very small Xo values. There is a threshold above which Y1 grows, followed by another threshold above which Y2 also grows.
Observing the system as it settles down to the FP3, one notices that there are damped oscillations, suggesting the FP3 attractor must have been a complex one, with a negative real part.
Below I show the phase space of the system in 3D for different values of Ef from 0 to 1, initial conditions and parameters being as follows. { k1 = 1; k1d = 0.1; d = 1}. ,{ k2 = 1; k2d = 0.1; d2 = 0.01}; E[0] = 1, Y1[0] = 0.0001, Y2[0] = 0.0001
As efficiency is decreased so does resting Y2 decrease.
Below I show the phase space of the system in 3D for different values of d2 from 0 to 1, initial conditions and parameters being as follows. { Ef = 0.3. k1 = 1; k1d = 0.1; d = 1}. ,{ k2 = 1; k2d = 0.1; E[0] = 1, Y1[0] = 0.0001, Y2[0] = 0.0001 . As d2 is increased from 0 to 1 (i.e. equal to d1), so resting Y2 increases.
In conclusion, scavenger sycles can co-exist with their host if X can be maintained above a crucial minimum value calculated above.
It is useful now to consider how other people have analysed predator-prey dynamics, and to compare our conclusions with theirs. One major difference between this and the standard predator-prey models is that this system is conservative of mass and so cannot experience undamped oscillations. Also, since reversible reaction exist, there is self-limitation. The existence of predator-prey relationships with self-limitation is expected to be stable for a wide range of replication rates k1, k2 and decay rates e1, e2, in the presence of high food concentrations, especially. We can look for such relationships in our generative chemistry. Also, as suggested by Jon Rowe, I could initialize the generative chemistry seimulator with such ecological interactions and observe their stability. This approach is analogous to community construction models.
3. (Constituent) Predator-Prey Cycles.
A predator-prey cycle is shown below. Y1 is food for Y2.
The above cycle is analogous to the scavenger cycle except that now e1 x Ef is replaced by the efficiency of the predator. Often this is quite low, e.g. 10%, the rest of the material being recycled by some other means. This is formally therefore equivalent to the above model of host-scavenger interaction. No extra analysis is required. Co-existence of different replicators occurs under the same conditions as the above case.
4. (AND-product) Predator-Prey Cycles.
A varient of the above predator-prey cycle is where a side-product of the autocatalyst is used, instead of the active product. An alternative is to pray on the inactive product (not shown here). Either A or B could be food to Y2. In either case, we must re-assure the recycling of the side-products also. Is this system equivalent to Y2 utilizing Y1 directly, if the rate of recycling of B to A is equal to the rate of recycling of Y1 to X1? If it is, then no further analysis is required, and co-existence arises with much the same ease as for the above case.
5. Waste-Product Autocatalysis.
The waste product of Y1 itself could be autocatalytic, using another recycling food source.
This seems to be a rather different case to those considered above, because the particles Y1 and W1 use different food sources. The details of the analysis can be found in the following mathematica file here. The cellerator code is shown below...
along with the differential equations...
We assume a slightly simplified system in which there is no abiotic flux of X1 or X2 to waste. We carry our the same procedure as before. The system is 4 dimensional. Solving the differential equations is hard.
6. AND-Product Autocatalysis.
The AND side-product of Y1 cycle could be itself autocatalytic, using another food source X2.
7. AND-product interconversion Mutualism.
There could be mutualism between cycles, i.e. the side-product of one cycle is food to another cycle, which re-produces another necessary food molecule of the first cycle. There are some thermodynamic constraints that must be considered carefully here.
The two autocatalytic cycles of Y1 and Y2 are linked here by a co-factor, A/B. A special case of this is where the co-factor is autocatalytic for its own production, i.e. where Y1 = B. In this case, Y1/B can be thought of as an energized version of A. In this case, Y1/B can drive many other peripheral reactions and couple them to the autocatalytic core. We choose to concentrate on the dynamics of that special case, and the ways in which it can be coupled with other chemical reactions. These are shown in the diagram below.
The boxes show different ways of coupling a system to the autocatalytic core. Effectively, the top box shows a mutualistic linkage between autocatalytic cycles. The second cycle may or may not have external food source; here it is shown without it. Y2 would act as a sink that might destroy the original cycle unless Y2 was somehow recycled back into Y1, e.g. unless it was self-limiting, i.e. if Y2 + Y2 --> Y1 + A1 + Y2. Alternatively Y2 might be waste, and recycled externally back to O1 molecules. The second box shows a simple coupling where Y1 is changed back into A, whilst activating C into D. The box on the bottom right shows the same reaction, except where D itself undergoes some reaction with one of the species already in the reactor.
We should also be aware of how this work relates to the existing work on mutualism between species. Mutualisms involve benefits and costs to the interacting species (Holland et al 2002). I am dealing with metabolic mutualisms, so for example, the benefit to cycle Y2 is that it has a type of food from Y1 that it would not be otherwise able to obtain. In Box 1, (Top) there is actually no benefit to Y1, since no extra Y1 molecules are created due to the interaction with Y2. For this to be considered a mutualism it may be that Y2 would have to confer some benefit to Y1, otherwise Y2 would merely be parasitic upon Y1. What of the costs in the above systems? Let us introduce a positive effect of Y2 on one of the reaction rates e.g. Y1 + X ---> 2 Y1. In this case there is a benefit to Y1, and a cost to Y2. The benefit is that Y1 production is accelerated, and the cost to Y2 is that it is has a lower concentration of active particles. Of-course, there is also an indirect benefit to Y2, since by increasing Y1 production, it increases production of itself. Another possibility would be that it catalyses some other reaction that benefits Y1 but does not benefit itself. Here we are entering the world of second order replicators, extensively studied by Eigen, Schuster, Stadler and others. But we are also entering the world of ecological modeling which provides a set of new methods and insights.
Is there a simple abstract formulation of such interactions that allows us to determine which species are stable, without extensive parameter explosions? Let us consider such an abstraction briefly. Consider a network of nodes, each node has a real value A, and influences other nodes in the network either positively or negatively, and either weakly or strongly. [See the work of Rene Thomas].
One useful method for describing interactions between elements is the concept of "functional response". The functional response links the growth rate of Y1, to the concentration of Y2. An example is a diagram from Holland et al 2002 below. .
C indicates costs, GB indicates gross benefits. NE indicates net effects. They find different stability properties depending on the functional response. For example, if the next effect is saturating, this limits the unlimited growth of mutualists. Similar methods for defining functional response can be used in chemical settings.
The cellerator files for varients of the mutualistic system above can be downloaded here. Consider first the simplest version of mutualistic coupling between two autocatalytic cycles shown in the AND-product cycle. I consider two cases, one where the species compete for the same food source, and one where they have independent food sources. These are shown in the diagram below.
The constituents of each autocatalytic cycle is food for the other autocatalytic cycle, in Ganti's terms they are AND - coupled. The systems can be simplified to exclude A and B, we just consider the conversion of Y1 and Y2 into X1, or X1 and X2 directly. The concentrations Y1 and Y2 cannot increase in this system, since however much Y2 is produced, is taken up in the production of Y1, and vice versa. This has given us little insight. It is perhaps better to look at the generative steps that could have lead to the mutualism, and whether each step would have been of benefit to the autocatalytic constituents.
Lets try again, we start with two independent cycles each capable of sustaining itself, and consider the costs and benefits of undertaking a mutualistic interaction. Such a process is shown below. Assume there is some environmental abiotic chemistry at equilibrium, A X, B, and that there exists some autotroph capable of utilizing an external source of energy (chemical unlimited or light unlimited) to shift the equilibrium of the environmental equilibrium (1). Now assume that there is some heterotroph that depends on the environmentally equilibriated high energy molecule A to power its autocatalytic cycle (along with X1) to produce the low energy molecule B. Both species are limited by fixed environmental capacities, the first by the environmentals capacity to turn A into B, and the second by the environments capacity to turn B back into A.
The equations describing these two species (1) and (2), and the fixed environmental capacities, are shown below.
This page continues here.
CONCLUSIONS.
Replicators can be demonstrated to co-exist by...
1. Visualizing the phase space and finding stable fixed points with positive values of each replicator. This demonstrates dynamic stability of coexistence, but says nothing about how the system will respond to the introduction of novel reactions, i.e. structural stability of coexistence.
2. Showing that both species can invade from rarity by linearizing at [Species] = 0, i.e. showing that the attractor in which the species concentration is zero is unstable. Again, this is a analytic version of the above approach, and shows dynamic stability of co-existence.
The conditions for co-existence in the re-cycling system are...
3. When non-reversible exponential autocatalytic replicators compete for the same food source in a recycling system then only the one with the lowest decay rate to growth rate ratio survives. This is not the case if exponential autocatalytic replicators are reversible, i.e. if they can decay back into their food particle in proportion to their concentration squared. In this case there is a minimum threshold of food above which unequal replicators can co-exist. This is analogous to the LV competition model.
4. For parabolic replicators with only active constituent decay, there is no co-existence, since the inactive constituent can replenish the active constituent as the active constituent decays, and so serve as a buffer store of active replicators, specific for the species concerned. If on the other hand there is also inactive constituent decay, there can be survival of everybody for a range of inactive constituent decay rates. Again, self-limitation is established.
5. Co-existence of predator-prey, or host-scavenger systems is obtained, for a wide range of host predator kinetics, especially if the concentration of food can be maintained at high values.
6. Mutualism
References
1. Caswell, H. 1978. Predator-mediated coexistence: a non-equilibrium model. American Naturalist 112:127-154.
2. Richard R. Vance. (1985). The Stable Coexistence of Two Competitors for One Resource. The American Naturalist . Vol. 126, No. 1 (Jul., 1985), pp. 72-86 [Co-existing competitors can exceed the number of resources they compete for if, i). they use different life stages of prey, or different body parts! ii). use of different growing seasons, and MANY other means of spatial and temporal heterogenaity. Interference competition
3. Bärbel M. R. Stadler and Peter F. Stadler. Molecular Replicator Dynamics. (2002).
5. Self-replicating systems [(1) Novick J.S., Feng Q., and Tjivikua T. Kinetic studies and modeling of a self-replicating system. J. Amer. Chem. Soc., 113, 8831{8838 (1991). (2) Tjivikua T., Ballester P., and Rebek Jr. J. A self-replicating system. J. Am. Chem. Soc., 112, 1249{1250 (1990). (3) Famulok M., Nowick J., and Rebek Jr. J. Self-Replicating Systems. Act. Chim. Scand., 46, 315{324 (1992).]
6. Molecular Ecologies[(1) Wlotzka B. and McCaskill J.S. A molecular predator and its prey: Coupled isothermal amplification of nucleic acids. Chemistry & Biology, 4, 25{33 (1997). (2). McCaskill J.S. Spatially resolved in vitro molecular ecology. Biophys Chem, 66, 145{158 (1997).]
7. Hofbauer J. On the occurrence of limit cycles in Volterra-Lotka equations. Nonlin. Anal., 5, 1003{1007 (1981).
8. Schuster P., Sigmund K., and Wol R. Mass action kinetics of selfreplication in Flow reactors. J. Math. Anal. Appl., 78, 88{112 (1980).
9. Definitions of Co-existence[ (1) Schuster P., Sigmund K., and Wol R. Dynamical systems under constant organization III: Cooperative and competitive behaviour of hypercycles. J. Di . Eqns., 32, 357{368 (1979). (2). Freedman H. and Waltman P. Mathematical analysis of some three-species foodchain models. Math. Biosc., 33, 257{276 (1977). (3). Freedman H. and Moson P. Persistence de nitions and their connections. Proc. Amer. Math. Soc., 109, 1025{1033 (1990). (4). Hutson V. and Schmitt K. Permanence and the dynamics of biological systems. Math. Biosc., 111, 1{71 (1992).]
10. Hofbauer J. Competitive exclusion of disjoint hypercycles. Z. Phys. Chem., 216, 35{39 (2002).
11. Hutson V. and Moran W. Persistence of species obeying difference equations. Bull. Math. Biol., 15, 203{212 (1982).
12. (1) Szathmary E. and Gladkih I. Sub-exponential growth and coexistence of nonenzymatically replicating templates. J. Theor. Biol., 138, 55{58 (1989). (2). Scheuring I. and Szathm ary E. Survival of replicators with parabolic growth tendency and exponential decay. J. Theor. Biol., 212, 99{105 (2001). 2 shows that autocatalysts with decay are equivalent to subexponentially growing replicators. (3) SZATHMARY , E. (1991). Simple growth laws and selection consequences. Trends Ecol. Evol. 6, 366}370. (4). VARGA, Z. & SZATHMAD RY, E. (1997). An extremum principle for parabolic competition. Bull. Math. Biol. 59, 1145}1154.]
13. von Kiedrowski G. Minimal replicator theory I: Parabolic versus exponential growth. In Bioorganic Chemistry Frontiers, Volume 3, pp. 115{146 (Springer-Verlag, Berlin, Heidelberg, 1993).
14.Tereshko V. Selection and coexistence by reaction-di usion dynamics in fitness landscapes. Phys. Let. A, 260, 522{527 (1999).
15. S. Lifson and H. Lifson. A Model of Prebiotic Replication: Survival of the Fittest versus Extinction of the Unfittest. J. Theor. Biol. 199,425-433 (1999).