Online Learning in a Chemical Perceptron

Chemistry provides many beneficial features that contribute to information processing, such as inherent parallelism, massive interactivity, redundancy, and asynchronicity [2, 12, 36]. Biomolecular systems have successfully tackled several computing problems, including the traveling salesman [1], 3-SAT [6], maximal clique [41], chess [17], and tic-tac-toe [49]. However, attempts to build a programmable molecular automaton, that is, an automaton with more than one (hard-wired) purpose, either failed or had limited scope and no reusability [4, 11, 42].

Our approach is to achieve programmability of a chemical system by learning and adaptation. Adaptation is one of the key aspects maintaining the functional, homeostatic closure of living systems [5], which are carried out by chemical processes. It enables individual organisms to adjust their decision-making schemes in constantly changing environments. Learning has been a vibrant topic in the ALife community for over two decades. It has been realized by means of neural networks [21, 44], various forms of evolutionary algorithms [38, 39], and reinforcement learning [50], in which agents learn from the consequences of their actions through rewards. The applications of learning include pathfinding problems [30], multi-agent systems [34], and robotics [8].

The idea of neural network computation in chemical systems is not new. Several theoretical or experimental DNA-based models [7, 9, 24, 25, 32, 36] have been proposed. Most recently, Qian et al. [52] demonstrated an experimental implementation of linear threshold circuits with DNA-strand-displacement seesaw gates and used these to construct a Hopfield network. The research in this area has mainly been limited to constructing logic gates and assembling them into circuits to compute custom Boolean functions. More important, the existing work does not dwell on the autonomous learning aspects of chemical neural networks. Typically [32, 52] the learning was performed by an external system that computes the weights for a formal neural network, before converting these to molecular concentrations to serve as parameters for the chemical implementation.

Spiking neural P systems [27, 28] are related types of systems that draw inspiration from neural network theory and incorporate membrane computing with a model of spiking neurons. Each neuron is wrapped in a membrane, where interneuron (intermembrane) communication is carried by the electrical impulses, called spikes. Some attempts were made to introduce learning into neural P systems [20]; there, however, similarly to DNA-strand implementations, learning has not been autonomous either. Also, P systems do not model reaction rates, and they have a different, more grammarlike update algorithm (kinetics) and manipulate discrete symbols.

Here, we model a two-input perceptron, a simple learning unit, as an artificial chemistry [13], where both linear integration and learning are implemented internally. To the best of our knowledge, it is the first such model. We restrict the interactions with the perceptron solely to injections of training instances consisting of two inputs and the desired output, and measurements of the concentrations of the output species. Learning requires a system to work continuously; therefore, a cleanup or a reset to a steady state is necessary. Again, we do not assume any manual reset. We show that the chemical perceptron can learn a logic function perfectly, and is also robust to perturbations of rate constants.

We simulate a general (artificial) chemistry based on Michaelis-Menten kinetics without any assumption about the molecular structure of the chemical species. Since the behavior of species is fully determined by reactions and rates, our model imposes basic constraints on how real molecular species should interact in order to adapt. This abstraction allows us to better understand the underlying principles, as well as the inherent challenges of adaptable artificial chemistries. After a conversion of Michaelis-Menten to mass-action kinetics, DNA strand displacement [46, 53] can potentially serve as a biochemical implementation of our chemical perceptron.

The contributions of this work are as follows.

• 1. 

  Our system is the first full-featured implementation of online learning in simulated artificial chemistry (Section 2), called the chemical perceptron (Section 4). Learning, as well as linear integration of weights and inputs, is handled internally.

• 2. 

  We present a systematic method that maps the variables of a formal two-input perceptron (Section 3) to species of the chemical perceptron (Sections 2.1, and 4.1).

• 3. 

  We implement two variants of a chemical perceptron, the weight-loop perceptron (Section 4.3) and the weight-race perceptron (Section 4.4), to demonstrate two qualitatively different approaches to linear integration, threshold comparison, and output production.

• 4. 

  The chemical perceptron is reusable, since it recovers its internal ready state after each processing (Sections 4.2, 4.3, and 4.4).

• 5. 

  The chemical perceptron learns perfectly all 14 linearly separable logic functions after 200 learning iterations (Section 7.1).

• 6. 

  The chemical perceptron is robust to perturbations of rate constants (Section 7.2). This property helps to substantially alleviate reaction-timing restrictions for real chemical implementations.

• 7. 

  We compare the weight-loop and weight-race perceptrons (Section 7.3) and discuss their potential biochemical implementation (Section 8).

Artificial chemistry (AC) is the standard framework for representing and simulating chemistry. An AC [13] consists of a set of molecular species or substances, S, and a set of reactions, R. Each reaction r ∈ R is an ordered pair X → Y, where both X and Y are multisets of species. Species from A are called reactants, and species from B are called products. Each reaction in our model is either a conversion (A → B or A + B → C + D), an annihilation (A + B → λ), or a decay (A → λ), where λ represents no species.

There are many types of ACs, which can be simulated with various techniques. In this article we employ macroscopic deterministic simulation, where species interact on the basis of their reactions with associated rate laws [15, 16, 23]. Each reaction has a rate, which defines the strength of the reaction's contribution to the production or consumption of particular species over time. An essential property of macro-chemistry is the absence of space. In a well-stirred tank, the probability that a molecule is involved in a reaction does not depend on its position, but on its type. Consequently, a multiset of species, where each species is characterized only by a concentration measured in moles per liter (M), describes the state of the system. For instance, a state of the AC with the species set S = {A0, A1, B} can be [A0] = 2 M, [A1] = 2 M, and [B] = 10 M.

By applying rate laws over all reactions, we obtain the change of concentration of molecular species as described by a system of ordinary differential equations (ODEs). Since it is, in general, impossible to find an analytical solution of such a system explicitly, we employ numerical integration of ODEs, which delivers an approximate solution [46]. Our AC numerically integrates concentration ODEs by the simple one-step Euler method [26, 48].

Deterministic simulation of rate laws is fast and provides a good approximation of real chemistry if the number of molecules in the reservoir is sufficiently high. Otherwise, a microchemistry based on stochastic collisions of individual molecules [18, 19, 29, 51] produces results with higher precision, but comes at a higher simulation cost.

An AC can represent a variable by one or several substances. In our chemical model we need to encode variables of two types: Boolean with values 0 and 1, and Real with values from ℝ. We transform variables to species in systematic fashion as follows.

A Boolean variable is represented by enumeration; one variable s requires two species, S⁰ and S¹, which are mutually exclusive. If the concentration of S⁰ is nonzero, then s = 0; analogously S¹ nonzero implies s = 1. A positive value of a Real variable directly corresponds to the concentration of a species. Since a concentration is never negative, a Real variable s needs to be represented by positive and negative species S^⊕ and S^⊖. If both variants S⁰ and S¹ for Boolean, or S^⊕ and S^⊖ for Real, are simultaneously present in the reservoir, they annihilate very rapidly.

Further, the zero concentration of species cannot represent a strict zero value. The problem is that intentional nothing cannot be distinguished from the state in which a system is still working and has not produced an output yet, or is still waiting for the input substances to enter the system.

The concept of AC actions or action series is an extension of the input configuration—species concentrations can be modified at times other than t₀. An AC action emulates a step in the execution of an experimental protocol, where at a certain time the person performing the chemical experiment mechanically injects or removes substances into or out from a tank. An action is modeled by instantaneously changing the concentration of a species. For iterative processes, such as learning, it is useful to define a repetitive AC action series, in which a sequence of actions repeat in a loop at predefined time intervals.

An AC translation, also performed from outside the system, is used to interpret concentrations of output species as the results of the chemical computation over a specified time interval. Translations operate on concentration ranges; hence they are defined by aggregate functions, such as max, min, and ∑. Similarly to AC action series, AC translation series can be repetitive.

Figure 1 presents a concentration plot of species S = {A0, A1, B} driven by reactions R0 and R1 as introduced in Section 2, using the rate constants k = 0.00325, k₀ = 0.025, and k₁ = 0.5. AC actions occur at times t₀, t₁₀₀, and t₂₀₀ as described in the caption. By applying the AC translation defined by max([A1]) > max([A0]) on the intervals t₀ – t₉₉, t₁₀₀ – t₁₉₉, and t₂₀₀ – t₂₉₉, we translate the output to the sequence 1, 0, 0.

Artificial neural networks [44] are inspired by the coarse-grained behavior of biological neurons in the brain. The perceptron is an early type of artificial neural network and one of the simplest systems capable of learning [45].

A perceptron is a single neuron that processes a vector of input signals x = (x₁,…; x_n), x_i ∈ ℝ, and produces one output y based on the setting of its weights w = (w₀, w₁,…, w_n) as shown in Figure 2. More precisely, a perceptron first calculates the linear integration (the dot product) of weights w and inputs x as , and then passes the result z to an activation function f : ℝ → [0, 1] or f : ℝ → [1, 1], which produces the final output y. Note that the weight w₀, called the bias or offset, always contributes to an output, since its associated input x₀ is a constant 1.

A perceptron can classify only linearly separable functions [37]—functions in which a straight line or, in the general case, a hyperplane can divide the inputs into two classes. By combining several perceptrons, we can construct a multilayer perceptron network, also known as a multilayer feed-forward network [21], that overcomes the linear separability problem and in fact becomes a universal approximator.

Perceptron learning [44] is a type of supervised Hebbian learning [22] where a training data set T = {(x₁, d₁),…, (x_m, d_m)} consisting of input-output pairs characterizes the target behavior of the system. During each step of the learning process, a perceptron absorbs one training sample (x, d). If there is a discrepancy between the actual output y and the desired output d, the error is fed back to the perceptron and triggers an adaptation of the participating weights. The adaptation of the weight w_i for the training sample (x, d) at time t is defined as w_i(t + 1) = w_i(t) + α(d − y(t))x_i, where the learning rate α ∈ (0, 1] represents the adaptation strength.

If an error is detected, that is, if |d − y| > 0, the weight w_i shifts toward the desired output if its input signal x_i is nonzero. Conversely, if an input x_i = 0, the weight w_i is not involved in the global output y and therefore stays unaltered. Initially, the weights are set to small random values. The process of weight adaptation continues until the cumulative error of several consecutive training samples drops below the error threshold, or alternatively a fixed number of iterations is reached.

In this article we model a specific type of perceptron: the two-input perceptron with binary inputs x₁ and x₂, and real-value weights w₀, w₁, and w₂. The activation function f = sgn (hard delimiter) outputs one if the dot product z = w₀ + w₁x₁ + x₂w₂ is positive, and zero otherwise. Therefore, the linear integration part is reduced to four cases as presented in Table 1a.

Now, we investigate whether all weights must support both positive and negative values. Assume w₀ is negative and w₁ and w₂ are positive. Then, for instance, the weights w₀ = −10, w₁ = 7, w₂ = 9 model the AND function, and the weights w₀ = −10, w₁ = 12, w₂ = 13 the OR function. However, no combination of negative w₀ and positive w₁ and w₂ weight values can represent (e.g.) the NAND function. Table 1b summarizes the limitations of all sign-weight combinations for a representation of logic functions. It shows that each of the weights w₀, w₁, and w₂ must support both positive and negative values to implement a perceptron that can encompass all 14 linearly separable binary functions, that is, all binary functions except XOR and NXOR.

In this section we describe the implementation of the two-input perceptron by means of an artificial chemistry. We want to emphasize that there are many ways to approach this problem. Here we present two models—the weight-loop perceptron and the weight-race perceptron—which represent two fundamental, substantially different techniques for a calculation of the weight sum and the zero threshold. Before we describe the weight-loop and weight-race perceptrons in detail, first we need to formalize the features shared by both models.

The representation of the formal perceptron's variables by chemical species follows the scheme presented in Section 2.1. Each variable of type Boolean needs 0-value and 1-value species, and each variable of type Real splits into ⊕ and ⊖ species variants. Table 2a presents the core species of a chemical perceptron, along with their mappings to the variables of the formal perceptron described in Section 3.

The inputs x₁ and x₂ of type Boolean trigger different processing paths in the perceptron, based on their values. Therefore, each input-variable value pair requires its own species: X₁⁰ and X₁¹ to represent x₁, and X₂⁰ and X₂¹ to represent x₂. Note that if a single species encoded an input variable, we would need to differentiate values 0 and 1 based on low versus high concentration, which would lead to a more complicated design. Similarly, variables y and d are also Booleans, so the same representation applies: Y⁰, Y¹ and D⁰, D¹. Weights of the two-input binary perceptron need to support both positive and negative values (Table 1b), so each weight w_i splits into two distinct positive and negative variants, species W_i^⊕ and W_i^⊖, where i ∈ {0, 1, 2}. For simplification, we will use the name of a group or a subgroup to refer to all associated species as defined in Table 2a. For example, the species subgroup X₁ includes both X₁⁰ and X₁¹, and the group X includes species X₁⁰, X₁¹, X₂⁰, and X₂¹.

The chemical perceptron can function in two modes: binary-function mode and learning mode. In the binary-function mode, the perceptron basically acts like a logic gate; it takes two inputs X₁ and X₂, and produces an output Y. The second learning mode is built on top of the output production, so again inputs must be present. The learning is triggered by the desired-output molecules D. Once output Y is produced, it is compared against D, and if they differ (i.e., Y⁰ is produced but D¹ is expected, or vice versa), positive- or negative-weight molecules are created and added to existing ones, to modify the weights.

After each learning iteration, the perceptron needs to recover to its steady state. Hence, transient species not consumed during the chemical computation must be removed by some cleanup reaction such as decay, because there is no external cleanup. Only the weight species form the persistent state of the chemical perceptron; hence, all other species are considered transient. Since decay and annihilation reactions would break the mass-conservation principle, we assume that instead of nothing (λ), an inert by-product is produced.

Due to the multiplicity of input, output, desired-output, and weight species, reactions of the same type are collected into groups, which simplifies the reasoning as well as the simulation of the chemical perceptron. Reactions belonging to the same group share common structural characteristics, catalysts, and inhibitors, as well as rate constants.

In the following sections we present the species and reactions of the weight-loop and weight-race perceptrons. The actual setting of the reaction rate constants, which is optimized by the genetic algorithm, is discussed separately in Section 6.

The weight-loop perceptron follows the formal perceptron definition from Section 3 quite rigorously. It consists of 21 species and 34 reactions. Apart from the core species defined in Section 4.1, the weight-loop perceptron requires the processed-weight species and , and fuel species E (Table 2b).

The weight-loop perceptron computes the weight sum directly by transforming weights W into output species Y. The problem is that the weights encode the state of the perceptron, so their concentration must be preserved. Therefore, apart from Y species, the perceptron must also create backup copies of the weights, . The perceptron can then restore its weights after the output production is over. A reaction followed by would break the mass-conservation law, so the perceptron needs to consume a fuel, species E, which is provided to the system at constant concentration 1 M. From a functional perspective, the perceptron sequentially processes an input, produces an output, recovers weights, and finally performs a cleanup (Figure 3).

The perceptron starts working when inputs X₁ and X₂ are injected into the system. The perceptron processes the weight W₁ on input X₁¹, the weight W₂ on input X₂¹, and the weight W₀, in parallel, producing Y⁰ and Y¹ molecules. Species X₁¹, formally encoding x₁ = 1, catalyzes ⊕ and ⊖ versions of reaction . Similarly, species X₂¹, which represents x₂ = 1, catalyzes . Since the weight W₀ always contributes to the sum regardless of input, each of the possible inputs X₁⁰, X₁¹, X₂⁰, and X₂¹ catalyzes . In order to determine whether the total concentration of Ys is above or below the zero threshold, we let Y⁰ annihilate with Y¹. If there are more Y⁰ molecules at the end, the output is 0; otherwise 1. Weights could alternate between the normal version W and the processed version , each time consuming a fuel E and producing new Y molecules. To prevent a continuous cycling of the weights, the weight-loop perceptron must ensure that there is no input present before it rolls the weights back. That is, input species must decay, and the processed weights roll back only when substantial amounts of inputs are gone—that is, inputs act as inhibitors on reactions. The output molecules Y are removed from the system by a decay. Table 3a presents the full set of reactions with associated catalysts and inhibitors.

Only eight reactions are needed for the learning part (Table 3a, groups 9–11). During each learning step, the actual output Y is compared against the desired output D. If Y matches D (i.e., the substances Y⁰ and D⁰, or Y¹ and D¹, are simultaneously present in the system), the output is correct. In this case, no learning is needed and the weights remain unaltered. Otherwise, the desired-output species D transforms to ⊕ or ⊖ versions of the weight species W, but only for those participating in an output production for current inputs. Thus, an input and an output together catalyze the D → W reactions, so they are dependent (and) catalysts (Table 3a, group 11).

For instance, if the perceptron produces the output Y⁰ for the input species X₁⁰ and X₂¹, but the desired output D¹ is injected, then reactions D¹ → W₂^⊕ and D¹ → W₀^⊕ are triggered, and weights W₂^⊕ and W₀^⊕ are produced and added to (or annihilate with) existing weights. The strength of the adaptation (i.e., the learning rate α) is incorporated into the concentration of the desired output species D. For example, if 10 M of D¹ is injected into the system in the previous example, that 10 M is distributed between production of W₀^⊕, and W₂^⊕ (a small amount of D¹ actually disappears because of a decay reaction).

Since the system is open and weights W can switch reversibly to , consuming a fuel E provided from outside, an infinite loop might emerge, in which the concentrations of Y molecules increase without bound. The correct timing of phases is crucial for avoiding this problem.

The functioning of the weight-loop perceptron is based on rather conservatively designed phases working in a sequence. This approach works well, since there is almost a one-to-one relation between the routines of the formal perceptron and those of the chemical perceptron. Nevertheless, the idea of direct calculation of the weight sum and recovering the original state seems unnecessarily cumbersome for a chemical system.

The weight-race perceptron improves on the weight-loop model by switching the chemical roles of inputs and weights—that is, instead of having inputs catalyzing a transformation of weights to a weight sum, which determines an output, weights simply catalyze the input-to-output reactions as presented in Figure 4. Thus, the perceptron does not compare weights directly, but lets them compete on input-to-output reactions as catalysts, so it basically implements a rate (derivation-based) comparison. For this to work, species Y⁰ must annihilate with Y¹ quickly, racing must be simultaneous, and the rate functions must have similar shapes.

The weight-race chemical perceptron consists of 14 species and 30 reactions. Input species X₁ and X₂ are transformed directly to Y⁰ or Y¹, depending on the signs of the currently present weights. Weight W₁ solely catalyzes reaction X₁¹ → Y¹ for variant W₁^⊕, and X₁¹ → Y⁰ for variant W₁^⊖. Similarly, weight W₂ drives two reactions, in which X₂¹ is the reactant. Since weight W₀ is always active, it drives all possible input-to-output reactions. Weights catalyze the reactions concurrently, so the one with the highest concentration consumes the largest portion of an input and therefore has the highest contribution in an output. Analogously to the weight-loop model, an annihilation of Y¹ and Y⁰ decides whether the concentrations of the Y's are above or below the zero threshold. The full set of reactions with associated catalysts is presented in Table 3b. The learning part is the same as in the weight-loop model.

If the weight-race perceptron is to treat all weights equally, it must ensure that the weight race is fair. Following the formal perceptron definition, the contribution of weights in the sum must be uniform, meaning there is no preference among weights. Apart from the concentrations of weights, the reaction rate constants determine the actual speed of the input consumption. Now, if all weights have the same sign, then it does not matter what rate constants are set, so let the qualitative state of the perceptron be W₀^⊕, W₁^⊖, and W₂^⊖ (Figure 5).

For inputs X₁⁰ and X₂⁰, only the weight W₀^⊕ is active, so there is no racing. The weight W₀^⊕ competes with W₁^⊖ and W₂^⊖ for inputs X₁¹ and X₂¹; however, W₀^⊕ is privileged, since it consumes X₁¹ and X₂¹ at the same time. Note that W₁^⊖ consumes just X₁¹, and W₂^⊖ just X₂¹. To avoid this problem we have to penalize W₀^⊕ by setting the rate constants of the reactions X₁¹ → Y¹ and X₂¹ → Y¹, both catalyzed by W₀^⊕, to δ, and those catalyzed by W₁^⊖ and W₂^⊖ to 2δ. As a result, the contribution preference of the weights is balanced.

A new problem emerges for inputs X₁¹ and X₂⁰, where W₀^⊕ drives two, but W₁^⊖ just one reaction. Even if the X₁¹ reactions follow the two-to-one rate constant ratio, whatever constant is assigned to the reaction X₂⁰ → Y¹ catalyzed by W₀^⊕ will result in an unfair advantage for W₀^⊕, since eventually all X₂⁰ molecules will change to Y¹. To balance the preference of W₀^⊕ we need to introduce a decay of X₂⁰, such that exactly two-thirds is taken away. In order to do that, W₀^⊕ must catalyze not just X₂⁰ → Y¹ with the rate constant δ, but also the decay X₂⁰ → λ with the rate constant 2δ. That is the reason why the reaction set of the weight-race perceptron contains a decay of the input species X₁⁰ and X₂⁰ (Table 3b, group 2).

The two-to-one ratio of rate constants must hold if the goal is to model the perceptron with no preference among weights. However, as we show in Section 6, if the weight-race perceptron does not obey this ratio and has a bias on a particular weight, it can still perform well.

Compared to the previous model, the weight-race perceptron is substantially simpler. It is minimal in number of species (only the core set is used), and it contains just 30 reactions, without any inhibition. Unlike the weight-loop perceptron, the system does not need any externally supplied fuel species. In fact, the input species adopt this role, so they are essentially an information and energy source.

We have introduced the chemical perceptron models structurally as collections of species and reactions with catalysts and inhibitors. Now, we shall specify the setting of executions (or chemical experiments) in terms of the action series and the translation series (Section 2.2).

Since the input processing and the output production take some time, we cannot inject input species X₁ and X₂ immediately after the previous pair, but we have to wait for a certain number of simulation steps, S = 5,000. In the binary-function mode, the first action at time t₀, handling an initialization sets the weights W according to the target logic function. Then every S steps we execute one of the four actions (one per input pattern) for the weight-loop perceptron:

• • 

  [X₁⁰] = 1 M, [X₂⁰] = 1 M,

• • 

  [X₁¹] = 1 M, [X₂⁰] = 1 M,

• • 

  [X₁⁰] = 1 M, [X₂¹] = 1 M,

• • 

  [X₁¹] = 1 M, [X₂¹] = 1 M.

In the weight-loop perceptron, inputs serve as catalysts of the weight-to-output reactions, which consume fuel species E. Since the weight-race perceptron directly transforms input species to output species, the concentration of inputs must be higher:

• • 

  [X₁⁰] = 2 M, [X₂⁰] = 2 M,

• • 

  [X₁¹] = 2 M, [X₂⁰] = 2 M,

• • 

  [X₁⁰] = 2 M, [X₂¹] = 2 M,

• • 

  [X₁¹] = 2 M, [X₂¹] = 2 M.

Figure 6 presents simulations of the weight-loop perceptron and the weight-race perceptron, each computing the NAND function on four consecutive inputs.

In learning mode, the initial concentrations of weights are generated randomly in the interval 2–10 M, with equal probability of selecting positive and negative variants. A learning rate α, which is constant throughout the whole training, is incorporated into the concentration of the desired output D.

The original definition of perceptron learning (Section 3 adjusts weight w_i by Δw_i = α(d − y)x_i for a given output y and desired output d. Assuming d ≠ y, each weight participating in the output production increments by Δw = α(d − y). Since the sign of the weight sum fully determines output, weight adaptation is stronger for inputs with the higher number of ones. For instance, the weight sum is adjusted by Δw for input (0, 0), but by 3Δw for input (1, 1), as shown in Table 4a.

In the chemical perceptron, the concentration of the desired output is divided between weights; hence the weight adaptation of the original perceptron would correspond to the concentration profile of the desired output for four consecutive inputs: [D] = α, [D] = 2α, [D] = 2α, [D] = 3α. Our simulations proved that this unfairness results in a split of performance for function pairs, such as CONST0, CONST1 or AND, NAND. Since these functions are the inverse of each other, the chemical perceptron should learn them with the same success.

Essentially, the uniform adaptation of individual weights causes a bias in the weight adaptation. To avoid that, we do not adapt individual weights, but the whole weight sum uniformly for all inputs, as presented in Table 4b. More specifically, the chemical perceptron divides Δw, the concentration of the desired output D among weights, so the whole weight sum is adjusted by Δw. Note that a small amount of the desired output disappears because of a decay.

The concentration of the desired-output species D is constant for all inputs in the chemical perceptron. By experiments we determined that the optimal concentration of D is 2 M. If the concentration was too high, the weights would oscillate and would not converge on a stable solution. Conversely, a low concentration of D prolongs the learning process and does not provide enough pressure to drive weights out of the zero region if their concentrations are very low.

As opposed to our chemical perceptron, the biased adaptation in the original formal perceptron does not cause substantial problems, because the weight sum is further processed by an activation function and the learning rate α decreases over time. As a result, small differences in the weight adaptation become unimportant.

Let f : {0, 1} × {0, 1} → {0, 1} be a target two-input logic function, which we wish to teach to the perceptron. By D^{f(x₁, x₂)} we denote species D⁰ if f(x₁, x₂) = 0 or D¹ if f(x₁, x₂) = 1. The training set of the weight-loop perceptron consists of four actions:

• • 

  [X₁⁰] = 1 M, [X₂⁰] = 1 M, [D^f(0,0)] = 2 M,

• • 

  [X₁¹] = 1 M, [X₂⁰] = 1 M, [D^f(1,0)] = 2 M,

• • 

  [X₁⁰] = 1 M, [X₂¹] = 1 M, [D^f(0,1)] = 2 M,

• • 

  [X₁¹] = 1 M, [X₂¹] = 1 M, [D^f(1,1)] = 2 M.

Similarly, the training set of the weight-race perceptron is:

• • 

  [X₁⁰] = 2 M, [X₂⁰] = 2 M, [D^f(0,0)] = 2 M,

• • 

  [X₁¹] = 2 M, [X₂⁰] = 2 M, [D^f(1,0)] = 2 M,

• • 

  [X₁⁰] = 2 M, [X₂¹] = 2 M, [D^f(0,1)] = 2 M,

• • 

  [X₁¹] = 2 M, [X₂¹] = 2 M, [D^f(1,1)] = 2 M.

Learning consists of a series of actions, each randomly chosen from the training set and performed every S steps. The total number of actions, L, per action series is L_f = 120 for the fitness evaluation (Section 6), and L_p = 200 for the learning performance and robustness analysis (Sections 7.1 and 7.2), so the perceptron runs for either S × L_f = 6 × 10⁵ or S × L_p = 10⁶ time steps.

If we injected inputs together with the desired output at the same time, the adaptation of the weights would start immediately, changing the actual output, so the actual output would differ from the one we would obtain by providing only input species. In the extreme case, the chemical perceptron could just copy the desired output to the actual output by having very low concentrations of weights. To prevent this, we inject an input, wait a certain number of steps, measure the output, and then provide the desired output. Note that a reaction D → W requires both catalysts, the input species X and the output species Y, to have a sufficient concentration at the moment of weight adaptation. Therefore, we must allow enough time for the output production, but we cannot postpone the injection of D for too long; if we did, the chemical perceptron would process or decay all input species. We found experimentally that this delay can be fairly short. More precisely, in our learning simulations we inject the desired output 100 steps after the input species.

Now, the only question is how to interpret the output, or in other words, how the concentrations of species Y⁰ and Y¹ translate to the value of the variable y from the formal definition. Since y is a Boolean variable, the translation compares the concentrations of value-0 species and value-1 species just before a desired output is injected. Hence, the value of the variable y is defined as [Y¹] > [Y⁰] at relative time step 99.

Since the number of rate constants of the chemical perceptrons is very large, it would be difficult and time-consuming to set them manually in a trial-and-error fashion or by exhaustive search. We therefore optimize the setting of the rate constants by employing a standard genetic algorithm (GA) [14, 43].

Because the reactions from the same group are structurally similar and in fact they share the same rate constants, the total number of representative rate constants constituting a possible solution, known as a chromosome, can be reduced substantially: from 68 to 21 for the weight-loop and from 52 to 14 for the weight-race perceptron. Figure 7 shows the chromosome structure for both perceptrons.

The fitness of a chromosome is evaluated as the performance of a chemical perceptron with the given rate constants (from the chromosome) on the binary-function learning task. Whether a chemical perceptron has learned a given function depends primarily on its final state; hence, the fitness embraces the performance for only the last steps. More precisely, the fitness is the Hamming distance between the last P = 20 translated actual outputs y_L−P,…, y_L and the desired outputs d_L−P,…, d_L. Because each learning action series starts with a random setting of weight concentrations, we calculate the fitness over 30 runs for each of the 14 functions.

Due to the nondeterministic nature of the action series, we find the fitness score more comparable across the current population if the chromosomes obtain the same problem instances by preevaluating action series. The GA combines elite selection with one point crossover and per-element (rate constant) mutation. Since only certain values of rate constants are plausible, we restrict their value ranges for the mutation and the generation of the initial population as shown in Table 5. The setting and constants of the GA are: population size M = 100, elite size E = 20, crossover probability p_c = 0.9, per-element mutation probability p_m = 0.2, uniform mutation strength s_m = 0.5, and generation limit G = 40.

The GA reaches solutions with fitness above 0.9 within a couple of generations, and then it continues at a slower pace toward the maximum fitness of 1, which is reached around generation 20 for both perceptrons (Figure 8). Since the learning action series are random and the simulation cost of a more precise fitness evaluation of each chromosome is too high, small fluctuations of the fitness occur even after the maximum has been reached. This is, however, not critical, since the best chromosomes found by the GA (Table 6) can learn all linearly separable logic functions with a 100% correct rate, as demonstrated in Section 7.1.

Overall, the fitness landscape of the rate constants for both perceptrons has the shape of a high plateau; hence, finding acceptable rate constants is not difficult. This demonstrates that our structural design of the chemical perceptrons in terms of species and reactions already provides correct behavior, and the perceptrons do not need to rely on the specific rate constants. In fact, we show in Section 7.2 that both models are extremely robust to the perturbation of rate constants.

We want to stress that in this article we consider the GA as just a tool, rather than an objective of our research. Since our GA setting produced satisfactory results, we have not explored optimization of the evolutionary process any further.

In Section 4.4 we demonstrated that the weight-race perceptron must follow the two-to-one rate constant ratio of X → Y reactions for the weights W₁ and W₂ versus the weight W₀, if the goal is to model the formal perceptron, where the weights contribute equally to the weight sum. Nonetheless, rate constants that do not obey this ratio can still perform well. Indeed, the best rate constants (chromosome) obtained by the GA (Table 6, right) makes the weight W₀ catalyze input-to-output reactions more strongly than the rest of the weights. Note that all reactions for weight species W₁ and W₂ are symmetric. The weight-race perceptron balances this preference by a lower consumption of desired output molecules D for an adaptation of weight W₀. Still, because of the linearity and monotonicity of adaptation, the feedback loop adjusts the weights correctly regardless of the weight preference.

We evaluated the performance of both the weight-loop perceptron and the weight-race perceptron with the best rate constants found by the GA. Similarly to the fitness evaluation, the perceptrons were run against 14 action series; however, in order to achieve higher precision, this process was repeated 10⁴ times as opposed to 30 times for the fitness evaluation, and each action series consisted of L_p = 200 as opposed to L_f = 120 actions (training iterations).

Figure 9 shows the weight-loop and the weight-race perceptron learning the NAND function. The training starts by the setting of weights as the constant 0 function. After 16 training iterations (8 × 10⁴ time steps), both perceptrons successfully adapt to the NAND's input-output characteristics. Furthermore, they smoothly handle the transition from W₀^⊖ to W₀^⊕ species and low concentration of the output species. An interesting feature is that a chemical perceptron continues to improve its weights even after the qualitative 1, 1, 1, 0 solution is found, trying to further strengthen the output signal.

Figure 10 presents the performance for all linearly separable binary functions averaged over 14 × 10⁴ runs. The results show that both chemical perceptrons can successfully learn 14 logic functions and reach the perfect score of 100% in all cases. The average error of each perceptron is around 6 × 10⁻³% after 120 training iterations. The error occurs when the concentrations of weights with opposite signs, such as W₀^⊕ and W₁^⊖, are (almost) identical. Then the perceptron does not produce enough output to catalyze weight adaptations and learning slows down. However, this situation is rather sporadic, and after running the perceptron for 80 more learning iterations, the error disappeared.

The chemical perceptrons perform almost identically. Moreover, their behaviors are also very similar under perturbation of the rate constants. The perturbation strength p defines how much the constants can change. Given p, each rate constant is perturbed uniformly over the interval (0, p). Thus, for a perturbation q ∈ (0, p) a rate constant γ changes to (1 ± q)γ, where increment and decrement have equal probability. The perturbation is not limited by the bounds used for the GA (Table 5); nonetheless, when p > 1, a perturbation to a negative value is replaced by zero, since the rate constants cannot be negative.

We analyzed the robustness of perceptrons with the best rate constants only. As presented in Figure 11, the chemical perceptrons maintain high robustness; even for p = 0.5, the mean correct rate after 200 learning iterations is 98.98% for the weight-loop perceptron and 99.34% for the weight-race perceptron. The main difference between the perceptrons in robustness occurs with high perturbation strengths, at which the weight-loop perceptron becomes slightly more vulnerable. For p = 2.0 the performance of the weight-loop perceptron even drops below 50%, which is a sign that the concentration of output species rises beyond the maximal limit in some simulations. Recall that for the weight-loop perceptron this situation might happen due to an open fuel influx.

Table 7 summarizes the attributes of our perceptron models. They perform almost identically and have very similar robustness characteristics; nonetheless, the weight-race perceptron is simpler, because it consists of fewer species and reactions and does not require any inhibition. On the other hand, the weight-loop perceptron is a closer match to the formal perceptron.

One potential framework for a biochemical implementation of our chemical models, especially the weight-race perceptron, is DNA strand displacement [53]. The use of DNA allows one to pick arbitrary sequences to stand for an arbitrary species from an artificial chemistry. In strand displacement systems, populations of these species are typically represented by the populations of single-stranded DNA molecules. These interact with double-stranded gate complexes, which mediate transformations between free signals. Soloveichik et al. [46] proved that a strand displacement circuit can approximate, with arbitrarily small error, any artificial chemistry based solely on mass-action kinetics. This offers a way to derive a DNA implementation of our perceptron models. First, we would need to translate the whole model into mass-action kinetics. The catalytic reactions that we previously modeled using Michaelis-Menten kinetics could be handled by expanding out the full set of enzymatic reactions. On the other side, the uncompetitive inhibitory reactions cannot be directly represented by mass action and therefore must first be adapted to the competitive version. Because of this complication, the weight-race perceptron, which does not use an inhibition, is more suitable for mass-action transformation and thus also for DNA strand displacement implementation.

An expanded mass-action version of the weight-race perceptron would undoubtedly be larger than the model presented in this article. We estimate it would consist of 40–50 species and 50–60 reactions, which is high but still manageable. Moreover, we might be able to simplify the encoding, for example by modeling the ⊕ and ⊖ and the 0 and 1 variants of species as two complementary strands, where annihilation is actually the production of an inert double-stranded complex. The state of the art in strand displacement circuits is a four-bit square root circuit [52] including 130 strands (74 non-input species). The number of species is higher than we need for our perceptron, but the square root circuit was built with seesaw gates, which are much simpler than the gates in Soloveichik's encoding. While the expanded model of a chemical perceptron would still be impractical for implementation in the laboratory in the near future, it would at least demonstrate that our model could in principle be implemented using real DNA strands. Moreover, the high robustness of the perceptron models to rate constant perturbations means that precisely reproducing the optimal rate constants might not be necessary in a strand displacement implementation.

There are several potential areas for improvement. First, instead of the fast but rather coarse-grained Euler method for the numerical integration of ODEs, the chemical model simulation could incorporate the more accurate Runge-Kutta approximation [48], or possibly some stochastic method [29, 51]. Similarly, lowering the time step would deliver higher-precision results, but at larger simulation costs.

All reactions in our model are irreversible, which simplifies reasoning about their causality. However, a reversible reaction with rates based on equilibrium of the forward and the backward part is more plausible and could cover an irreversible reaction as a special case. Also, the structural aspects of molecular species as an intermediate level between artificial and real chemistry would be worth investigation.

Since the weight-race perceptron excludes inhibition and in it all catalytic reactions with the exception of weight adaptation require only one catalyst, it corresponds more closely to real chemistry. We argue that after small adjustments, the reactions with two dependent catalysts, in which an input species and an output species together catalyze the weight adaptation, can broaden to the standard one-catalyst Michaelis-Menten kinetics. Then, the model can be further adapted, so it follows pure mass-action kinetics, where a catalysis is expanded to partial association and disassociation reactions.

Our AC model employs three concentration thresholds. The lower threshold avoids numerical instability by cutting the concentration to zero if it goes below 10⁻³ M. The upper threshold of 50 M detects whether the concentration of a species diverges beyond bounds (explodes). The difference threshold, 10⁻⁵ M, defines the minimal measurable change in the concentration of a species. If the concentration of a species increases or decreases very slowly, it is considered as a constant. Then if the concentration of every species is constant, the system is at the fixed point, and hence the AC run can proceed directly to the next input (action), and the AC simulation cost can be reduced.

To make the system more biologically realistic, we can wrap a perceptron in a permeable compartment, which provides integrity, controls input-output streams, and avoids duplication of species in case the chemical system consists of multiple perceptrons. Then an efflux of obsolete, noninteracting product might serve as a decay reaction.

The chemical perceptrons do not age; hence, the learning rate (or the weight-adaptation strength) remains constant through the whole learning process. A decrease of the learning rate over time usually improves convergence. In a chemical system we cannot assume that the perceptron knows when training starts and ends; therefore, if we assumed that the perceptron ages, this process would need an inner quality of the perceptron that degrades (decays) over time.

Our perceptron setup might present a modularity problem in the case of connecting multiple perceptrons together—the output from a perceptron in the first layer would not necessarily be at an appropriate level to feed into a perceptron in the second layer, since the input concentrations are fixed but the output concentrations are not. As a consequence, the execution setup with fixed interval length between actions and with constant concentration of input and desired-output species could have variations in order to obtain perceptrons with more generality, or a smaller bias to the specific concentration range. Also, it would be interesting to address the learning capabilities of the system in the face of noise in the training instances provided. We could introduce noise into reaction rates similarly to that in our robustness analysis, but at each step during the simulation rather than just initially. Finally, it would be worthwhile to explore repair capabilities of a chemical perceptron from the perspective of systems biology or relational biology. Decay of currently stable weight species would compel a chemical perceptron to repair its inner state perpetually by consuming an external resource.

We have demonstrated a proof of concept of learning and adaptation in an artificial chemical system. Artificial chemistry provided a conceptual framework, using sets of species, reactions, and reaction rates for modeling the functional characteristics of a two-input binary perceptron. After each processing, the perceptron recovers its default ready state, so it is reusable without outside intervention. We have introduced two models—the weight-loop perceptron and the weight-race perceptron. The weight-loop perceptron calculates the weight sum directly by transforming the weights to output and then recovering them to the original concentration. In the weight-race perceptron, the weights catalyze the input-to-output reactions; hence, the weight sum (the output) is calculated indirectly by weight competition (rate comparison).

Both chemical perceptrons can be trained perfectly to a desired logic function by providing sufficient input-output training pairs. Furthermore, they maintain high robustness to the perturbation of rate constants and achieve a 99% success rate for a 50% perturbation. Therefore, an implementation of our perceptrons in real chemistry would not have to use our specific rate constants in order to perform well. Since the weight-race perceptron's reactions are simpler and directly transformable to pure mass-action model, it is more practical and suitable for real chemical experimentation.

By implementing the mass-action model as a DNA-strand-displacement circuit, the chemical perceptron can serve as the basis of a new abstract layer for biochemical computing. Such a programming interface would hide inner chemical “hardware,” so we could program it to the desired function without knowing anything about chemistry.

This research was funded by NSF grant 1028120.