Exploring the programmability of autocatalytic chemical reaction networks

Influencing the rate of autocatalysis

We first examined the influence of the metal ions (Ca²⁺, La³⁺ and Nd³⁺) on the rate of trypsin autocatalysis under batch conditions. Tg (100 μM), Tr (1 μM) and Ca²⁺ (20 mM) were initially mixed, and we monitored the change in the concentration of Tr, [Tr], using a standard assay with Na-Benzoyl-DL-arginine 4-nitroanilide, Nα-Benzoyl-DL-arginine 4-nitroanilide hydrochloride (BAPNA). Figure 2a depicts a typical sigmoidal curve—a sharp transition that is characteristic to an autocatalytic conversion—can be initiated earlier when La³⁺ or Nd³⁺ was added to the mixture. In greater detail, Fig. 2b–d depicts the [Tr] at three specific time points of a set of experiments wherein we varied the concentration of the three metal ions (Supplementary Fig. 3). Two important effects were observed in the examined range of concentrations of the ions: (I) La³⁺ can accelerate the rate of autocatalysis, and the acceleration is much stronger than when Ca²⁺ was used. Figure 2b, c show that [Tr] increases until it saturates as the initial concentration of Ca²⁺ and La³⁺ ([Ca²⁺]₀ and [La³⁺]₀) are increased. (II) Nd³⁺, in contrast, can accelerate and decelerate the rate of autocatalysis. Figure 2d shows that the [Tr] increases when the initial concentration of Nd³⁺ ([Nd³⁺]₀) increases from 0 to 0.3 mM. Upon further increase of [Nd³⁺]₀, however, the [Tr] decreases. This pattern was established in any of the three time points. Overall, the observations suggest that all examined ions change the apparent rate of autocatalysis, affecting the rate to a different extent.

a Autocatalytic conversion of Tg into Tr in the presence of Ca²⁺. The addition of La³⁺ and Nd³⁺ significantly increase the rate of Tr autocatalysis. b–d Concentration of trypsin, [Tr], as a function of the metal ion and its concentration at fixed time after start of the reaction. Depicted are data at three time points (for extended data see Supplementary Fig. 3). SD bars are based on independent measurements (n = 2). Initial conditions: [Tg]₀ = 100 µM and [Tr]₀ = 1 µM. Initial conditions for a–d: [TRIS-HCl] = 100 mM (pH 7.8), T = 21–23 ^oC.

How the metal ions (X) affect Tr autocatalysis is proposed in Fig. 3. Based on the batch experiments, we assume that (I) autocatalysis cannot occur without any of the metal ions; (II) binding of the ions to Tr (depicted as [TrX]) can ‘activate’ autocatalysis, changing the apparent rate constant for the conversion of Tg into Tr (indicated by k_x); (III) binding of Nd³⁺ to Tg (depicted as [Nd³⁺Tg]) can ‘de-activate’ autocatalysis, slowing down the conversion of Tg into Tr. Notably, this effect was not observed when La³⁺ or Ca²⁺ were used. We performed control experiments using gel electrophoresis (i.e., SDS-PAGE) for samples from the reaction mixtures to examine the product distribution as a function of [Nd³⁺]₀ (Supplementary Fig. 4). The gels revealed that, only at optimal [Nd³⁺]₀ (0.3–0.4 mM), Tg is fully hydrolysed after 16 min of the reaction.

Both Tr and Tg form 1:1 complexes with metal ions but only Tg has a second binding step for Nd³⁺. For a detailed description of the mechanism, see Supplementary Methods Mathematical modelling.

Intrigued by this observation, we investigated if the effect of [Nd]₀ on autocatalysis can be explained using a Michaelis–Menten model. Supplementary Fig. 6, however, shows that both kinetic parameters (K_M and k_cat) change as a function of [Nd³⁺]₀, and led us to the conclusion that the Michaelis–Menten approach could not provide a straightforward explanation for the observed nonlinear effect. Alternatively, a first-order degradation of Tr was also investigated but Supplementary Fig. 7 showed that such degradation step could not explain the nonlinear effect either. Considering these limitations, we developed a mathematical model based on the set of ordinary differential equations according to the processes in Fig. 3. Arithmetic equations are incorporated to account for the binding of the ions to the proteins. Details of the mathematical model is appended to Supplementary Methods Mathematical modelling section.

Nonlinear effect of Nd³⁺ in flow

Next, we establish control over trypsin concentration in time under flow conditions. Figure 4a depicts a flow experiment wherein we changed the [Nd³⁺]₀ (from 0 to 0.7 mM and back in steps of 0.05 mM per 30 min) and monitored [Tr] continuously. Essentially, the residence time in the CSTR, τ (determined by the volume of the CSTR divided over the total flow rate), was comparable to our batch experiments with a reaction time of 6 min. In accordance with our batch experiments, a maximum in [Tr] could indeed be reached when the input satisfied the condition 0.20 ≤ [Nd³⁺]₀ ≤ 0.40 mM. Outside this range, we found that [Tr] was significantly lower. Figure 4b summarises how the use of τ as a flow parameter could tune the behaviour of the system, characterised by an apparent non-equilibrium steady state, [Tr]*, which we determined as the average [Tr] over the last 5 min before [Nd³⁺]₀ was changed. The use of a higher residence time (τ = 10 min) could maintain the nonlinear response but with an elevate values of [Tr]. Hence, while a change in the flow parameter changes the absolute concentrations in [Tr]*, the nonlinear effect of Nd³⁺ on Tr autocatalysis remains.

a Influence of [Nd³⁺]₀ on the autocatalytic conversion of Tg into Tr in a CSTR. [Nd]₀ is changed at regular intervals, programmed by changing [Nd³⁺]_feed (see Methods section). b Nonlinear control over Tr autocatalysis using [Nd³⁺]₀. [Tr]* is the average [Tr] determined based on the time series in a. Error bars are based on the average [Tr] over the last 5 min before [Nd³⁺]₀ was changed. Time series based on τ = 10 min is appended to Supplementary Fig. 13. c Influence of [Nd³⁺]₀ on the autocatalytic conversion of Tg into Tr in a CSTR, in the presence of an additional metal ion (Ca²⁺ and La³⁺). d Effect of additional metal ions on the [Nd³⁺]₀-induced nonlinearity, shown by mapping the dependence of [Tr]* onto polynomial functions with different degrees (indicated with n). Simulated data (bold lines) are appended to Supplementary Fig. 8. Initial conditions: [TRIS-HCl]₀ = 60–100 mM (pH 7.8), T = 21–23 ^oC, [Tg]₀ = 100 µM, [Tr]₀ = 1 µM. For a residence time, τ = 10 min. For c, d residence time, τ = 5 min.

We increase the level of control by forcing competitive binding^32,33,34 among several ions. Figure 4c illustrates our approach wherein five different values of [Nd³⁺]₀ was used to change the rate of Tr autocatalysis in the presence of a second ion (La³⁺ or Ca²⁺). The addition of a slower catalysing ion [Ca²⁺]₀ (20 mM) resulted in the same nonlinear response as in the absence of a second ion (black line) but changed in the position of its local maximum (dashed line). That is, the maximum [Tr]* remained around 35 µM but moved towards slightly higher [Nd³⁺]₀ compared to the systems without Ca²⁺. Similarly, the addition of the faster acting catalyst [La³⁺] also resulted in the same nonlinear response but led to a change in the position of its local maximum into the direction of a lower [Nd³⁺]₀ (dotted line). We used our model to confirm the observation that the local maximum can shift in different directions in the presence of the metal ions (Supplementary Fig. 8).

The opposing effects that the mixtures of ions induce can be mapped onto polynomial functions with different degrees (Fig. 4d). The relation between [Tr]* and [Nd³⁺]₀ can be approximated with a quadratic function in the absence of an additional ion, and the degree of the polynomial function is therefore 2. In the presence of Ca²⁺, this relation changes but only at higher concentrations of [Nd³⁺]₀ and therefore the nature of the polynomial function, its degree, does not change and remains 2. In the presence of La³⁺, however, the behaviour can change from a quadratic to a linear relation, and the degree of the polynomial function becomes 1. An experiment with both La³⁺ and Ca²⁺ shows that the competition for different binding sites allows for a combination of two opposing effects, creating a mechanism to maintain [Tr]* within the narrow range of [Nd³⁺]₀. In the presence of Ca³⁺ and La³⁺, the behaviour became close to constant (i.e., a polynomial function with the degree 0). The mapping of a chemical output onto polynomial functions as demonstrated here underscores the flexibility of using metal ions in controlling the autocatalytic network.

Nd³⁺-induced temporal logic operations

Next, we developed a procedure to obtain Boolean functions for [Tr]* by varying only [Nd³⁺]₀. Figure 5a illustrates our concept wherein three different values for [Nd³⁺]₀ of a given range (a minimum, maximum and a value in between) correspond to four binary combinations (0|0; 0|1; 1|1; 1|0). The notation [Nd³⁺]_A|B is used to represent the input sequence [Nd³⁺]_0|0; [Nd³⁺]_1|0; [Nd³⁺]_1|1; [Nd³⁺]_0|1. As an example, Fig. 5b shows the output sequences for the gates NAND and XOR. We chose the not-AND (NAND) and exclusive-OR (XOR) gates because they are exemplary for functional completeness³⁵ and nonlinear classification³⁶, and are considered difficult to construct using chemical systems^37,38. The differences between the two gates reside on the range of [Nd³⁺]₀ that was applied: the maximum [Nd³⁺]₀ in both cases is 0.6 mM but minimum [Nd³⁺]₀ is 0.2 and 0 mM, respectively. The top graph in Fig. 5b, shows five conditions wherein [Nd³⁺]₀ is changed in steps of 45 min. The first four conditions represent the input 0|0; 0|1; 1|1; 1|0, with the corresponding response in Tr for the two sequences demonstrated in the below graph. A threshold [Tr] with an arbitrary value (20 µM), was used to determine if the observed values in the sequence were true, ‘1’ or false, ‘0’. For the first sequence, [Tr]* exceeded the threshold (red line) at all conditions, except for [Nd³⁺]_1|1. Hence, the input combination resulted in an output sequence representing the NAND-gate: [Tr]_0|0 = 1; [Tr]_0|1 = 1; [Tr]_1|0 = 1; [Tr]_1|1 = 0. For XOR, [Tr]* exceeded the threshold at two conditions, namely for [Nd³⁺]_0|1 and [Nd³⁺]_1|0, yielding the sequence representing the XOR gate: [Tr]_0|0 = 0; [Tr]_0|1 = 1; [Tr]_1|0 = 1; [Tr]_1|1 = 0. The initial value was recovered in the final condition ([Tr]* at [Nd³⁺]_0|0), demonstrating that the output signal can be reset.

a Conceptual scheme for implementing Boolean logic functions. b Experimental time series of a NAND and a XOR gate. Top panel: input sequence encoded in [Nd³⁺]_A|B. Bottom panel: System’s response observed as [Tr](t). The time series for XOR is based on an average from a triplicated experiment (n = 3). c Boolean functions determined as a function of an input sequence of [Nd³⁺]₀ at τ=5 min. [Tr]* is the average [Tr] determined based on the time series in Supplementary Fig. 9. Data points, and error bars, are based on the average over the last 5 min of each condition. d Simulated phase space wherein logic gates can be defined. The ratio between steady states of [Tr]_0|0, [Tr]_0|1 and [Tr]_1|1 determines the type of gate (Supplementary Methods Criteria for distinguishing logic gates section): AND (white circle), OR (white square), XOR (white star), NAND (white hexagon) and NOR (white triangle). The nature of the gate changes when τ=5 min is changed to 10 min (for time series, see Supplementary Fig. 11): from AND to OR (red circle); from XOR to OR (red star). Initial conditions for b and c: [TRIS-HCl]₀ = 100 mM (pH 7.8), T = 21–23 ^oC, [Tg]₀ = 100 µM, [Tr]₀ = 1 µM, residence time, τ = 5 min.

Other gates (AND, OR, NOR) can be developed using a similar procedure (Supplementary Fig. 9). Figure 5c summarises these experiments by plotting the output concentration as a function of [Nd³⁺]₀. The results underline that a change in the minimum and maximum of [Nd³⁺]₀ is sufficient to create different logic operations. We, thus, showed that the autocatalytic reaction can accept instructions (in this case, encoded as an input sequence based on [Nd³⁺]) to perform a range of logic operations.

We used our mathematical model to simulate the response of this autocatalytic network based on changes in the residence time and the concentration of the metal ion. A procedure was developed to discriminate between the different types of logic. Briefly, the ratio [Tr]_1|1:[Tr]_1|0 was used to distinguish AND, OR and XOR and the ratio ([Tr]_1|1 + [Tr]_1|0):[Tr]_0|0 was used to distinguish NAND and NOR (Supplementary Methods Criteria for distinguishing logic gates section). Figure 5d depicts the simulated phase space of the autocatalytic network and shows that the regions of AND, NAND, NOR, XOR and OR can be found in a narrow range of parameters and depending on the [Nd³⁺]_0|0 the regions can overlap. Each grid in the phase space represents a predicted response in [Tr]_ss as a function of input sequence (0|0; 0|1; 1|1; 1|0), and [Nd³⁺]_1|1 is placed on the Y axis. [Nd³⁺]_0|0 on Fig. 5d for AND, OR, XOR equals 0 mM and for NAND and NOR equals 0.2 mM. The phase plot shows that a parameter window exists to maintain each logic operation, ensuring that the desired functions of the system are robust³⁹. This robustness, however, differs per operation. A comparison between the simulated areas for, and experimentally observed, logic gates (open symbols) demonstrates that the model is in good agreement with the experiments. The method of creating logic gates is not restricted to Nd³⁺, and ions can be combined to change the phase space (Supplementary Fig. 10). More importantly, we show that lowering the residence time during the experiment (indicated with the red symbols) could change the output sequence from a XOR or an AND into an OR gate, demonstrating that the Boolean functions can be readily tuned (Supplementary Fig. 11).

Increasing the control over temporal logic operations

The use of an external component such as an inhibitor for thresholding provides further control over the network properties (Fig. 6a). Without the inhibitor, the low steady state was stabilised by continuously depleting Tr using flow, which leads to an unstable steady state. The incorporation of the inhibitor enables a stable and low steady state (as Tr becomes inactive) but does not change the stability of the high steady state, which is supported by autocatalytic production of Tr³⁰. Figure 6b shows a sequence wherein the [Nd³⁺]₀ was changed from a low to a high value, and vice versa, in the presence of a trypsin inhibitor (I, at a concentration of 8, 4, and 2 µM). Overall, [Tr] changes from a high to a low concentration when [Nd³⁺]₀ was changed from 0.3 to 0.8 mM. At an inhibitor concentration that appeared to be too high ([I]₀ = 8 µM), we observed that original [Tr]* cannot be recovered and instead remains low when [Nd³⁺]₀ was changed in the opposite direction (from 0.8 to 0.3 mM). The observation that a pulse in [Nd³⁺]₀ (i.e., the sequence 0.3–0.8–0.3 mM) enabled a retention of the low [Tr]* demonstrates that the autocatalytic network, in the presence of the inhibitor, is capable of a so-called long-term depression of the input. This effect disappeared when the inhibitor concentration was lowered to [I]₀ = 4 µM but notably [Tr]* remains to depend on the amount of inhibitor present in its preceding state, I.e., the response is history dependent. Validations in Supplementary Fig. 12 supports the observation that the XOR gate is (and perhaps other logic gates are) history dependent. The history dependency was lost when [I]₀ = 2 µM was applied and, instead, the response became bistable. That the signal could fully recover, however, shows that a trace of inhibitor allows for an increased robustness of the output signal. Hence, while the metal ion provides control for tuning and switching of the logic operation, an inhibitor was introduced as an additional control parameter for strengthening or weakening its history dependency.

a Schematic representation of autocatalysis controlled by the metal ion, Nd³⁺, and Soybean Trypsin inhibitor, I. b Experimental time series with [I]₀ = 8, 4 and 2 µM. Initial conditions: [TRIS-HCl]₀ = 100 mM (pH 7.8), T = 21–23 ^oC, [Tg]₀ = 100 µM, [Tr]₀ = 1 µM, 20 mM CaCl₂, τ = 4.5 min.