Evolution of associative learning in chemical networks
Networks were initialised as follows. A small number of “seed” chemicals (by default, 4) with distinct formulas of length 3 were added to the network. New chemical species, whether generated at initialisation or due to adding a new reaction to the network during initialisation or mutation, were initialised with uniformly random parameters in the following ranges: potential [0–7.5], initial concentration [0–2], food inflow [0–1], decay [0–1]. The function to add a new reaction was called 20 times, thereby adding an unpredictable number of new chemicals to the network. New reactions, whether generated during initialisation or mutation, were initialised with a uniformly random “favoured reaction constant” in the range [0–0.1]. The input-output list for the network was shuffled fairly.
The networks were evolved using a non-generational genetic algorithm (GA) similar to the Microbial GA [41]. A genetic algorithm is the natural selection algorithm run in a computer [10], specifically it is artificial selection in which an explicit fitness function (phenotypic target) is defined, rather than allowing fitness to emerge as the result of ecological interactions. In our case the fitness function rewards chemical networks capable of the kind of associative learning we require. The basic algorithm was as follows:
Pick two different networks from the population (for spatial evolution, choose two neighbours)
Replace the worse-performing network with a mutated copy of the better-performing network
All reactions were modelled using reversible deterministic mass action kinetics (apart from the implicit decay reactions which are irreversible). It is clearest to explain this scheme by example.
The rate at which a reaction takes place, in our simulation, is set equal to the product of the concentrations of those species on its left-hand side, multiplied by its rate constant. The reaction consumes its reactants at this rate and generates its products at this rate. The overall rate of change of a species' concentration due to explicitly-modelled reactions is equal to the sum of the rates at which it is generated (over all reactions) minus the sum of the rates at which it is consumed (over all reactions). Spontaneous decay (at a rate for chemical X) contributes an additional term to this sum, and inflow (at rate for the food chemical F) contributes an additional term to F's rate of change. Hence, a system consisting only of the reaction
For computational efficiency, simulations during evolution used Euler integration with a step size of 0.01 time units. Input boluses were modelled as discontinuous jumps in concentration at the given time steps. These simulations were qualitatively validated after evolution using a variable step-size Runge-Kutta ODE solver.
Networks were simulated on chemical protocols, with each protocol consisting of a time series of input boluses, and a time series of target values for the network output. Note that for most time steps, the input bolus values were zero and the target output values were “don't care”. The exact details of the protocol inputs and targets varied from task to task.
For every task, networks were simulated on a number of protocols, and the (instantaneous) concentration of the designated network output chemical compared to the protocol target for every time step. The fitness of a network was set equal to the negative mean square difference between these two quantities averaged over all protocols and all time steps (ignoring time steps where a “don't care” target was specified). In order to provide a reliable fitness comparison, when two networks were chosen for competition during evolution, they were evaluated on the same set of protocols. Additionally, the protocols for different experimental conditions within the same task were deliberately matched to be similar, so that network response to the experimental condition could be measured as directly as possible.
Initial experiments indicated that randomly generating protocols during evolution results in very noisy fitness comparisons, with little fitness gradient for evolution to climb. To avoid this problem, for each task we generated fixed “training data” and saved it to file. Networks were evaluated during evolution on their performance on the training data set. For most tasks, the training data set was a file consisting of 10 randomly generated protocols. A number of tasks were devised requiring the detection of different environmental features by the networks. Some of these tasks were “clocked”, i.e. pulses were constrained to only occur at predetermined regular “clock tick” times, and some were not.
This task constrained B boluses to a regular “clock tick” schedule every 100 time steps and had two experimental conditions. There was only a 0.5 probability of a chemical B bolus on a given clock tick. In the “associated” condition, a chemical B bolus was always followed 20 time steps later by a chemical A bolus. In the “unassociated” condition, chemical A boluses never occurred. A single protocol featured both experimental conditions, with identical B boluses in each condition. The desired behaviour for the network was: upon receiving a pulse of chemical B, output either zero (in the “unassociated” condition) or one (in the “associated” condition) for 20 time steps afterwards.
This was identical to the previously described task except that there was a small (p = 0.1) probability of “noise” occurring at each time step with a chemical B bolus. Noise consisted of a B bolus being followed by an A bolus in the “unassociated” condition or a B bolus followed by no A bolus in the “associated” condition. Within a single protocol, the occurrence of noise was matched between experimental conditions.
This task had two experimental conditions and involved boluses at random intervals. In both conditions, pulses of chemical B occurred at random intervals uniformly in the range [100, 300]. In the first (“associated”) condition, a pulse of chemical B was followed shortly afterwards (20 time steps) by a pulse of chemical A. In the second (“unassociated”) condition, pulses of chemical A occurred independently of B, at random intervals uniformly in the range [100, 300]. Within a single protocol, pulses of chemical B were identical.
This task, featuring two experimental conditions, was specifically designed to involve a non-trivial accumulation of information. Within this task, input “events” occurred randomly at a low rate (0.025 per time step) with a refractory period of 50 time steps between events, over a total period of 2000 time steps. Each event consisted of either a pulse of chemical A followed closely (20 time steps later) by a pulse of chemical B, or vice versa. In the first experimental condition (“A→B”), events were 75% likely to be “A→B” pulses and 25% likely to be “B→A” pulses, and vice versa for the second (“B→A”) experimental condition. The desired output behaviour was to respond to a “B” pulse with a low output in “A→B” environments and a high output in “B→A” environments. Note that this task was both noisier than the other tasks and involved a longer evaluation period (to allow the noise some time to average out).
Unlike the other tasks, every environment in this task was designed by hand. The intention was to construct a range of radically different environments such that both short- and medium- term network memory-traces would be required to attain maximum fitness. The inspiration was loosely drawn from the concept of the “radical envelope of noise” [Jakobi, 1998]. Input pulses (boluses) in this task always occurred in closely-separated pairs, although the second bolus in a pair did not have to contain the same chemical as the first bolus. The pulse pairs occurred at regular intervals of 100 time units each. Each experimental condition was characterised by a “typical” pulse pair (A→A, A→B, B→A or B→B). In addition to the “typical” pulse pair corresponding to the experimental condition, every protocol for this task also had a “noise” pulse pair. There were in total 4 protocols (one for each pulse pair type), each containing 4 experimental conditions, for a total of 16 different input series. A single input series had the following structure: First, a pulse pair corresponding to the protocol's “noise” pair. Next, three “signal” pulse pairs all of the “typical” type for that experimental condition. Next, a “probe” pulse pair (see below). Next, another “noise” pulse pair of the protocol's “noise” type. Last, a final “probe” pulse pair. “Probe” pulse pairs consisted of a pulse of “B” chemical followed by either a pulse of “A” chemical (in the B→A environment) or a pulse of “B” chemical (in other environments). The desired network behaviour was to produce a low output for 10 time steps prior to each “probe” pulse pair, followed by either a high output (in the B→A environment) or a low output (in other environments) for 20 time steps. Errors in the B→A environment were weighted three times as heavily as errors in the three other environments.
We calculate the number of reactions per effective chemical species in a network by first excluding any species which do not take part in reactions (this is possible if all reactions featuring a particular species are lost from a network by structural mutation). We then simply calculate the mean number of distinct reactions each remaining species is involved in.
To investigate the effects of different genetic encoding factors on network connection density, we conducted 10 evolutionary runs on the 2-bit environment problem in each of 4 encoding variations. These were:
A benchmark case with maximum formula length 4, 2 symbols in the chemical alphabet, and the aggregation chemistry.
For all these runs, we recorded the effect of every mutation on both fitness and also the number of reactions per chemical species.
Our method is as follows. We imagine an ideal Bayesian reasoner, equipped with knowledge of the statistics of the different network task environments. For each input train, at each point in time, we calculate what subjective probability the reasoner should assign to the possibility that the input train up to that point came from an “associated” environment. This establishes what the ideal Bayesian posterior would be at each point in time for each input train. If a network's chemical concentrations somehow encode this time-varying Bayesian posterior in all environments, then it would seem reasonable to attribute a Bayesian interpretation to the network. For the purposes of this paper, we will skirt over the complexities introduced by the non-dissipation of information in smooth continuous dynamical systems. In principle, the state of our simulated networks will usually contain all information about their historical inputs, because information can be stored in arbitrarily small differences in concentrations. However, in practice this information will be destroyed by noise.
Calculation of the ideal posteriors for our environments is straightforward. A random variable X will represent the type of environment: either 1 (“associated”) or 0 (“unassociated”). Another random variable Y(t) will represent the train of input boluses up to time t. The ideal posterior probability of being in the “associated” environment after observing an input train y is
The prior P(X = 1) was set equal to the proportion of “associated” environments in the network training sets, i.e. 0.5. The probabilities P(Y(t) = y|X) were calculated as follows. In the environments we analysed in this way, Input trains were organised into “events”: a bolus of one or other input chemical, followed possibly at a set short interval by another bolus. The time between input events always exceeded the inter-spike interval within an event. The timing of events in an input train provided no information about the type of environment. Hence, events can be extracted from an input train and treated as a discrete process. “Associated” and “unassociated” environments correspond to Bernoulli processes with “associated” (B→A) and “unassociated” (B alone, for the noisy clocked task, or A→B for the AB-BA task) events. Thus, the posterior P(X|Y(t)) can be calculated by counting the total number of associated and unassociated events in Y(t).
where n and m are the number of “associated” and “unassociated” events in Y(t), and p1 and p0 are the probabilities of an “associated” event in the “associated” and “unassociated” environments respectively. This gives
which can be attached to a time series at the appropriate points after events have occurred.
We use a straightforward logistic regression model to match network concentrations to Bayesian posteriors. Given a concentration vector x, a weight vector w and a bias value b, the model is
where is the scalar product of x and w. Note that the output of this function is bounded between 0 and 1 like a probability value (this would not be the case for a linear regression model). The idea is that we can investigate the degree to which a rational Bayesian belief is encoded transparently in the network's state. We expect that at time t, having observed an input history Y(t),
To determine appropriate model parameters, we randomly generate 200 environments (100 “associated” and 100 “unassociated”), and run the evolved network in those environments. Weights w and bias b are set to minimise mean square error over all environments and time steps, using Levenburg-Marquadt optimisation. No attempt was made to regularise the parameters or otherwise avoid overfitting, since the model has relatively few parameters. For comparison, 200 random networks (produced by random initialisation followed by 200 mutations) were tested in the same way.
The work is funded by an FP-7 FET OPEN Grant, E-FLUX, the Hungarian NAP project, and the Templeton Grant FQEB entitled “Bayes, Darwin and Hebb”. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.