A thermodynamically consistent chemical spiking neuron capable of autonomous Hebbian learning
AA thermodynamically consistent chemical spiking neuroncapable of autonomous Hebbian learning
Jakub Fil and Dominique ChuCEMS, School of Computing, University of Kent, CT2 7NF, Canterbury, UK
Abstract
We propose a fully autonomous, thermodynamically consistent set of chemical reactions that implements a spiking neuron.This chemical neuron is able to learn input patterns in a Hebbian fashion. The system is scalable to arbitrarily manyinput channels. We demonstrate its performance in learning frequency biases in the input as well as correlations betweendifferent input channels. Efficient computation of time-correlations requires a highly non-linear activation function. Theresource requirements of a non-linear activation function are discussed. In addition to the thermodynamically consistentmodel of the CN, we also propose a biologically plausible version that could be engineered in a synthetic biology context.
Beside electronic devices there are a number of substrates that can be used to realise computations [Adamatzkyet al., 2019]. A particularly important class of what is sometimes called “unconventional computers” arechemical systems and especially biochemical systems [Amos, 2004]. These underpin the remarkable abilities ofunicellular organisms to adapt to changes in their environment in the absence of a nervous system. Examplesof biochemical information processors include sensing [Alon, 2019, Govern and ten Wolde, 2014a,b], chemotaxis[Hoffer et al., 2001, Yi et al., 2000] or diauxic growth [Chu, 2017, 2018, Chu and Barnes, 2016]. Understandingthe function of biochemical systems has so far been the domain of molecular biology. With the advent ofsynthetic biology it is now possible to build custom-made biochemical systems and hence to engineer wet-waremolecular computers. There are a number of possible applications of such computers. Particularly promisingareas include personalised precision medicine, especially targeted drug delivery but also environmental clean-upor sensing. Before any such computers can be built, it is first necessary to understand how to program them.Here we propose a chemical reaction network — henceforth referred to as the chemical neuron (CN) —that behaves like an artificial spiking neuron (SN). SNs are a type of learning machine that is used widely inArtificial Intelligence as a component of neural networks [Afshar et al., 2020, Sengupta et al., 2019]. It is nowalso well established that a single SN has significant learning capabilities [Fil and Chu, 2020, G¨utig, 2016]. Incomputer science they have been applied to a number of real-world tasks including principal component analysis[Oja, 1982], recognition of handwriting [Diehl and Cook, 2015] or classification of fighter-planes [Afshar et al.,2019]. There are many different models of SNs in the literature. Commonly a SN has an internal state whichis typically represented by a positive number. The input state may decay, that means that it reduces over timewith some rate. The input can also increase when the SN receives a stimulus (an input spike ) via one of its N input channels. Importantly, these input channels are weighted such that the increase of the internal statevariable is different for the different input channels. The higher the weight of a channel, the stronger a stimuluswill affect the internal state.A central problem of research in neural networks is to find ways to to determine the weights of a neuron suchthat it performs a particular task. The process of finding these weights is usually referred to as the “training”or alternatively as the learning phase. A common method to train SNs is Hebbian learning . The basic idea ofHebbian learning is that the weights relating the input to the output are plastic, i.e. they change over time.In Hebbian learning, the weights will be strengthen by some amount if the input to the neuron coincides withthe neuron generating an output event. A well-known special case of Hebbian learning is associative learning.Assume that A and A are both inputs to B. Assume further that the connection between A is sufficientlystrong such that its firing can trigger a firing event of B. If A fires usually at around the same time as A , thenits weights will be strengthened. Eventually, the rates of the second channel will have adjusted sufficiently suchthat firing of A alone will be sufficient to trigger a response on its own.There have been a number of previous attempts to design chemical implementations of learning machines.Amongst the earliest is perhaps work in the 1980s by Okamoto and collaborators [Okamoto et al., 1988] whoshowed that certain biochemical systems implement the McCulloch-Pitts neuronic equations. A chemical neuronhas been implemented by [Hjelmfelt et al., 1991]; however, this neuron could only detect stable concentrationsof input molecules and no learning was shown. Fernando and coworkers [Fernando et al., 2009] showed that1 a r X i v : . [ c s . N E ] S e p Results Fig. 1:
Example simulation showing the core of the CN dynamics. The graphs show the internal state B , thelearning signal E and the weight H for a single channel. We assume a bolus provided at time t = 0 . t = 0 .
03 and consequently the weight is increase by (in this case) 15 molecules of H .a relatively simple set of molecular interactions can implement associative learning. Their model is fully au-tonomous, but it is also inflexible. Association is learned after just a single coincidence, and hence not fit todetect statistical correlations robustly. Moreover, the system cannot forget the association between the inputs.[McGregor et al., 2012] improved on that with systems that were found by evolutionary processes. A biochem-ically plausible implementation of the associative learning based on genetic regulatory networks was proposedby [Macia and Sole, 2014, Macia et al., 2017a,b]. They imagine their system not to be within a single cell, butdistributed across multiple compartments. This has the advantage that is reduces crosstalk between molecularspecies and allows more complicated networks to be implemented. Another model close to biology is by Nesbeth[Nesbeth et al., 2016]. A perceptron formulated in an artificial chemistry was proposed by Banda et al. [Bandaet al., 2014]; later this was extended to a fully fledged feed-forward neural network [Blount et al., 2017] whichcould solve the XOR problem. Their model still requires regular intervention by outside operators. Besidesthese simulation studies, there have also been attempts to implement learning in vivo [Chen and Xu, 2015,Nesbeth et al., 2016, Racovita and Jaramillo, 2020, Shirakawa and Sato, 2013].The main contribution of this article is a design for a novel chemical neuron (CN), that is a family ofmolecular reaction networks. The CN implements a fully autonomous spiking neuron capable of full Hebbianlearning. A tuneable parameter of the design of the CN is the non-linearity of the activation function . Thisactivation function determines the probability that the CN generates a learning signal, depending on the internalstate.We demonstrate that the CN is capable of Hebbian learning, in the sense that its internal molecular abun-dances will reflect statistical biases of its inputs. More specifically, via each its input channels the CN can accept boli of some chemical species. Over time these boli will lead to a reconfiguration of the steady state abundancesof the CN. These abundances will then be reflective of statistical biases of the input boli. In particular weconsider two types of biases: ( i ) Frequency biases (FB), where the average waiting time between two boli isdifferent for different channels. Here the input times are independent but not identically distributed randomvariables. ( ii ) Time Correlations (TC), where the frequency of boli may be the same for each input channel, butinputs for one channel tend to happen before or after the inputs of one or more other channels. Here, input boliare identically distributed, but not independently drawn from one another. We find that the CN can identifyboth types of biases. However, we also show that the FB task is trivial for chemical systems, whereas the TCtask can only be solved when the non-linearity is high. We find that a high non-linearity comes with minimalresource requirements.In this contribution, we will give two versions of the CN. The first (basic) version will be a minimal set ofreactions. It is thermodynamically consistent in that it consists only of micro-reversible reactions with mass-action kinetics. This basic version cannot be realised easily though. Therefore, we shall propose a second versionof the model, which is not thermodynamically explicit, but biologically plausible and can likely be engineeredin a synthetic biology framework.
Overview
We implement the CN as a set of micro-reversible elementary chemical reactions obeying mass-acting kinetics. Micro-reversibility makes the model thermodynamically consistent. See table 1 for the list ofreactions. The system is best understood by thinking of the molecular species A i as the input to the system via Results Input A n k AB −−− (cid:42)(cid:41) −−− k BA BActivation function B + E i k + −− (cid:42)(cid:41) −− k – E i +1 , i < m −
1B + E m –1 k + −−− (cid:42)(cid:41) −−− k –last E Learning A n + E k AE −−− (cid:42)(cid:41) −−− k EA A E n k EH −−− (cid:42)(cid:41) −−− k HE H n + E A n + H n k AH −−− (cid:42)(cid:41) −−− k HA AH n k HB −−− (cid:42)(cid:41) −−− k BH B + H n Leak H n k H ∅ −−−→ ∅ B k B ∅ −−−→ ∅ Tab. 1:
List of chemical reactions in a single CN.input channel i . The weight equivalent of the i -th input channel of the CN are the abundances of the species H i . The species E is the activated form of E and takes the dual role as the learning signal and the output ofthe CN. The internal state of the CN is represented by the abundance of the molecular species B . Input-Output relation
We assume here that the CN has N different species of input molecules A , . . . , A N .These represent the N input channels and are associated with the corresponding weights H , . . . , H N . Theinput is always provided as an exponentially decaying bolus at a particular time t si . Concretely, this means thattime t = t si the CN is brought in contact with a reservoir consisting of β (unmodelled) precursor molecules thatthen decay into A i molecules with a rate constant κ >
0. A particular consequence of this is that the A i arenot added instantaneously, but will enter the system over a certain time. This particular procedure is a modelchoice that has been made for convenience. Different choices are possible and would not impact on the resultsto be presented. The important point is that there is a sense in which inputs can be provided to particularinput channels i by providing a bolus of A i at a particular time t . Finally, note that the CN has a dissipativecomponent. Both B and H i molecules decay with a fixed rate constant. This enables the system to reach asteady state provided that the input is stationary.The basic idea of the CN is that input boli A i are converted into internal state molecules B . The speed ofconversion depends on the amount of weight H i . If at any one time there is enough of B in the system thenthe learning signal E is created by activating E molecules. Once the learning signal is present then some of the A i are converted into weight molecules, such that the weight of the particular input channel increases. Thisrealises Hebbian learning in the sense that the coincidence of inputs B and output E should activate learningfollowing the well known Hebbian tenet “What fires together, wires together”. Activation function
Precisely what constitutes “enough” B molecules to trigger a learning signal dependson the precise tuning of the system and will be discussed extensively below. The link between the internalstate molecules B and the learning signal is often called the activation function . More specifically, this is thefunctional dependence of the probability to observe the activated form E and the abundance of B . As we willshow below, this activation function determines the ability of the CN to learn various kinds of input biases. Inthe simplest, albeit hard to realise case, the activation is a step function. This means that the learning signalis only generated if the abundance of B molecules is greater than some threshold ϑ .In reality, the activation function will at best approximate a step function. In the CN it is realised as follows:In all simulations presented below, we start at time t = 0 with a fixed number of E molecules. Each of the E molecule has m binding sites to which the internal state molecules B can bind. Once all m binding sitesare occupied, then E is converted into its active form E . We make here the simplifying assumption that theconversion from E to E is instantaneous once the last B binds. Similarly, if a B unbinds, then the E changesimmediately to E . In this model, the balance between E and E molecules depends on the binding and unbindingrates of B . We assume that there is a cooperative interaction between the B molecules such that unbinding of B from E is much slower than unbinding from E . With an appropriate choice of rate constants, this systemis known to implement ultrasensitivity, i.e. the probability for the fully occupied form of the ligand chain ( E )to exist transitions rapidly from close to 0 to close to 1 as the concentration of ligands approaches a thresholdvalue ϑ ≈ k + /k − . Such systems are often approximated by the so-called Hill kinetics. It can be shown thatthe maximal Hill exponent that can be achieved by such a system is m [Chu et al., 2009]. This means that thechain-length m , which we henceforth shall refer to as the “nonlinearity”, controls the steepness of the activationfunction of E . In the limiting case of m = ∞ , this will be a step function, whereby the probability to observe E is 0 if the abundance of B is below a threshold and 1 otherwise. We are limited here to finite values of m > m = 1. The parameter m and Results A A B H Weight 1Weight 2
Fig. 2:
Associative learning in CN. The first two graphs show the inputs A and A . Clearly, a single A doesnot lead to a sufficient increase of the internal state B , such that no learning signal is triggered. After afew coincidences of A and A the weights H (last graph) have increased sufficiently for A to triggera signal in its own at time t = 0 .
8. Note the increase in weights for the second channel after eachcoincidence. On the contrary, the weight associated with the first input channel experiences a marginalincrease. This shows that the system implements a form of local weight normalisation mechanism similarto the one proposed by Oja [1982].hence the steepness of the activation function will turn out to be a crucial factor determining the computationalproperties of the CN.The functional dependence of E on B is the analogue of an activation function in SNs. However, note thatunlike standard practice in SNs, the activation function determines the probability of a learning signal beingpresent as a function of the abundance of B molecules, i.e. it is a stochastic function. Learning
The weights H i of a particular input channel i can only be increased if two conditions are fulfilled:( i ) the learning signal E is present and ( ii ) there are still input molecules A i in the system. In short, learningcan only happen if input and output coincide, which is precisely the idea of Hebbian learning. Throughout this paper we will assume that the dynamics of
A, B and E is fast compared to the change inconcentration of H . This is a crucial assumption to allow the weights to capture long-term properties of inputs;in particular, the weights should not be influenced by high frequency noise present in the system. Furthermore,we also assume that the lifetime of E is short. For details of the parameters used see table 3.We first demonstrate that the CN is capable of associative learning; see fig. 2. To do this, we generate a CNwith N = 2 input channels. Then, we initialise the CN with a high weight for the first channel ( H = 100) anda low weight for the second channel ( H = 0). Furthermore, we set the parameters of the model such that abolus of A is sufficient to trigger an output, but a bolus of A , corresponding to stimulating the second channelis not; see fig. 2 for details. This also means that presenting simultaneously both A and A triggers a learningsignal and increases H and H . If A and A coincide a few times, then the weights of A have increasedsufficiently so that a bolus of A can push the internal state of the system over the threshold on its own; seefig. 2. This demonstrates associative learning. Note that unlike some previous molecular implementations ofassociative learning (e.g. [Fernando et al., 2009]), the system requires several coincidences before it learns theassociation. It is thus robust against noise. It will also readily unlearn the correlation if input patterns change. We now show that the ability of the CN to learn extends to full Hebbian learning with an arbitrary number of N input channels. First we consider the FB task. To do this, we provide random boli to each of the N inputchannels. Random here means that the waiting time between two successive boli of A i is distributed accordingto an exponential distribution with parameter 1 /f i , where f i is the frequency of the input boli to channel i .The CN should then detect the difference in frequencies f i between input channels. We consider the FB taskas solved if (after a transient period) the ordering of the abundances of weights reflects the input frequencies,i.e. the number of H i should be higher than the number of H j if f i > f j . Results TC 2 TC 3 TC 4 FB 2 FB 3 FB 40.000.050.100.150.200.250.300.35 N o r m a li s e d w e i g h t s a t s t e a d y s t a t e H1 H2 H3 H4 H5
Fig. 3:
Normalised weights for a variety of TC and FB tasks. The first (blue) bar refers to the first weight,the second (orange) to the weight for the second channel and so on. Each value represents the averageover 300 time units of a single simulation. Data was only collected after the weights reached the steadystate (after 700 time units). In all experiments we set the number of E molecules at the start of thesimulation to 40. The nonlinearity was set to m = 5 for the TC, and m = 1 for FB.In order to test this, we consider 3 variants of the FB task: We use a CN with N = 5 and m = 1 andassume that boli to the first two channels input come at a frequency of 4Hz whereas channels 3, 4 and 5 fireat a frequency of 2 Hz; we call this variant FB 2. Similarly, for FB 3 and FB 4 the first 3 and 4 channelsrespectively fire at the high frequency. Fig. 3 shows the steady state weights for each of the three tasks. Asexpected, in each of the experiments the weights of the high-frequency inputs are higher when compared to thelow frequency inputs. We conclude, that the CN can work as a frequency detector.A more interesting scenario is the direct generalisation of the associative learning task to an arbitrarynumber of input channels. For this problem we assume that all input frequencies are the same, i.e. f i = f j for all i, j ≤ N . Instead of differences in frequency, we allow temporal correlations between input boli of somechannels. This means that there are pairs of channels i, j and time windows ∆ τ such that the probability toobserve a bolus of j between time t and t + ∆ τ — where t is the time where input A i fired — is greater thanthe probability to observe a bolus between t − ∆ τ and time t . In practice, we implement such correlations asfollows: If A and A are temporally correlated then each bolus of A is followed by a bolus of A after a timeperiod of δ + ξ ; with δ being a fixed number and ξ is a random variable drawn from normal distribution with µ = 0 and σ = 0 . A i tends to precede A j , then theabundance of weight H i should be higher than the abundance of H j . Furthermore, if A i is correlated with someother channel k , but A j is not, then the abundance of H i must be greater than that of H j . In order to testthe ability of the system to detect TC biases, we initialised the CN with N = 5 input channels. Initially theweight molecules H i = 0. We then determined the steady state weights. We compared four different scenarios:( i ) there are correlations between A and A (TC 2), ( ii ) A , A , A (TC 3), and ( iii ) A , A , A , A (TC 4).The temporal order is always in ascending order of the index, such that in the last example, A occurs before A , which in turn occurs before A . We find that the behaviour of the CN is as expected; see fig. 3. At steadystates the weights reflect the correlation between input channels, including the temporal ordering. We now examine how the nonlinearity parameter m impacts on the ability of the system to compute. Weconsider two extreme cases: Firstly, the case of minimal non-linearity (i.e. m = 1) and the secondly, thelimiting case of maximal nonlinearity (i.e. m = ∞ ). This latter case would correspond to an activation functionthat is a step function. Strictly speaking, the CN cannot realise a pure step function, but it still providesvaluable insight to think through this limiting case.We consider first this latter scenario with a CN with two inputs A and A . In this case, there will bea learning signal E in the CN if the abundance of B crosses the threshold ϑ . Let us now assume that theparameters are set such that a single bolus of either A or A is not sufficient to push the abundance of B overthe threshold, but a coincidence of both is: Results • A single bolus of A will not lead to a threshold crossing. No learning signal is generated and weights arenot increased. • If a bolus of A coincides with a bolus of A then this may lead to a crossing of the threshold of theinternal state. A learning signal is generated. Weights for both input channels 1 and 2 are increased(although typically not by equal amounts).Next consider an activation function tuned to the opposite extreme, i.e. m = 1. It will still be true that both A and A are required to push the abundance of B across the threshold. However, the learning behaviour ofthe CN will be different: • A single bolus of A will not lead to a threshold crossing. A learning signal may still be generated evenbelow the threshold because the activation function is not a strict step function. The weight H willincrease by some amount, depending on the bolus size. • If a bolus of A coincides with a bolus of A then this will lead to more learning signal being generatedthan in the case of A only. As a result, the weights for both input channel 1 and 2 are increased by morethan if they had occurred separately.These two extreme cases illustrate how the CN integrates over input. In the case of low non-linearity theweights of a channel will be a weighted sum over all input events of this channel. The weights will be higher forchannels whose boli coincide often. On the other hand, a step-like activation function will integrate only overthose events where the threshold was crossed, thus specifically detect coincidences. From this we can derive twoconjectures: • The higher the nonlinearity, the better the CN at detecting coincidences. Low non-linearity still allowscoincidence detection, but in a much weaker form. • As the bolus size increases, the CN will lose its ability to detect coincidences, especially when the bolussize is so large that a single bolus is sufficient to push the abundance of B over the threshold. In thiscase, a single input spike can saturate the activation function, thus undermining the ability of the systemto detect coincidences effectively.In order to check these conjectures, we simulated a CN with 3 inputs, where A and A are correlated and A fires at twice the frequency of A and A . We considered the minimally nonlinear case ( m = 1) and a moderatenon-linearity ( m = 4); Fig. 4 shows the weights as a function of the bolus size. The minimal non-linear CNdetects both coincidences and frequency differences, but loses its ability to detect coincidences as the bolus sizeincreases. This is consistent with the above formulated hypothesis. In contrast, for the non-linear CN andmoderately low bolus-sizes the weights indicate the coincidences strongly (i.e. the weights H are highest). Yet,as the bolus size increases the non-linear CN also loses its ability to detect coincidences.Next, we check how the coincidence detection depends on the time-delay between the correlated signals.To do this we created a scenario where we provided two boli to the system. The first bolus A comes at afixed time and the second one a fixed time period δ thereafter. We then vary the length of δ and record theaccumulation of weights H as a fraction of total weight accumulation. Fig. 5 confirms that the CN with thelow non-linearity is less sensitive to short coincidences, but can detect coincidences over a longer period of time.In contrast, as the non-linearity increases, the differential weight update becomes more specific, but also morelimited in its ability to detect coincidences that are far apart. In the particular case, for a δ > . m > m = 1 shows some differential weightupdate throughout.Next, we tested the conjecture that the TC can be solved more effectively by the CN when the nonlinearityis higher. To do this, we generated a CN with N = 5 input channels on the TC 2 task. We then trainedthe CN for nonlinearities m = 1 , . . . ,
10. As a measure of the ability of the system to distinguish the weights,we used the index of dispersion, i.e. the standard deviation divided by the mean of the weights. A higherindex of dispersion indicates more heterogeneity of the weights and hence a better ability of the system todiscriminate between the different frequencies. Consistently with our hypothesis we found that the ability todiscriminate between frequencies increases with the nonlinearity. However, it does so only up to a point (theoptimal nonlinearity), beyond which the index of dispersion reduces again; see fig. 6. Increasing the bolus size,i.e. increasing the number of A i that are contained within a single bolus, shifts the optimal nonlinearity to theright. This suggests that the decline in the performance of the CN for higher chain lengths is due to a resourcestarvation. The realisation of the sigmoidal function, i.e. the thresholding reactions in table 1, withdraws m molecules of B from the system. As a consequence, the CN is no longer able to represent its internal stateefficiently and the activation function is distorted. If the total abundance of B is high compared to E , then thiseffect is negligible. We interpret this as a resource cost of computing non-linearity. The higher m , the higherthe bolus size required to faithfully implement the activation function. Results N o r m a li s e d w e i g h t s a t s t e a d y s t a t e m=1 H1 H2 H30.0 2.5 5.0 7.5 10.00.10.20.30.40.50.6 m=4
Bolus size
Fig. 4:
The weights as a function of bolus size for a CN with 3 inputs. Input A (green) is provided at 4Hz, A and A are correlated with a δ = 0 . m = 2 thesystem detects the higher frequency of A as indicated by its high weight throughout. The weight of A (orange) is only slightly higher than the weight of A , indicating that the CN detects the coincidence onlyto some limited extent. For increased non-linearity, the CN indicates the correlation because weights areonly increased when coincidences occur. If the bolus size is increased, then a single bolus is sufficient tocross the threshold and the frequency bias is recognised again. Δ H Δ H Δ Δ H m=1m=2m=3m=4 0.00 0.05 0.10 0.15 0.20δ0.40.50.60.70.8 Δ H Δ H Δ Δ H k B∅ = 0.05k B∅ = 0.1k B∅ = 0.15 Fig. 5:
Differential weight increase for different chain lengths. A value of 1 means that only the first inputchannel received weight accumulation. A value of 0 . B (right). The faster the removal, themore specific the coincidence detection, i.e. inputs need to occur within a narrower window. For bothgraph points were computed as follows: We simulated a CN with two input channels only and an initialcondition of H , H = 0. At time t = 0 we provided a bolus of A and after a variable time we providedthe bolus A . We then continued the simulation for another 0 . y -axis records therelative increase of H averaged over 1000 repetitions. Results I n d e x o f d i s p e r s i o n β/ϑ = 0.5 β/ϑ = 1.0 β/ϑ = 1.52 4 6 8non-linearity (m)050100150 Fig. 6:
Index of dispersion for different bolus sizes β expressed as a fraction of the threshold ϑ . We show thatfor TC 2 (left) and FB 2 (right) tasks. The index of dispersion measures how different the steady stateweights are from one another, and hence indicates how well the CN distinguished between input channels.Completely unbiased input would give an index of dispersion of ≈
0. The graph shows that for the TCtask, there is an optimal nonlinearity. Increasing the bolus size, increases the optimal non-linearity,which is consistent with the fact that the optimum is due to resource starvation.
We have established that the CN is a coincidence detector, which is why it is able to solve the TC task. Weshowed above that the CN can also perform the FB problem, i.e. it can record frequencies of input signals. TheFB task is fundamentally about integrating over input, which can be done naturally in chemical systems. Indeedit can be done by systems that are much simpler than the CN, for example: A i d −−→ ∅ . For appropriatelychosen values of d , the steady state value of A i would then reflect the input frequency. To understand this, notethat the input frequency determines the rate of increase of A i . This rate, divided by the decay rate constant d then determines the steady state abundance of A i , such that A i trivially records its own frequency. This systemis the minimal and ideal frequency detector.The CN itself is not an ideal frequency detector because all weight updates are mediated by the internalstate B . Hence, the weights are always convolutions over all inputs. The weights thus reflect both frequencybias and temporal correlations. In many applications this may be desired, but sometimes it may not be. Wenow consider the conditions necessary to turn the CN into a pure frequency detector, i.e a system that indicatesonly FB, but not TC. One possibility is to set the parameters such that the CN approximates the minimalsystem. This could be achieved by setting k BA (cid:28) k AB and all other rate constants very high in comparison to k AB . The second possibility is to tune the CN such that a single bolus saturates the threshold. In this case,the strength of the learning signal does not depend on the number of boli that are active at any one time. Asingle bolus will trigger the maximal learning signal. This is confirmed by fig. 4, which shows that as the bolussize increases, the system becomes increasingly unable to detect temporal correlations, but remains sensitive tofrequency differences. The basic model of the CN, as presented in table 1 is thermodynamically plausible and has the benefit of beingeasy to simulate and analyse. However, it is biologically implausible. As written in table 1 the molecular species A i , H i and B would have to be interpreted as conformations of the same molecule with different energy levels.Additionally, we require that these different conformations have specific enzymatic properties. Molecules withthe required properties are not known currently, and it is unlikely that they will be discovered or engineered inthe near future.As we will show now it is possible to re-interpret the basic model of the CN (table 1) so as to get a modelwhose elements are easily recognisable by computational biologists as common mechanistic motifs. This requiresonly relatively minor adjustments of the reactions themselves, but a fundamental re-interpretation of what thereactions mean. For the list of modified reactions see table 2 and fig. 7 for a graphical illustration of theintuition behind the model.The main difference between the basic CN and the biological version is that the latter is compartmentalised.While in the basic model the index (as in A i and H i ) referred to different species that exist in the same volume, Discussion
21 4 B o l u s i n j e c t i on ( A ) Lea k ( B ) h hh h B A*HA
Fig. 7:
Graphical representation of a compartment-version of the CN. A i and A j are the same molecular speciesbut contained in different compartments i and j respectively. We allow for an activated form of A ,denoted by A ∗ , which binds to the promotor site of h and activates its expression. H is an activetransporter molecules for A . The internal state molecule B is any of the A i when in the outermostcompartment. We assume that each compartment has a trans-membrane protein E with m extra-cellular binding sites. If all m binding sites are occupied by B , then the internal site becomes active(indicated by green) and can catalyse the activation of A .it should now be interpreted as indicating different compartments that separate the same molecular species.These compartments, which are themselves enveloped in a further compartment (the “extra-cellular space”),could be thought of as individual bacterial cells or else artificially generated membranes with a minimal genome.Input is provided by boli of the molecular species A into the compartment i . There is now also an activatedform of A , denoted by A ∗ . The conversion from A to A ∗ is catalysed by the learning signal E . Also new isthat each compartment contains a gene h that codes for the molecule H (we suppress the index indicating thecompartment). Expression of the gene is activated by A ∗ binding to the promoter site of h . We also allow a lowleak expression by the unactivated gene (denoted as h in table 2). Gene activation of this type is frequentlymodelled using Michaelis-Menten kinetics, thus reproducing in good approximation the corresponding enzymekinetics in the basic model of the CN. The molecules of type H are now transporters for A . We then interpretthe conversion of A i to B as export of A from compartment i to the extra-cellular space. In this interpretationthe molecular species B is then the same as A but contained directly in the outer compartment. The rate ofexport of A is specific to each compartment in that it depends on the abundance of H in this compartment.Finally, we interpret the E molecules as transmembrane proteins that are embedded in the membrane of eachcompartment. Their extra-cellular site has m binding sites for B molecules which bind cooperatively. When allsites are occupied then the intra-cellular part is activated, i.e. becomes E . In its activated form it can convert A to A ∗ .Another difference between the two versions of the models is that the molecule E is now specific to eachmembrane. The minimum number of copies of E is thus N whereas in the basic model a single copy of E at time t = 0 could be sufficient. This has two consequences. Firstly, at any particular time the number of occupiedbinding sites will typically be different across the different N compartments. This is a source of additionalvariability. Moreover, since the number of copies of E is higher than in the basic CN, the model in table 2 ismore susceptible to starvation of B as a result of the extra-cellular binding sites withdrawing molecules fromthe outer compartment. Both of these potential problems can be overcome by tuning the model such that theabundance of B molecules is high in comparison to E molecules.This highlights that the difference between the basic CN and its biological interpretation are deeper than thelist of reaction suggests. We checked that the biological interpretation still allows the same phenomenologicalbehaviour, provided that the parameters are adjusted. Fig. 8 confirms that both associative learning and fullHebbian learning (with 5 input channels) is possible with this model. To our knowledge the CN presented here is the first fully autonomous design for a chemical systems capable offull Hebbian learning. The model is, at least in principle, fully scalable. Previous attempts (e.g. [Fernando et al.,
Discussion A A B H Weight 1Weight 2
TC 2 TC 3 TC 4 FB 2 FB 3 FB 40.000.050.100.150.200.25 N o r m a li s e d w e i g h t s a t s t e a d y s t a t e H1 H2 H3 H4 H5
Fig. 8:
Associative (top) and Hebbian learning (bottom) for the biologically plausible version of the CN. Theexperimental setup is as in fig. 2, and fig. 3 respectively. For the parameters used see table 4. Bothexperiments approximated the ligand dynamics by a Hill function, in order to speed up the simulations.Input I n k IA −− (cid:42)(cid:41) −− k AI AActivation function B + E i k + −− (cid:42)(cid:41) −− k – E i +1 , i < m −
1B + E m –1 k + −−− (cid:42)(cid:41) −−− k –last E Learning E + A k AE −−− (cid:42)(cid:41) −−− k EA E A k A*E −−− (cid:42)(cid:41) −−− k EA* E + A * A * k A*A −−− (cid:42)(cid:41) −−− k AA* A h + A ∗ A · h −−− (cid:42)(cid:41) −−− hA · h h k leak −−−→ H n + h h k h −−→ H n + hA + H k AH −−− (cid:42)(cid:41) −−− k HA AH k HB −−− (cid:42)(cid:41) −−− k BH B + HLeak H k H ∅ −−−→ ∅ B k B ∅ −−−→ ∅ Tab. 2:
List of chemical reactions in a single CN unit interpreted as a cell. Molecular species
A, E, E , h , h and H are compartmentalised. Each compartment has a gene h which when activated by A ∗ can expressa transporter H . Discussion leaky integrate and fire (LIF) neuron which is a commonly used continuoustime SN [Fil and Chu, 2020, Gerstner and Kistler, 2002]. The LIF neuron is a powerful computational unit[G¨utig, 2014, G¨utig and Sompolinsky, 2006, Maass, 1997, Maass and Natschlaeger, 1998]. The CN inherits this.One of the attractive features of the (basic) CN neuron is that it is a micro-reversible model and as such itis a thermodynamically plausible model. While a thorough analysis of the energy requirements of the systemis beyond the scope of this article, we note that the physical plausibility of the model has highlighted resourcerequirements of the computation. In particular, we found that increasing the non-linearity comes at an additionalcost in resource. The CN suffers from starvation of B molecules as the chain-length of realising the activationfunction increases. For a sufficiently high number of m , this leads to a breakdown of the mechanisms and thesystem loses its ability to detect coincidences, as illustrated in fig. 6. The CN is better at solving the TC taskwhen the nonlinearity is higher. At the same time, the non-linearity requires resources that cannot be met anymore by the system. This “starvation” effect can be alleviated by increasing the bolus size (while keeping thethreshold fixed); see fig. 6. The minimal amount of work required to add particles to the CN is proportional tothe number of particles, i.e. the bolus size. Less visibly, there is also an additional expenditure of energy due tothe fact that the B and H molecules need to be removed from the system, which again requires a concentrationgradient to be maintained. Hence, at least in this model the non-linearity comes directly at a thermodynamicalcost. To the best of our knowledge, it is an open question whether the computation of non-linearities necessarilycomes at a higher thermodynamic cost or whether this is just a side effect of the computational medium usedhere, i.e. chemical reactions. References
A. Adamatzky, C. Fullarton, N. Phillips, B. De Lacy Costello, and T. C. Draper. Thermal switch of oscillationfrequency in belousov–zhabotinsky liquid marbles.
Royal Society Open Science , 6(4):190078, April 2019. doi:10.1098/rsos.190078. URL https://doi.org/10.1098/rsos.190078 .S. Afshar, R. Hamilton, J. Tapson, A. vanSchaik, and G. Cohen. Investigation of event-based surfaces for high-speed detection, unsupervised feature extraction, and object recognition.
Frontiers in Neuroscience , 12:1047,2019. ISSN 1662-453X. doi: 10.3389/fnins.2018.01047. URL .S. Afshar, N. Ralph, Y. Xu, J. Tapson, A. van Schaik, and G. Cohen. Event-based feature extraction usingadaptive selection thresholds.
Sensors , 20(6), 2020. ISSN 1424-8220. doi: 10.3390/s20061600. URL .U. Alon.
An Introduction to Systems Biology: Design Principles of Biological Circuits . Chapman & Hall/CRCComputational Biology Series. CRC Press LLC, 2019. ISBN 9781439837177. URL https://books.google.co.uk/books?id=MWXdQgAACAAJ .M. Amos.
Cellular Computing . Oxford University Press, 2004.P. Banda, C. Teuscher, and D. Stefanovic. Training an asymmetric signal perceptron through reinforcement inan artificial chemistry.
Journal of The Royal Society Interface , 11(93):20131100, April 2014. doi: 10.1098/rsif.2013.1100.D. Blount, P. Banda, C. Teuscher, and D. Stefanovic. Feedforward chemical neural network: An in silicochemical system that learns xor.
Artificial Life , 23(3):295–317, August 2017.M. Chen and J. Xu. Construction of a genetic conditional learning system in escherichia coli.
Science ChinaInformation Sciences , 58(11):1–6, September 2015. doi: 10.1007/s11432-015-5308-8. URL https://doi.org/10.1007/s11432-015-5308-8 .D. Chu. Limited by sensing - a minimal stochastic model of the lag-phase during diauxic growth.
Journal ofTheoretical Biology , 414:137–146, feb 2017. doi: 10.1016/j.jtbi.2016.10.019.D. Chu. Performance limits and trade-offs in entropy-driven biochemical computers.
Journal of TheoreticalBiology , 443:1–9, 2018. ISSN 0022-5193. doi: https://doi.org/10.1016/j.jtbi.2018.01.022. URL . Discussion D. Chu and D. Barnes. The lag-phase during diauxic growth is a trade-off between fast adaptation and highgrowth rate.
Scientific Reports , 6:25191, 2016. doi: 10.1038/srep25191. URL http://dx.doi.org/10.1038/srep25191 .D. Chu, N. Zabet, and B. Mitavskiy. Models of transcription factor binding: Sensitivity of activation func-tions to model assumptions.
Journal of Theoretical Biology , 257(3):419 – 429, 2009. ISSN 0022-5193. doi:https://doi.org/10.1016/j.jtbi.2008.11.026. URL .P. Diehl and M. Cook. Unsupervised learning of digit recognition using spike-timing-dependent plasticity.
Frontiers in Computational Neuroscience , 9:99, 2015. ISSN 1662-5188. doi: 10.3389/fncom.2015.00099. URL .C. Fernando, A. Liekens, L. Bingle, C. Beck, T. Lenser, D. Stekel, and J. Rowe. Molecular circuits for associativelearning in single-celled organisms.
Journal of the Royal Society Interface , (October 2008):463–469, 2009. doi:10.1098/rsif.2008.0344.J. Fil and D. Chu. Minimal spiking neuron for solving multilabel classification tasks.
Neural Computation , 32(7):1408–1429, 2020.W. Gerstner and W. M. Kistler.
Spiking Neuron Models: Single Neurons, Populations, Plasticity . CambridgeUniversity Press, 2002. doi: 10.1017/CBO9780511815706.C. Govern and P. ten Wolde. Energy dissipation and noise correlations in biochemical sensing.
Physical ReviewLetters , 113(25):258102, Dec 2014a. doi: 10.1103/PhysRevLett.113.258102. URL http://dx.doi.org/10.1103/PhysRevLett.113.258102 .C. Govern and P. ten Wolde. Optimal resource allocation in cellular sensing systems.
Proceedings of theNational Academy of Science USA , 111(49):17486–17491, Dec 2014b. doi: 10.1073/pnas.1411524111. URL http://dx.doi.org/10.1073/pnas.1411524111 .R. G¨utig. To spike, or when to spike?
Current Opinion in Neurobiology , 25:134–139, 2014. ISSN 09594388.doi: 10.1016/j.conb.2014.01.004.R. G¨utig. Spiking neurons can discover predictive features by aggregate-label learning.
Science , 351(6277):aab4113–aab4113, 2016. ISSN 0036-8075. doi: 10.1126/science.aab4113.R. G¨utig and H. Sompolinsky. The tempotron: a neuron that learns spike timing-based decisions.
NatureNeuroscience , 9(3):420–428, 2006. ISSN 1097-6256. doi: 10.1038/nn1643.A. Hjelmfelt, E. D. Weinberger, and J. Ross. Chemical implementation of neural networks and turing machines.
Proceedings of the National Academy of Sciences , 88(24):10983–10987, 1991. ISSN 0027-8424. doi: 10.1073/pnas.88.24.10983. URL .S. M. Hoffer, H. V. Westerhoff, K. J. Hellingwerf, P. W. Postma, and J. Tommassen. Autoamplification of atwo-component regulatory system results in ”learning” behavior.
Journal of bacteriology , 183(16):4914—4917,August 2001. ISSN 0021-9193. doi: 10.1128/jb.183.16.4914-4917.2001. URL https://europepmc.org/articles/PMC99548 .W. Maass. Networks of spiking neurons: The third generation of neural network models.
Neural Networks , 10(9):1659–1671, 1997. doi: 10.1016/s0893-6080(97)00011-7.W. Maass and T. Natschlaeger. Associative memory with networks of spiking neurons in temporal coding.
Progress in Neural Processing , 1998.J. Macia and R. Sole. How to make a synthetic multicellular computer.
PLOS ONE , 9(2):1–13, 02 2014. doi:10.1371/journal.pone.0081248. URL https://doi.org/10.1371/journal.pone.0081248 .J. Macia, B. Vidiella, and R. V. Sol´e. Synthetic associative learning in engineered multicellular consortia.
Journal of The Royal Society Interface , 14(129):20170158, April 2017a.J. Macia, B. Vidiella, and R. V. Sol´e. Synthetic associative learning in engineered multicellular consortia.
Journal of The Royal Society Interface , 14(129):20170158, 2017b. doi: 10.1098/rsif.2017.0158. URL https://royalsocietypublishing.org/doi/abs/10.1098/rsif.2017.0158 .S. McGregor, V. Vasas, P. Husbands, and C. Fernando. Evolution of associative learning in chemical networks.
PLOS Computational Biology , 8(11):1–19, 11 2012. doi: 10.1371/journal.pcbi.1002739. URL https://doi.org/10.1371/journal.pcbi.1002739 . CN chemical reaction network details D. N. Nesbeth, A. Zaikin, Y. Saka, M. C. Romano, C. V. Giuraniuc, O. Kanakov, and T. Laptyeva. Syntheticbiology routes to bio-artificial intelligence.
Essays in Biochemistry , 60(4):381–391, November 2016.E. Oja. Simplified neuron model as a principal component analyzer.
Journal of Mathematical Biology , 15(3):267–273, 1982. doi: 10.1007/BF00275687. URL https://doi.org/10.1007/BF00275687 .M. Okamoto, T. Sakai, and K. Hayashi. Biochemical switching device realizing mcculloch-pitts type equation.
Biol. Cybern. , 58(5):296–299, April 1988. ISSN 0340-1200.A. Racovita and A. Jaramillo. Reinforcement learning in synthetic gene circuits.
Biochemical Society Transac-tions , 48(4):1637–1643, August 2020.A. Sengupta, Y. Ye, R. Wang, C. Liu, and K. Roy. Going deeper in spiking neural networks: Vgg and residualarchitectures.
Frontiers in Neuroscience , 13:95, 2019. ISSN 1662-453X. doi: 10.3389/fnins.2019.00095. URL .T. Shirakawa and H. Sato. Construction of a molecular learning network.
Journal of Advanced ComputationalIntelligence and Intelligent Informatics , 17(6):913–918, 2013. doi: 10.20965/jaciii.2013.p0913.T. Yi, Y. Huang, M. Simon, and J. Doyle. Robust perfect adaptation in bacterial chemotaxis through integralfeedback control.
Proceedings of the National Academy of Sciences , 97(9):4649–4653, 2000. ISSN 0027-8424.doi: 10.1073/pnas.97.9.4649. URL . A CN chemical reaction network details
Input k IA = 10, k AI = 0 . k AB = 0 . k BA = 0 . k + = 1, k − = 5 k − last = 0 . k AE = 0 . k EA = 0 . k EH = 100, k HE = 0 . k AH = 0 . k HA = 0 . k HB = 100, k BH = 0 . k H ∅ = 0 . k B ∅ = 0 . Tab. 3: