Control Experiments to Test the Functioning of a Reaction-Diffusion Model.

 

Download the Simulator:

Newest version here (19th July 2006). Autocatalytic Spot Production.

The initial C++ code for the model can be downloaded here (version 1) for XCode on Mac OS X. It is recommended that you download the update that was used to obtain the particular result shown below. This will be indicated. You will need to also download the chemical network file available here. New chemical network files will be provided throughout the page as new chemical networks are introduced.

and then simply set the line below in world.cpp to the directory in which you put this file,

//1. *********************FIRST INITIALIZE THE CHEMISTRY AS WAS DONE ABOVE.*****

srand(time(NULL));
ifstream myfile ("/Users/ctf20/Documents/MetabolismEvolver/RecyclingEvolver14thJuly/Release2/chemicalReloadable3.data");

Shown in red is that part you need to modify in world.cpp so that it can find the chemical network file.

Versions.

1. The Von Neumann neighbourhood was extended to the Moore neighbourhood. This removes the 1 patch period crosses that are stable otherwise. This modified code is downloadable here. (18th July 2006).

2. Simple autocatalytic reaction with decay ona surface. Download here (19th July 2006).

3. Autocatalytic reaction of A that repels the side-reactant R, producing patches of A with low steady-state light absorption. (19th July 2006). Download here.

Diffusion Method:

The surface consists of patches. Each patch runs an ODE of chemical reactions to calculate within patch update at each timestep, and contains a continuous vector of chemical concentrations. Diffusion is determined as follows. A physical potential is determined for each chemical species i on each patch. The potential of species i on a central patch is the sum over the product of the repulsionMatrix[i][j] and the concentration of species [j], for the central and the adjacent patches. Once the potential for species i has been calculated for all patches, the potential difference is used to calculate the propensity of diffusion between patches according to the equation [i] . dP/(e^dP - 1), where dP is the potential difference between the central patch and the patch to which we are calculating the extent of diffusion, and [i] is the concentration of species i on the central patch. This is determined with each patch simultaneously and deterministically.

Results:

1. Results obtained with version 1 (Phase seperation of A and B). The initial simulation results below test the diffusion part of the simulation for a system consisting of two non-reactive molecules a and b. Tims-step = 0.0001. Diffusion constant = 10 for both A and B. The repulsion matrix is shown below.

	A	B
A	0.0001	0.01
B	0.01	0.0001

A colorMap is used to visualize the concentrations on the surface. Red represents [A], and Green represents [B]. The surface is initalized with random concentrations of A and B, uniformly distributed from 0 to 100. The following series of pictures shows the process of phase seperation taking place on a 20x20 toroidal grid. This takes about 1 minute between pictures on a iBook G4.

There is a process of amalgamation of initially fine grained phase islands and channels into islands and channels of larger scale. Notice also the dark region between the red and blue regions which is an area with fewer molecules in total. This is because at the interface the potential for diffusion is high for both A and B compared to within the center of a phase seperated patch. The speed at which further amalgamation occurs decreases as the blob size increases since the characturistic diffusion scale dictates the flow. More complex structural forces arising from macrostructure in the system are not explicitly considered.

2. Results obtained with Version 2. Introduction of Autocatalytic Chemical Reactions. The following reactions are introduced. They are all bimolecular reactions, i.e. reactions capable of exergonic-endogonic coupling.

X + A <---> A + A

A + R <----> W + W

W + W ---> X + R (The extent of recycling of the system can be measured by looking at the steady state rate of this reaction, since it is the only irreversible reaction in the system.)

The first reaction is autocatalytic growth of A, by using food X. The second is decay of A by combining with the side-product R to produce two waste molecules W. W can be converted by light back into X and R. Download version 2 that runs the above system here. The new network file is required also, here.That file specifies the network and is written by the void Patch::saveReloadableChemistryAndPatch(void) function in chemistry.cpp.

The chemicalnetwork file contains the following information

Reaction, CatalyticFactor, (Name Number EnergyDefined FreeEnergy AbsorbsLight) for reactant 1, reactant 2, product 1 and product 2, followed by, ReactionAbsorbsLight, HeatReleased, Kf, Kb.

R 0.01 ba 0 1 5 0 ab 1 1 1 0 ab 1 1 1 0 ab 1 1 1 0 0 4 0.049715 0.01
R 1 ab 1 1 1 0 abbb 2 1 0.1 0 abb 3 1 0.01 1 abb 3 1 0.01 0 0 0.01 0.0154188 0.01
R 1 abb 3 1 0.01 1 abb 3 1 0.01 1 ba 0 1 5 0 abbb 2 1 0.01 0 1 5 1 0

The catalytic rate of the autocatalytic reaction is set to 0.01, so we expect a lower steady state concentration of autocatalyst.

The physical properties of the network are defined in the chemistry constructor in chemistry.cpp

//HAND=DESIGNED REPULSION MATRIX INSERTED BELOW.
repulsionMatrix[0][0] = 0.001; //A-A
repulsionMatrix[0][1] = 0.001; //A-B
repulsionMatrix[0][2] = 0.001; //A-C
repulsionMatrix[0][3] = 0.001; //A-D

repulsionMatrix[1][0] = 0.001; //B-A
repulsionMatrix[1][1] = 0.001; //B-B
repulsionMatrix[1][2] = 0.001; //B-C
repulsionMatrix[1][3] = 0.001; //B-D

repulsionMatrix[2][0] = 0.001; //C-A
repulsionMatrix[2][1] = 0.001; //C-B
repulsionMatrix[2][2] = 0.001; //C-C
repulsionMatrix[2][3] = 0.001; //C-D

repulsionMatrix[3][0] = 0.001; //D-A
repulsionMatrix[3][1] = 0.001; //D-B
repulsionMatrix[3][2] = 0.001; //D-C
repulsionMatrix[3][3] = 0.001; //D-D

The above system has no tendency for phase seperation since all combinations of species pairs are equally repulsive. [A] is the extent of red, [R] is the extent of green, and [W] is the extent of blue. W is present in lower quantities than A and R and so there is not much blue. The first control is to visualize the behaviour without diffusion, see below (part 1). Initial concentration of X, A, and R is set to a uniformly distributed random number from 0 to 100. Initial [W] = 0. [A] settles down to a low value, which does not go to zero due to the continued presence of A in some patches by the end of the experiment. This is presumably because some patches start with very low concentrations of R.

In part 2 diffusion with a diffusion constant of 0.01 is introduced. In this case [A] decays completely.

Part 3 shows the case where A and R repel each other 100 times more strongly than the repulsion of any other pair. By doing this I hope to reduce the effective concentration of A and R and hence establish islands of A that can maintain themselves in an excluded sea of R. This is indeed what we observe. Although A initially decays, the establishment of phase seperation allows it to grow in total concentration. Larger and larger patches form. There is a very interesting bump in the concentratiion of [W], which appears to be associated with the aggregation of two 'spots' into a larger spot, or due to the destruction of a spot. You can see from the surface picture that there is a border of low concentration around the spot of A, and a region of high R concentration surrounding that.

The code and chemical network file that produced the self-maintaining spots can be downloaded here and here!

I would now like to visualize how the rate of light absorption varies over the surface. This is given by the concentration of [W] since light is only absorbed by this species. The low concentration of W means that steady-state light absorption is low.

Above the color map is altered so that W is represented more strongly. We see that W concentration decreases steadily over time. [W] is high in the interfaces between red A spots and the R sea. It is often present in high concentration in the region between A spots and the R sea. This is to be expected since it is only the interaction between A and R that can produce W, and this interaction only occurs at the boundary. It is at the boundary that the reaction responsible for the conversion of external light energy to chemical energy takes place in this system. The system tends to a state with one large patch of A. This allows the smallest boundary which reduces the physical potential between A and R. We would expect the system to tend to this state in time. In that state, the extent of light absorption is determined by the negligable diffusion that occurs between the inside and the outside of the spot.

We confirm the cause of the blip in [W] by observing that the amalgamation of two patches is associated with a temporary increase in W production, see below...

But, previously we have been interested in maximizing the extent of light absorption, i.e. maximizing the steady state concentration of W. How can this be done in the surface system? One option is to use a GA to find solutions to this problem. Another is to think about it. What happens if the rate of diffusion is increased? This might be expected to increase the extent of A and R that reacts. What physical parameters need to be altered in order to maximize the rate of W production? A and R must react at the greatest rate possible. This cannot happen in a perfectly mixed system since R will destroy A, and it cannot happen in a perfectly phase seperated system since R cannot react with A except at the boundary.

Below I examine the state of the system after 20000 time steps with different diffusion rates. Is there an optimium diffusion rate at which light absorption is maximized?

Dif	[W] after 20000 ts.	[A]
0	9.48545	568.872
0.001	2.81521	1596.98
0.01	0.719684	14398.1
0.1	0.580884	19779.2
0.5	0.555288	21171

The high diffusion rates permit the establishment of a phase seperation that reduces the rate at which chemical equilibrium can be reached. With low diffusion rates the chemical reactions happen in seperate well mixed reactors, and there is no possibility of phase seperation contributing to reduced decay rate of A. The code used to produce the results below is downloadable here. The repulsion rates can be set in chemistry.cpp in the chemistry constructor.

Below I examine the state of the system after 20000 time steps with different repulsions between A and R at diffusion rate 0.5.

Repulsion	[W] after 20000 ts.	[A]	Structure	Total [W]
0.001	0.529737	0.0000102761	No spots
0.01	0.529238	0.0000098502	No spots
0.02	0.623618	0.0000141692	Spots form but are lost.
0.0210	~0.6	~0	Spots form and decay "brightly"
0.0215	1.34281	0.000060776	Spots form and decay	23903.7
0.0215	3.56277	0.000473748	Spots form and decay "brightly" slowly.
0.0217	86.2571	211.32	Single spot coheres and decays bright
0.0218	2.187	0.000182222	Spots form and decay.
0.0219	134.528	17704.1	Large spot formed, growth!	37633.0
0.0219	146.473	13015.3	Large spot formed, growth!	37268.4
0.0219	109.126	24867.8	Large spot formed, no decay.
0.0220	120.173	24614.1	Large spot formed, slow decay.
0.0225	113.405	25666.5	Spots form and cohere.
0.025	81.9911	27519.2	Large spots form, decay slowly?
0.03	41.0018	28789.9	Spots form and cohere
0.04	13.1913	27849.9	Spots form and cohere
0.05	4.72786	27029.8	Spots form and cohere
0.05	4.68326	28377	Spots form and cohere into mesh	1331.55
0.1	0.555288	21171	Thin stripes form and cohere very slowly.
1.0
10

It seems there is a threshold repulsion between A and R above which spots are stable so that [A] is high and [R] is also high. However, as the repulsion is increased further beyond this threshold both the value of [A] and [W] decrease. Behaviour at the threshold is very interesting. Just below the threshold spots form, but are lost in a process of slow decay that is 'bright', i.e. involves the production of a high concentration of W, which absorbs a lot of light. Above the threshold spots form but decay more and more slowly. This is associated with a lower rate of light absorption.

Since we are interested in systems capable of high rates of light absorption, we should measure the total amount of light absorption that took place on the surface over 20000 timesteps for each of the above parameter values,. This is proportional to the integral of W. These results are shown on the right of the table above. You can see that maximum light absorption over the course of the experiment occurs in the 0.0219 regime where stable existence of a spot depends on constant re-cycling of W into A and decay of A into W. In the regimes of repulsion above and below this, A either decays to zero, or persists in a non-dissipitive manner at high concentration. This is the spatial analogue of the optimal decay rate that was observed in the well-mixed case.

Spots that are too small decay in the 0.0219 regime above. It seems there is a characturistic spot size at which spots can be persistant at steady-state in the 0.0219 regime.

No phase seperation has been considered for the W and X particles. What physical properties of these particles might allow the maximization of light absorption rate?

Conclusions

What is the significance of these results? If it is the case that there is an autocatalytic particle existing in the 0.0219 regime, i.e. a regime where the extent to which A is capable of repelling R is precisely such that the rate of light absorption is maximized, then several properties suggest that this may have been an important stepping stone towards the origin of life.

0. If the existence of some other particle P requires its constant re-synthesis due to the fact that it is susceptable to decay, then it is only chemical organizations capable of either 1. reducing P's decay rate or 2. increasing P's rate of production, that can maintain a steady-state concentration of P. Let us imagine that a autocatalytic reaction whose constituents take part in many side-reactions (e.g. the formose reaction) may produce a large number of such molecules P1, P2... etc.... each with a completely random set of chemical and physical properties. We can add such particles to the spatial system one by one and see whether they are stable, or whether they destroy the system. All that we require is an underlying very simple autocatalytic system that can persist under a broard range of conditions, and whose probabilty of spontaneous formation is high given a set of environmental conditions available on the early Earth. Once a high probability simple autocatalytic system exists as a recycling system capable of existence at steady-state, then as long as it is capable of producing chemical variety, I propose that through mere random chemical novelty, it is possible to obtain increasing complexity leading to life. There is a variety of such simple autocatalytic systems to choose from, the spot model above is one such example. Some may more suitably allow the recursive generation of functional constraints. Examples of mechanisms for the persistence of P molecules are given below...

1. If W has certain properties such as for example, the capacity for membrane formation, then the individuality of the spot can be reinforced. How can we introduce into the model the W property that it forms membranes? This can be defined by specifying an anisotrophic potential that requires a high A concentration on one side and a high R concentration on the other side of a W molecule, but not all around the W molecule! Several other methods exist for defining anisotropic molecule behaviour.

2. If it is the case that A reacts not just with R to produce W, but also reacts to produce other molecules that might have some catalytic effect on the rate of A production then if good molecules that increase the extent of A production can be trapped within the spot, and bad molecules that reduce the extent of A production can be forced to leave the spot, then the spot can become more stable.

3. If it is the case that several recycling systems are required to maintain A, not just the recycling of W into A, then if W traps the other molecules that require recycling then the spot can become more stable.

4. So far there is no concept of HOW chemicals should change their properties. For now we are merely interested in the means by which high rates of recycling can be obtained, at increasing levels of complexity, e.g. increasingly complex physical properties posessed by chemicals, not how these mechanisms can arise.

From the above considerations, the following experiment suggests itself.

1. Introduce random new reactions that produce new particles one at a time onto the surface. Novel particles are produced in proportion to the extent to which A undergoes new reactions. New reactions are undertaken by A with other particles, or can be re-arrangements of A. New reactions are most likely to occur between or within particles that have high free energies of formation. A new reaction can basically be considered to be a random "mutation", i.e. the addition of a new random reaction to the chemical network. The unit of selection is the supramolecular assembly in which that new reaction occurs. Greater search through the chemical space of new reactions is undertaken by assemblies with high free energy molecules. Such assemblies must be capable of re-cycling, or steady state energy absorption. This is a natural fitness function, only the supramolecular organizations capable of maintaining high energy molecules undertake significant search in the chemical space and so reach regions of greater dynamic stability. Some 'mutations' will take the system to regions of static stability, i.e. not dissipitive stability, but these will no longer be capable of further change, other 'mutations will take the system towards regions of instability, where the organization is destroyed.

This scenario can be explored by applying the same GA approach used to optimize flux in the well mixed reactor. The algorithm is as follows.

0. Initialize a population of surfaces. Each surface has its own chemistry with initially 0.0219 region spots. In addition, each surface is initialized with a random new reaction. The new reaction is produced according to a generative regime in the same manner as previously, e.g. species with higher free energies of formation are more likely to undergo a bimolecular reaction. Any new species are assigned completely random physical properties! Some of these bimolecular reactions may absorb light, some may release heat. In the extreme case a 1:1 GA can be used, i.e. we just have one surface, we mutate it, if the mutant is betten than the parent, we replace the parent, if not, then we mutate the original parent.

1.Each fitness trial involves assessing the integral of light absorption in a spatial system (perhaps a 1D spatial system for the sake of spead, similar dynamics should be expected on a line as on a 2D surface).

The physical properties of species will also have to be stored in the file. These can be stored in a seperate file as a matrix. Ideally we should be able to visualize the end state of each fitness assessment, i.e. it should be stored in a file so that a pastiche of end-states can be viewed all at once and analysed in detail afterwards.The code should run in the command line so that it can be ported to Clusters and run remotely. However, during development it would be useful to be able to visualize the simulations in real time.

Experiments with the GA acting in the spatial system are discussed in the following page here.

 

Implementation of the Grey-Scott Model

The Grey-Scott reactions are as follows.

 

How to Model Attraction as Well as Repulsion?

 

How to Model Anisotropic Membranes?

 

How to Model Selective Permeability?

 

 

 

 

 

Return to Evolution of Metabolism Part 4: