Autocatalytic, Bistable, Oscillatory Networks of Biologically Relevant Organic Reactions
Sergey N. Semenov,a Lewis J. Kraft,a Alar Ainla,a Mengxia Zhao,a Mostafa Baghbanzadeh,a Victoria E. Campbell,a Kyungtae Kang,a Jerome M. Fox,a and George M. Whitesidesa,b,c* a Department of Chemistry and Chemical Biology, Harvard University 12 Oxford Street, Cambridge, MA 02138 b Kavli Institute for Bionano Inspired Science and Technology, School of Engineering and Applied Sciences, Harvard University 29 Oxford Street, Cambridge, MA 02138 c Wyss Institute for Biologically Inspired Engineering, Harvard University 60 Oxford Street, Cambridge, MA 02138
* Author to whom correspondence should be addressed.
	 1	

Networks of organic chemical reactions are centrally important in life, and were likely to have played a central role in its origins.1-3 Network dynamics regulate cell division,4-6 circadian rhythms,7 nerve impulses,8 chemotaxis,9 and guide development of organisms.10 Although out-of-equilibrium networks of chemical reactions have the potential to display emergent network dynamics11 such as spontaneous pattern formation, bistability, and periodic oscillations,12-14 the principles that enable networks of organic reactions to develop complex behaviors are incompletely understood. Here we describe a network of biologically relevant organic reactions (amide formation, thiolate-thioester exchange, thiolate-disulfide interchange, and conjugate addition) that displays bistability and oscillations in concentrations of organic thiols and amides. Oscillations arise from the interaction between three subcomponents of the network: (i) an autocatalytic cycle that generates thiols and amides from thioesters and dialkyl disulfides; (ii) a trigger that controls autocatalytic growth; and (iii) inhibitory processes that remove activating thiol species produced during the autocatalytic cycle. In contrast to previous studies demonstrating oscillations and bistability using highly evolved biomolecules (i.e., enzymes15 and DNA16,17) or inorganic molecules of questionable biochemical relevance (e.g. those used in Belousov-Zhabotinsky-type reactions),18,19 the organic molecules used in our network are relevant to current metabolism and similar to those that might have existed on early Earth. By using small organic molecules to build a network of organic reactions with autocatalytic, bistable, and oscillatory behavior, we identified principles that clarify how dynamic networks relevant to life might possibly have developed. In the future, modifications of this network will clarify the influence of molecular structure on the dynamics of reaction networks, and may enable the design of biomimetic networks, and of synthetic self-regulating and evolving chemical systems.
	 2	

Figure 1 summarizes the network of organic reactions that we used to assemble our model system. All of these reactions are nearly quantitative, and the structure of their reactants can be varied by synthesis to control rates of reactions. Thiols and thioesters, which are central to these reactions, are important in many biological processes (e.g., the formation of disulfide bonds in proteins, transformations involving coenzyme-A, the reduction of oxidized molecules by glutathione,20 the synthesis of polyketides,21 and the nonribosomal synthesis of peptides21), and thus, might represent reactions that enabled life to emerge on early Earth.22
To control the dynamics of these processes, we constructed a modular chemical reaction network (CRN) shown schematically in control-theoretic terms23 in Fig. 1a. A “trigger” produces an initial chemical signal, and an “auto-amplifier” amplifies this signal, which may or may not be suppressed by inhibition. To keep the reactions out of equilibrium—and thus, to enable the self-organization of reactions by communication through concentrations of reactants and products — we used a continuous-stirred tank reactor (CSTR) to mix reactants and products, while allowing a flux of species into and out of the network over time. A biological cell has some conceptual analogies to a micron-scale, diffusively mixed, tank reactor. The dynamic behavior of this system — especially bistability and oscillations — reflects the balances of triggering, auto-amplification, and inhibition.
We first constructed a chemical network capable of auto-amplification using thiols and thioesters (Fig. 1b). The starting components of the network are cystamine (CSSC, 3) and L-alanine ethyl thioester (AlaSEt, 2). Trace amounts of cysteamine (CSH, 1) are generated as follows: AlaSEt slowly hydrolyzes, generating alanine (8) and ethanethiol (ESH, 4); EtSH then reacts with CSSC via thiolate-disulfide interchange,24 yielding disulfide 6 and CSH. With CSH present, self-amplification occurs through two steps: (i) CSH reacts with AlaSEt rapidly by thiolate-thioester exchange and intramolecular rearrangement (Kent
	 3	

ligation),25 yielding two thiols, EtSH, and L-alanine mercaptoethyl amide (5); and (ii) these two thiols then undergo thiolate-disulfide interchange with CSSC, yielding two molecules of CSH (and disulfides 6 and 7). Because one molecule of CSH generates two molecules of CSH (a 2X multiplication), this network proceeds autocatalytically (hereafter, we refer to it as a “autocatalytic thiol network”).
We tested this autocatalytic thiol network in a batch reactor (in 500-mM phosphate buffer at pH 7.5). The curve showing the formation of alanine amides (5, 7, 11, 12, 13) from the reaction of AlaSEt and CSSC has two attributes that are consistent with autocatalysis (Fig. 2a): (i) it is sigmoidal in shape, and (ii) the addition of β-mercaptoethanol eliminates the lag period. The results of numerical simulations agreed well with experimental results (Extended data Fig. 1e and Methods). Production of thiols (CSH, 5, and EtSH) also follows autocatalytic kinetics (Extended Data Fig. 1f). Thus, this system provides a kinetically simple, well-defined autocatalytic reaction—with reactants that can be easily varied structurally by synthesis—based on thiols. Autocatalytic reactions are important in life and play an important role in evolution.26
With an autocatalytic reaction system in hand, we focused on developing a set of reactions to control the onset of exponential growth (the timing of the trigger). The mechanism outlined in Fig. 1b suggests that a molecule capable of rapidly reacting with—and thus, rapidly removing—free thiols might delay autocatalytic growth until that molecule is itself depleted by reaction. We, therefore, hypothesized that the addition of maleimide, which reacts extremely rapidly with thiolates (in comparison to all other rates of reaction in the autocatalytic thiol network), would prevent the formation of 5—and thus prevent autoamplification of thiols—until it is completely consumed. Therefore, in a reaction of AlaSEt, CSSC, and maleimide, molecules of EtSH, a free thiol generated by the hydrolysis of
	 4	

AlaSEt, will react with—and, eventually, consume—maleimide, at which point the autocatalytic formation of CSH can proceed (Fig. 1b).
To test the ability of maleimide to delay autocatalytic growth, we added 5 mol % of this molecule to a mixture of CSSC and AlaSEt (in batch), and used proton nuclear magnetic resonance (1H NMR) to monitor both the production of alanine amides and the disappearance of maleimide by (Fig. 2b and Extended Data Fig. 3a-c). The production of alanine amides does not proceed until maleimide is consumed, at which point it proceeds autocatalytically.
In a CSTR, autocatalytic reactions can become bistable.27 In flow systems, bistability often manifests itself through hysteresis; concentrations of reaction intermediates, which are determined by flow rates, do not return to their original values when flow rates are changed and then changed back. To test our reaction network under flow conditions, we constructed a small CSTR and connected it to a spectrophotometer to monitor the total concentration of thiols in the system continuously (Fig. 3, Extended Data Fig. 4). The molecules fed into this CSTR included AlaSEt, CSSC, and maleimide.
Numerical simulations suggested that space velocities (SV = flow rate/reactor volume, s-1) in the range of 0.0001-0.01 s-1 would produce hysteresis (see Extended Data Fig. 5a-c). To test the result of these simulations, we monitored the total concentration of thiols during stepwise changes in SV within the predicted range. We started from a low flow rate, ramped up to a high flow rate, and then returned to the low flow rate (Fig. 4a). High thiol concentrations generated through self-amplification (of CSH) require the SV to be lowered to 0.0005 s-1 or below to activate the autocatalytic pathway. The system transitions out of the self-amplifying state when the SV is ≥ 0.006 s-1 (Fig. 4b). These limits make qualitative sense: self-amplification requires maleimide to be removed from the CSTR (via reactions with free thiols – more exactly, thiolate anions – or transport out of the CSTR through the outlet port) more rapidly than it is added through the inlet port; the termination of self-
	 5	

amplification, once it has started, requires free thiols to be removed from the CSTR by transport out of the outlet port more rapidly than they are produced (via autocatalysis). An increase in maleimide concentration lowers the limiting flow rates for bistability (Fig. 4c; Extended Data Fig. 5d-f), as the model predicts. This chemical reaction network illustrates a general way to convert any quadratic autocatalytic system (and perhaps, any system capable of exponential amplification) into a bistable switch by combining this autocatalysis with linear generation and quantitative sequestration of the self-amplifying species.
Epstein, De Kepper, and coworkers showed that bistable systems can give rise to oscillations in the presence of an inhibition reaction that slowly (relative to a production reaction) removes the self-amplifying species.13,27,28 In our system, we chose acrylamide as an inhibitor, because this molecule reacts more slowly with thiols than any other molecule in the reaction network. When we tested this system with acrylamide in batch, it exhibited a single oscillation (i.e., one peak) in the concentration of free thiols (Fig. 2c; Extended Data Fig. 3d-e). NMR analysis (Fig. 2d) shows that the oscillation is not triggered until the maleimide is removed. In addition, it is apparent that AlaSEt is depleted during the oscillation.
Using a combination of numerical simulations and experiments in the CSTR under different flow rates (inputs were AlaSEt, maleimide, and acrylamide), we identified conditions under which the addition of acrylamide, present in excess, would produce sustained oscillations in RSH (Fig. 4d; Extended Data Fig. 6). Each oscillation arises from three steps: (i) the trigger sets up a delay for auto-amplification of the autocatalytic thiol species; (ii) autocatalytic growth competes with inhibition of thiols by reaction with acrylamide, thus depleting the thioester, AlaSEt; (iii) the thioester, AlaSEt, (and other reagents) is recharged by flow into the CSTR; (Fig. 4e and Methods).
	 6	

To determine how changes in flow rate effect oscillations, we examined the influence of flow rate on the stability, period, and amplitude of oscillations. Sustained oscillations occurred for values of SV between 1.4 x 10-4 and 2.7 x 10-4 s-1 (Fig. 4f; Extended Data Fig. 7). Increasing the SV to 2.8 x 10-4 s-1 resulted in a sudden loss of oscillations, and decreasing the SV to 1.2 x 10-4 s-1 resulted in damped oscillations. The period increased nonlinearly with SV (Fig. 4f); the amplitude, by contrast, increased linearly (Fig. 4g).
To explain the trends in period and amplitude of the oscillations, and the nature of the bifurcations at low and high limiting values of SV, we constructed a simple kinetic model (eq. 1-3, and diagram on Fig. 4h) that enables qualitative analysis of dynamic behaviors in our oscillating network. This model simplifies the autocatalytic thiol network to bimolecular autocatalytic production of thiols from thioester, and considers [CSSC] and [acrylamide] as constants.
In this model, A corresponds to [RSH], I to [Maleimide], S to [AlaSEt], k1-4 to rate constants, and k0 to space velocity. Linear stability analysis13 carried out with this model shows that increasing k0 from low to high values causes two transitions: first, the system transitions from one having a stable focus (damped oscillations) to one having a stable orbit (sustained oscillations); this transition marks an Andronov-Hopf bifurcation.29 Second, the system transitions from one having a stable orbit to one having a single stable equilibrium “node” (an equilibrium in which the system, when perturbed, returns directly to the stable state, as opposed to orbiting around it); this transition marks a saddle-node or fold bifurcation29 (see Extended Data Fig. 8a and Methods). Further analysis of this model shows that lowering the concentration of maleimide (I0) fed into the reactor results in a transition from oscillatory to bistable behavior (see phase plot in Extended Data Fig. 8b). Indeed, experi	 7	

ments with [Mal]0 = 4 mM (all other parameters are the same as those in experiments on oscillations) demonstrated the predicted bistable behavior of the system at space velocities from 0.0015 to 0.005 s-1 (Extended Data Fig. 8c).
Our experiments show that we have successfully converted a bistable organic reaction network into an oscillatory network with two interesting, experimentally accessible—and biologically relevant—classes of bifurcations. Systematically designed chemical oscillators are rare. Notable examples are enzymatic,15 pH,30 DNA,16,17 and redox oscillators;28 here, we have demonstrated a new class of chemical oscillators based on thiols. This system is unique in that it has characteristics that make it relevant to the discussion of the spontaneous formation of networks from types of molecules relevant to metabolism, to other dynamic networks of organic reactions, and perhaps to the origin of life.
This network, which requires only four components (a thioester, a disulfide of a βaminothiol, a kinetically rapid and irreversible thiol inhibitor, and a slow and irreversible thiol inhibitor), can be varied through structural modifications of the constituent organic molecules; these flexibilities in structures offer broad potential to explore the evolution and functional design of dynamic networks. The complex behavior of this simple network of organic reactions illustrates four ideas important in considering the emergence of complexity in dynamic networks of chemical reactions: (i) They demonstrate the ability of a simple set of reactants to produce emergent behaviors when feedback is present. (ii) They show “memory,” in the sense that the system is hysteretic, and thus codes information about its past in the concentrations of molecules in the network; “memory” in discussions of origin of life, is usually attributed to sequences of bases in replicating molecules of RNA, or its precursors, and alternative forms of memory in dynamic networks broaden its possible origins. (iii) It is compatible with the sorts of behaviors envisioned by De Duve in his proposal for a “thioester” basis for the origin of life, and this compatibility reinforces the relevance of this
	 8	

proposal; (iv) It suggests that the complex regulatory dynamics generally associated with cellular function (e.g. threshold-based behaviors, circadian rhythms) may, in principle, have arisen spontaneously before cellular life and—perhaps, most importantly—before enzymatic catalysis.
Online Content Methods, along with any additional Extended Data display items, are available in the online version of the paper; references unique to these sections appear only in the online paper.
Acknowledgements This work was funded in part by the Simons Foundation under award 290364, and the
Templeton Foundation under award 48423.
Author Information The authors declare no competing financial interests.
Author Contributions S.N.S. and G.M.W. conceived the research and designed the experiments. S.N.S.,
L.J.K., A.A., M.Z., V.E.C., M.B., and K.K. performed experiments and analyzed data. S.N.S., L.J.K., A.A., and J.M.F. performed computational simulations. S.N.S., L.J.K., J.M.F, V.E.C., and G.M.W. wrote the manuscript.
	 9	

References
1. Dyson, F. J. A model for the origin of life. J. Mol. Evol. 18, 344-350 (1982). 2. Nghe, P. et al. Prebiotic network evolution: Six key parameters. Mol. Biosyst. 11,
3206-3217 (2015). 3. Patel, B. H., Percivalle, C., Ritson, D. J., Duffy, C. D. & Sutherland, J. D. Common
origins of RNA, protein and lipid precursors in a cyanosulfidic protometabolism. Nat. Chem. 7, 301-307 (2015). 4. Tyson, J. J., Chen, K. C. & Novak, B. Sniffers, buzzers, toggles and blinkers: Dynamics of regulatory and signaling pathways in the cell. Curr. Opin. Cell. Biol. 15, 221-231 (2003). 5. Ferrell, J. E., Tsai, T. Y. C. & Yang, Q. O. Modeling the cell cycle: Why do certain circuits oscillate? Cell 144, 874-885 (2011). 6. Tyson, J. J. Modeling the cell-division cycle - cdc2 and cyclin interactions. P. Natl. Acad. Sci. USA 88, 7328-7332 (1991). 7. Goldbeter, A. A model for circadian oscillations in the drosophila period protein (PER). Proc. Biol. Sci. 261, 319-324 (1995). 8. Fitzhugh, R. Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445-466 (1961). 9. Laub, M. T. & Loomis, W. F. A molecular network that produces spontaneous oscillations in excitable cells of dictyostelium. Mol. Biol. Cell 9, 3521-3532 (1998). 10. Lander, A. D. Pattern, growth, and control. Cell 144, 955-969 (2011). 11. Hazen, R. M., Griffin, P. L., Carothers, J. M. & Szostak, J. W. Functional information and the emergence of biocomplexity. P. Natl. Acad. Sci. USA 104, 8574-8581 (2007). 12. Whitesides, G. M. & Ismagilov, R. F. Complexity in chemistry. Science 284, 89-92 (1999). 13. Epstein, I. R. & Pojman, J. A. An introduction to nonlinear chemical dynamics : Oscillations, waves, patterns, and chaos. (Oxford University Press, 1998). 14. Grzybowski, B. A. Chemistry in motion : Reaction-diffusion systems for micro- and nanotechnology. (Wiley, 2009). 15. Semenov, S. N. et al. Rational design of functional and tunable oscillating enzymatic networks. Nat. Chem. 7, 160-165 (2015). 16. Kim, J. & Winfree, E. Synthetic in vitro transcriptional oscillators. Mol. Syst. Biol. 7, 465 (2011). 17. Montagne, K., Plasson, R., Sakai, Y., Fujii, T. & Rondelez, Y. Programming an in vitro DNA oscillator using a molecular networking strategy. Mol. Syst. Biol. 7, 466 (2011). 18. Belousov, B. P. Periodically acting reaction and its mechanism. Sbornik Referatov po Radiatsionni Meditsine, 145 (1958). 19. Gyorgyi, L., Turanyi, T. & Field, R. J. Mechanistic details of the oscillatory belousovzhabotinskii reaction. J. Phys. Chem.-US 94, 7162-7170 (1990). 20. Meister, A. & Anderson, M. E. Glutathione. Annu. Rev. Biochem. 52, 711-760 (1983). 21. Fischbach, M. A. & Walsh, C. T. Assembly-line enzymology for polyketide and nonribosomal peptide antibiotics: Logic, machinery, and mechanisms. Chem. Rev. 106, 3468-3496 (2006). 22. De Duve, C. Singularities landmarks on the pathways of life. (Cambridge University Press, 2005).
	 10	

23. LeDuc, P. R., Messner, W. C. & Wikswo, J. P. How do control-based approaches enter into biology? Annu. Rev. Biomed. Eng. 13, 369-396 (2011).
24. Houk, J. & Whitesides, G. M. Structure reactivity relations for thiol disulfide interchange. J. Am. Chem. Soc. 109, 6825-6836 (1987).
25. Dawson, P. E., Muir, T. W., Clarklewis, I. & Kent, S. B. H. Synthesis of proteins by native chemical ligation. Science 266, 776-779 (1994).
26. Bissette, A. J. & Fletcher, S. P. Mechanisms of autocatalysis. Angew. Chem. Int. Ed. 52, 12800-12826 (2013).
27. Boissonade, J. & De Kepper, P. Transitions from bistability to limit-cycle oscillations - theoretical-analysis and experimental-evidence in an open chemical-system. J. Phys. Chem.-US 84, 501-506 (1980).
28. De Kepper, P., Epstein, I. R. & Kustin, K. A systematically designed homogeneous oscillating reaction - the arsenite-iodate-chlorite system. J. Am. Chem. Soc. 103, 2133-2134 (1981).
29. Kuznetsov, I. U. A. Elements of applied bifurcation theory. 3rd edn, (Springer, 2004).
30. Kovacs, K., McIlwaine, R. E., Scott, S. K. & Taylor, A. F. An organic-based pH oscillator. J. Phys. Chem. A 111, 549-551 (2007).
	 11	

Figure 1 | Overview of the network of organic reactions presented in this study. a, A diagram showing inputs and outputs from a continuously stirred tank reactor (CSTR). Input parameters include reactor space velocity (SV), pH, concentrations of chemical species, buffer composition, and temperature. The output is thiol concentration ([RSH]), which changes with time. Inside the CSTR, the trigger is a balance between the generation and the destruction of RSH, which regulates the autoamplification of RSH. Inhibition and washout remove RSH. b, Reaction schemes describing the elements of the diagram in a. 1. Trigger: slow production of RSH by hydrolysis of the thioester, in combination with removal of RSH by reaction with maleimide. 2. Autoamplification by autocatalysis: RSH provides a positive feedback loop in our system by thiol-disulfide exchange with CSSC, followed by native chemical ligation with AlaSEt (direct attack of CSSC on AlaSEt is negligible; see Extended Data Fig. 1a-d). 3. Inhibition: removal of RSH by conjugate addition to an excess of acrylamide provides a negative feedback loop. 4. Washout and refill of CSTR: AlaSEt, CSSC, maleimide, and acrylamide, are supplied from syringes; all compounds are removed as effluent (see Extended Data Fig. 2 for omitted exchange reactions and Extended Data Table 1. for details of protonation equilibria).
	 12	

Figure 2 | Summary of the batch kinetic studies of the organic reaction network. a, Experiments showing elimination of the lag period in the autocatalytic reaction (AlaSEt with CSSC) by addition of β-mercaptoethanol (ME). To quantify the kinetics of this reaction, we followed the production of amide derivatives of alanine using 1H NMR spectroscopy; the chemical shift of a proton linked to the α-carbon of L-alanine is different for AlaSEt (4.05 ppm) and alanine amides (5, 7, 11, 12, 13); (the chemical shifts, 3.85 ppm, of these five compounds are indistinguishable). Reaction conditions: D2O, 500-mM phosphate buffer pH 7.5, 25 oC. b, To quantify the kinetics of the trigger (the reaction of AlaSEt with CSSC in the presence of maleimide), we again followed the production of amide derivatives of alanine and maleimide by 1H NMR. Region I, a period of delay regulated by the trigger; region II, a period of autocatalytic production of RSH; and region III, AlaSEt is depleted and thus no further reactions occur. See Extended Data Fig. 3a-c for additional experiments. Reaction conditions: same as in a. c, In order to quantify the kinetics of a single oscillation in [RSH] (the reaction of AlaSEt, CSSC, maleimide, and acrylamide), we monitored [RSH] using Ellman’s reagent (Extended Data Fig. 4). Regions I and II correspond to the trigger and autocatalytic regimes, respectively; region III corresponds to a period of inhibition by conjugate addition of RSH to acrylamide. Reaction conditions: H2O, 1 M phosphate buffer pH 8.0, 25 °C. d, Time dependence of concentrations of AlaSEt, CSH, EtSH, and maleimide during a single oscillation by 1H NMR. Reaction conditions: D2O, 1 M phosphate buffer pH 8.0, 25 °C, [AlaSEt] = 47 mM; [CSSC] = 92 mM; [Maleimide] = 10 mM; [Acrylamide] = 320 mM.
	 13	

Figure 3 | Schematic representation of the CSTR experimental setup. In this setup, syringe pumps feed reactants to the inlet ports of the CSTR, where they mix and react. A tube connected to the outlet port of the CSTR takes the effluent (products and reactants) to a detection chip, where RSH is derivatized and detected by UV-Vis spectroscopy.
	 14	

Figure 4 | The network of organic reactions displays bistability and sustained oscillations under flow. a, Representative experimental trace showing changes in [RSH] over time as a result of changes in SV (dashed lines). b, Hysteresis curve derived from steady-state [RSH] in the trace presented in a (see Methods for additional conditions). Error bars correspond to standard deviations that were calculated from data points (n > 100) in the time intervals that were used to determine steady state values of [RSH] for each SV. c, Phase plot reconstructed from three separate experiments using different maleimide concentrations. Experimental conditions: H2O, 500 mM phosphate buffer pH 7.5, 25 oC; [AlaSEt] = 46 mM; [CSSC] = 92 mM; [Maleimide] = 1.16, 2.31, and 3.47 mM; SV = 0.25 – 6.10-3 s-1. Gray shading in a-c indicates regions of bistability. (d) Experimental data showing sustained oscillations at two different SV values, 2.7.10-4 s-1 (red line) and 2.0.10-4 s-1 (black line), using the following reactants: [AlaSEt] = 47 mM; [CSSC] = 92 mM; [Maleimide] = 10 mM; [Acrylamide] = 320 mM in 1000 mM phosphate buffer pH 8.0, 28 oC. (e) Simulations that show an example of the evolution of the concentration of key species (RSH, maleimide, AlaSEt) during oscillations (see Extended Data Fig. 6a for modeling parameters). (f) Influence of the SV on the stability and period of the oscillations. Open circles represent experimental data. The solid line represents a hyperbolic fit to the data. Subpanels show examples of the dynamics of [RSH] in the indicated regions. (g) The influence of SV on the amplitude of oscillations. Open circles represent experimental data. The solid line represents a linear fit to the data. Dashed lines represent estimated borders for sustained oscillations. Error bars correspond to 90 % confidence intervals (3 independent experiments for each SV). Separate plots for each experimental examination of oscillations can be found in Extended Data Figs. 7. (h) Simplified depiction of the overall network of organic reactions. Arrow heads and blunt arrow heads represent activation and inhibition, respectively.
	 15	

Methods
	 16	

Synthesis General: Unless otherwise stated, all starting materials were obtained commercially and
were used without further purification. Nuclear Magnetic Resonance (NMR) spectra were measured on a Varian INOVA I-500 spectrometer at 500 MHz for 1H and 125.8 MHz for 13C(1H) or on Varian INOVA I-400 spectrometer at 400 MHz for 1H. The chemical shifts for 1H and 13C are given in parts per million (ppm) relative to TMS, and calibrated using the residual 1H peak of the solvent; δ = 7.26 for CDCl3 and δ = 4.79 for D2O in 1H NMR and δ = 77.2 for CDCl3 in 13C NMR. Fourier transform infrared spectroscopy (FT-IR) spectra were recorded on a Bruker TENSOR 27 spectrometer fitted with an attenuated total reflectance (ATR) cell.
N-Boc-L-alanine ethyl thioester (BocAlaSEt): N-(3-Dimethylaminopropyl)-N′ethylcarbodiimide hydrochloride (EDC•HCl), (5.57 g, 29.0 mmol) was added to a stirred solution of N-Boc-L-alanine (4.00 g, 26.4 mmol) and ethanethiol (2.45 mL, 34.4 mmol) dissolved in 25 mL of DMF at 0°C under an argon atmosphere. After allowing the reaction to proceed for 1 hour at 0°C, the mixture was left at room temperature overnight. DMF was removed in vacuo. The residue was dissolved in ethyl acetate (50 mL), washed with DI water (2 x 50 mL), dried over Na2SO4, and concentrated on the rotary evaporator. The crude product was purified by column chromatography (SiO2, CH2Cl2/ EtOAc 95:5, Rf ~ 0.5). Yield: 63 % (3.90 g, 16.74 mmol) 1H NMR: (400 MHz, CDCl3) δ = 4.38 (m, 1H, NHCH), 2.87 (q, 3JH-H = 7.4 Hz, 2H, SCH2), 1.45 (s, 9H, C(CH3)3), 1.37 (d, 3H, 3JH-H = 7.2 Hz, NHCHCH3) 1.25 (t, 3JH-H = 7.4 Hz, 2H, SCH2CH3).
L-alanine ethyl thioester trifluoroacetic salt (AlaSEt•TFA): Trifluoroacetic acid (TFA) (12 mL) was added dropwise to a solution of N-Boc-L-alanine ethyl thioester (1.25 g 5.4 mmol) in dichloromethane (CH2Cl2) (12 mL) at 0°C. The reaction was stirred at 0°C for
	 17	

15 min and then at room temperature for 45 min. The solvents were removed in vacuo, and the residue was dissolved in toluene (5 mL). Removal of toluene in vacuo yielded pure Lalanine ethyl thioester trifluoroacetic salt (AlaSEt•TFA) as a white solid (96% pure, TFA is a major admixture). Its purity was determined by NMR with L-alanine as an internal standard. Yield: 97 % (1.30 g, 5.3 mmol) 1H NMR: (500 MHz, D2O) δ = 4.18 (q, 3JH-H = 7.5 Hz, NHCH), 2.85 (m, 2H, SCH2), 1.43 (d, 3H, 3JH-H = 7 Hz, NHCHCH3) 1.11 (t, 3JH-H = 7.5 Hz, 3H, SCH2CH3); Batch kinetic experiments
General: The kinetic experiments were performed at 25°C in a 500 - mM phosphate buffer (pH 7.5) solution in D2O by monitoring the change in the 1H NMR spectra. We calculated the progress of the Kent ligation reaction by integrating the 1H NMR signal of the proton linked to the α-carbon of L-alanine derivatives at δ = 3.85 ppm (α-protons of 5, 7, 11, 12, 13). The concentrations of the starting materials were calculated by setting the total integral of alanine α-protons equal to the starting concentration of AlaSEt. Total thiol concentration ([CSH]+ [EtSH]+[5]) in the batch experiments was determined with Ellman’s reagent (5,5’-dithiobis-2-nitrobenzoic acid) (0.5 mM, phosphate buffer, pH 7, 100 mM). A 7.5 - µL sample was added to the Ellman’s solution and the absorbance of this solution was measured by UV-visible spectroscopy after 1 min incubation. We calculated thiol concentration from absorbance using a calibration curve (Extended Data Fig. 4a). In all plots, we show total concentration of all protonated/deprotonated forms, see Extended Data Table 1. for details of protonation equilibria. 31-33
Reaction between AlaSEt and cystamine (CSSC). First, we dissolved cystamine hydrochloride (CSSC•2HCl) (40.2, 80.5, or 161.0-mM solution) in 0.4 mL of 500-mM phosphate buffer, pH 7.5, in D2O and filled the NMR tube with this solution. Second, we dissolved AlaSEt•TFA (107 mM) in 0.3 mL of 500-mM phosphate buffer, pH 7.5, in D2O
	 18	

and mixed it with the CSSC•2HCl solution. Third, we recorded the 1H NMR spectrum every minute (72 s time intervals considering 12 s of acquisition time).
Reaction between AlaSEt and CSSC in the presence of maleimide. We performed experiments as described above. We dissolved CSSC•2HCl (80.5 mM) and maleimide (2.0, 4.0, or 6.1 mM) in 0.4 mL of 500-mM phosphate buffer, pH 7.5, in D2O and filled the NMR tube with this solution. Thereafter, we dissolved AlaSEt•TFA (107 mM) in 0.3 mL of 500 mM phosphate buffer, pH 7.5, in D2O and mixed it with the first solution. In addition to the progress of Kent ligation, we monitored the disappearance of maleimide by following the disappearance of maleimide proton signals in the 1H NMR spectra. The interval between measurements was 2 min (132 s time intervals considering 12 s of acquisition time).
Reaction between AlaSEt and CSSC in the presence of maleimide and acrylamide. First, we weighted 8.1 mg of AlaSEt•TFA in a 2 mL glass vial containing a stirring bar. Second, we dissolved CSSC•2HCl (92 mM), of maleimide (10 mM), and acrylamide (160 or 320 mM) in 0.7 mL of 1 M phosphate buffer, pH 8, in H2O. Third, we quickly added the solution to the alanine thioester, shook vigorously, and put it on the magnetic stir plate. Aliquots 7.5 µL of this mixture were used to determine the amount of free thiols at each time point.
Flow experiments We used two types of flow systems: (i) a flow system with the CSTR, mixing, and
detection chip made completely from polydimethylsiloxane (PDMS) (Setup #1), and (ii) a flow system with a hybrid PDMS/glass CSTR, and detection chip with mechanical mixing (Setup #2). Setup #1 is simple to fabricate, and was used for the bistability experiments. Setup #2 is harder to fabricate than Setup #1, but it gave a more stable signal than Setup #1, therefore, it was used for the oscillatory experiments.
	 19	

Setup #1 (for bistability experiments). The CSTR had a volume of 250 µL, was made of PDMS, used a built-in temperature control channel, and had three separate inlets. We mounted three syringes on a Harvard PhD 2000 syringe pump and fed the CSTR with reactants via the inlet tubing. The temperature control channel was connected to a thermostatic water bath. The outlet of the CSTR was connected to the microfluidic detection chip. This detection chip had two inlets, a mixing channel, and a detection chamber that was 150 µm in height. It served two functions: (i) mixing outgoing flow of the CSTR with Ellman’s reagent (5,5’-dithiobis-(2-nitrobenzoic acid)), and (ii) providing an optical window for the UV-Vis absorption measurements. The syringe pump supplied the Ellman’s reagent via the second inlet (detection chip). Ellman’s reagent quantitatively and irreversibly reacted with the free thiols in the reaction mixture releasing an orange colored 2-nitro-5-thiobenzoate that was detected by UV-Vis absorption. To perform UV-Vis absorption measurements, the detection chamber was aligned with collimated light from the Nikon Intensilight C-HGFI and an optic fiber going to the Ocean Optics HP 4000 UV-Vis absorption spectrometer. We detected the intensity of 412 nm light over course of the experiment. We converted absorbance into the thiol concentration using a calibration curve (Extended Data Fig. 4b). Artifacts from air bubbles passing detection window were filtered for clarity. All of the pumps were controlled by a MATLAB program that was specifically designed for these experiments (the program can be found in the supplementary files).
Setup #2 (for oscillatory experiments). Detection of the oscillations in the oscillatory experiments required three changes in the fluidic system: using glass instead of PDMS as the material for the CSTR; using active instead of passive mixing during the derivatization step, and using a glass flow cell coupled to a Cary 60 UV-vis spectrometer. The glass CSTR ensured no leakage of EtSH by diffusion through PDMS, and also prevented gas bubble accumulation. The active mixing decreased the retention time during the derivatization step,
	 20	

thus reducing the effect of a side reaction between alanine ethyl thioester and Ellman’s reagent (we determined experimentally that, in this setup, the background reaction with AlaSEt will generate a signal corresponding to an apparent thiol concentration in the range of 0.4 – 0.8 mM for the flow rates and concentrations of reagents that were used in our experiments; apparent rate constant is about 4.104 times slower than the reported apparent second-order rate constant of the reaction between Ellman’s reagent and β-mercaptoethanol at pH 7).34 The use of a UV-vis spectrometer eliminated any intensity fluctuations from the light source. Absorbance data were converted into total concentration of free thiols using a calibration curve (Extended Data Fig. 4 c, d).
The procedure of making a CSTR with a glass body consisted of three steps. In the first step, we fabricated a glass body. We took a Pasteur pipette, fixed it upside-down in a drill press and heated it 13 mm from the bottom end using a hand hold burner (MT-51). Interplay between gravity and surface tension resulted in the formation of a glass neck with thick walls. The neck had an internal diameter 0.9-1 mm. We cut the pipette 8 mm above the neck. In order to make an outlet connection for tubing, we put the newly made glass chamber upsidedown in a petri dish filled with PDMS, such that PDMS was up to the middle of the neck, and inserted a needle (0.9 mm outer diameter) in to the neck from the top. A well-defined channel was formed after the PDMS was cured. In the second step, we fabricated a PDMS body with all of the inlet connections and temperature controller. We placed a brass cylinder (d = 6.82 mm; h = 6.82 mm) and a copper tube (one turn around the cylinder) in the petri dish, and filled the dish with PDMS up to the top of the cylinder such that the top of the cylinder was not covered by PDMS. After curing, we pulled the PDMS out of the dish and removed the cylinder. We cut the inlet channels on the bottom of the PDMS body, and punched holes for the inlet tubing. On the final, and third step, we plasma-bonded the PDMS body to a glass slide, inserted the glass body into the PDMS body (hole from the brass cylinder), sealed the
	 21	

top with additional PDMS, and inserted tubing. Note, it is important to insert the glass body as close as possible to the glass slide. The elasticity of the PDMS would not, however, allow for full contact between the glass body and glass slide, which left a gap (0.1-0.2 mm) that was exploited to construct inlet channels for the reactor.
To make a fast in-flow mixing chip, we designed the microfluidic device with mixing chamber (diameter: 3mm, height: 1mm, volume: 7 µL) in AutoCAD and 3D printed molding master with an Objet Connex 500 3D printer (Stratasys). After printing, the master was heated in oven at 70°C overnight to complete the photo-curing process (otherwise, unreacted components inhibit the PDMS curing in a later step). Using syringe needles, we defined liquid ports (0.9 mm OD), and the placement for the impeller shaft (0.6 mm OD). Thereafter, we used a standard PDMS molding procedure to cast the chip body. After removal from the master, a syringe needle (0.9 mm OD) was inserted into the mixing chamber, and a small impeller was grafted into its end. Then the PDMS body was plasma bonded to a glass substrate, and the syringe needle was sealed with small drop of PDMS. We connected the shaft of the impeller to a 5V micro-DC motor running at a speed of about 600 rpm. Efficient mixing in the outlet channel of the mixing chip was qualitatively confirmed by visual inspection using water and dye.
To construct the glass flow cell, we drilled two holes, 18 mm apart, in the glass slide. Next, we cut two cover slips №2 (220 mm thick) in half and glued pairs of these pieces together using photocurable glue between two thicker pieces. We glued these double thick pieces of glass on the glass slide leaving a 1 mm gap between them (the holes should, of course, point to the gap). Next, we attached the top glass slide (without holes) using glue to create the channel. We removed any leftover glue using nitric acid and sodium hydroxide. We then blocked the ends of the channel using epoxy glue, and we made two rectangular shaped PDMS pieces with holes for tubing in them. Finally, we plasma-bonded these PDMS
	 22	

pieces on top of the holes to make connections for the channel inlet and outlet tubing. We noted that air bubbles could get stuck in the channel when this chip was used at low flow rates (< 1000 µL/h), and thus to prevent bubbles from getting stuck, we stuffed the inlet with crushed glass capillaries—the broken glass pieces created a hydrophilic net that prevented the continuity of a bubble from PDMS part of the inlet to the channel.
Bistability experiments. Four syringes were used: Syringe 1 (connected to CSTR), AlaSEt•TFA (92 mM) in water; Syringe 2 (connected to CSTR), CSSC•2HCl (184 mM) in 1.127 M phosphate buffer pH 7.5. Syringe 3 (connected to CSTR), maleimide (19.3, 38.7, or 58 mM) in water: Syringe 4 (connected to Detection chip), 16.6 mM Ellman’s reagent in mixture of MeOH/(NaH2PO4 208 mM) 2/3. Note: Ellman’s reagent should be dissolved in MeOH before adding the solution of NaH2PO4. We used the following flow rates (µL/min) to construct hysteresis plots. Syringes 1, 2: (i) 1.7626; (ii) 3.525; (iii) 7.05; (iv) 14.1; (v) 21.15; (vi) 28.2; (vii) 35.25; (viii) 42.3; (ix) 35.25; (x) 28.2; (xi) 21.15; (xii) 14.1; (xiii) 7.05; (xiv) 3.525; (xv) 1.7626. Syringe 3: (i) 0.225; (ii) 0.45; (iii) 0.9; (iv) 1.8; (v) 2.7; (vi) 3.6; (vii) 4.5; (viii) 5.4; (ix) 4.5; (x) 3.6; (xi) 2.7; (xii) 1.8; (xiii) 0.9; (xiv) 0.45; (xv) 0.225. Syringe 4: (i) 11.25; (ii) 22.5; (iii) 45; (iv) 90; (v) 135; (vi) 180; (vii) 225; (viii) 270; (ix) 225; (x) 180; (xi) 135; (xii) 90; (xiii) 45; (xiv) 22.5; (xv) 11.25.
Oscillatory experiments. Four syringes were used. Syringe 1 (connected to the CSTR), AlaSEt•TFA (140 mM) in water. To consider the volume of AlaSEt•TFA, 104 mg of AlaSEt•TFA were dissolved in 2.90 mL of Milli-Q water. Syringe 2 (connected to the CSTR), CSSC•2HCl (276 mM) in 3 M phosphate buffer, pH 8 (K2HPO4 and KH2PO4 have to be used to get a 3 M solution). To consider the volume of CSSC•2HCl, 193.2 mg of CSSC•2HCl was dissolved in 2.85 mL of the buffer. Syringe 3 (connected to the CSTR), Maleimide (31 mM) and Acrylamide (962 mM) in water. To consider the volume of acrylamide, 9 mg of maleimide and 205 mg of acrylamide were dissolved in 2.85 mL of
	 23	

Milli-Q water: Syringe 4 (connected to the mixing chip), 16.6 mM Ellman’s reagent in a mixture of MeOH/(NaH2PO4 208 mM) 2/3. Note: Ellman’s reagent should be dissolved in MeOH before adding the solution of NaH2PO4. We used the following flow rates (SV (s-1), Syringes 1, 2, 3 flow rate (µL/min), Syringe 4 flow rate (µL/min)): 1.2.10-4, 36, 324; 1.4.10-4, 42, 378; 1.6.10-4, 48, 432; 2.0.10-4, 60, 540; 2.5.10-4, 75, 675; 2.7.10-4, 81, 729; 2.8.10-4, 84, 756.
Numerical Modeling Numerical modeling was done by COPASI35 and MATLAB programs. MATLAB numerical models for the oscillatory and bistable networks are described in what follows. The architectures of our bistable and oscillatory networks are outlined in Extended Data Fig. 9a. Under flow conditions, the evolution of the species in the bistable network over time is governed by the following set of ordinary differential equations 1a-1h:
(1a)
(1b)
(1c)
(1d)
	 24	

(1e) (1f) (1g) (1h) Under flow conditions, the evolution of the species in the oscillatory network over time is governed by the following set of ordinary differential equations 2a-2i: (2a) (2b)
(2c)
(2d)
(2e)
	 25	

(2f)
(2g)
(2h)
(2i) We numerically solved these two systems of equations (equations 1 and 2) for the bistable network and oscillatory network for a variety of flow rates, FvV, using a custom code that relies on the MATLAB ode45 algorithm. Selected results of the numerical simulations are shown in Extended Data Fig5a-c, Fig.6a. Our MATLAB scripts are available as supplementary files. Linear stability analysis of the three variable model of the oscillator The full oscillator model can be reduced to a three variable model by applying three approximations: (i) the positive feedback loop is described by simple quadratic autocatalysis with rate constant k1; (ii) the negative feedback loop (reaction with acrylamide) is described by a first order reaction with rate constant k3; and (iii) we neglect all of the end products. The system of reactions will be described by:
	 26	

Where A corresponds to [EtSH] + [5] + [CSH]; I to [Maleimide]; S to [AlaSEt]; k1-4 are rate constants; and k0 is space velocity; stands for washout or for the formation of inactive products. The system of kinetic equations takes the following form:
(3a) (3b) (3c) To represent the experimental system, we used the following values: S0 = 0.05 M, I0 = 0.01 M, k1 = 0.25 s-1M-1, k2 = 300 s-1M-1, k3 = 0.0035 s-1, k4 = 7.10-5 s-1. Using these parameters, our numerical model demonstrates sustained oscillations with k0 from 1.5.10-4 to 3.10-4 s-1; the region of the sustained oscillations is similar to the region where sustained oscillations were observed in the experiment. To describe the experimentally observed dependence of the stability of oscillations on space velocity, we analyzed the stability of the steady states derived from the following considerations. The steady-states are defined by:
	 27	

equal to:

(4)

(5)

And the dependence on k0 for the given rate constants, and starting concentrations will be:
(6)
The dependences of the roots of this cubic equation from k0 for the region of interest (1.10-4 < k0 < 4.10-4) are shown on the graph in Extended Data Fig. 8a, b. An analysis of the discriminant of the cubic equation written above shows that solutions of this equation experience a Fold bifurcation at k0 = 3.08266724620437.10-4 s-1 where two positive solutions collide and disappear. To analyze the stability of the steady-states we wrote a Jacobian matrix (J) as follows:
(7)

We used a Mathematica script (supplementary file, x ≡ k0; h0 ≡ I0) to compute the stability
of the steady states. For k0 = 3.08266724620437.10-4 s-1 and [A]ss = 0.0000263986, eigenvalues defined by equation det[J] = 0 are λ1 = -0.11456; λ2 = -0.000247982; λ3 = - 1.59534.10-10 , thus λ3 approaches 0 as expected for a Fold bifurcation. For k0 higher than 3.08266724620437.10-4 s-1, the lowest steady state is the only stable steady state. For example, for k0 = 4.10-4 s-1 :
(i) [A]ss = 3.99596.10-6 { λ1 = -0.745063; λ2 = -0.000453008; λ3 = -0.000407547} – all eigenvalues are negative (stable state)
	 28	

(ii) [A]ss = 0.000242641 { λ1 = -0.088553; λ2 = -0.000319249; λ3 = 0.00427642} – λ3 > 0 (unstable state)
(iii) [A]ss = 0.00197459 { λ1 = -0.594797; λ2 = 0.000159607 - 0.00128566i; λ3 = 0.000159607 + 0.00128566i } – the real part of λ2 and λ3 is higher than zero (unstable state)
For k0 lower than 3.08266724620437.10-4 s-1, there is only one physically meaningful (positive and real) steady-state. This steady-state undergoes a transition from an unstable focus (eigenvalues are complex conjugates with a positive real part) to a stable focus (eigenvalues are complex conjugates with a negative real part) through Andronov-Hopf bifurcation (referred to as a Hopf bifurcation in the text). At the Hopf bifurcation, the real part of the eigenvalues should go through zero. Computations show that this is indeed the case (see also Extended Data Fig.8a):
(i) for k0 = 2.10-4 s-1 ; [A]ss = 0.00122403{ λ1 = -0.369038; λ2 = 0.0000302409 -0.00112125i; λ3 = 0.0000302409 +0.00112125i } - real part of λ2 and λ3 is positive, unstable focus resulting in a stable orbit and sustained oscillations.
(ii) for k0 = 1.7680658.10-4 s-1 ; [A]ss = 0.00111479{ λ1 = -0.336194; λ2 = 3.55662.10-13 0.00108786i; λ3 = 3.55662.10-13 + 0.00108786i } – real part of λ2 and λ3 approaches zero (within the limits of the precisions of the computation). This point is a point of the Hopf bifurcation.
(iii) for k0 = 1.4.10-4 s-1 ; [A]ss = 0.00093034{ λ1 = -0.280744; λ2 = -0.00006549760.00102364i; λ3 = -0.0000654976+0.00102364i } - real part of λ2 and λ3 is negative, stable focus resulting in damped oscillations.
Supplementary discussions Bifurcations.
A bifurcation point represents a set of conditions where a smooth change in a control parameter causes a sudden shift in system behavior. For example, a system might transition
	 29	

from having multiple steady states to having just one, or vice versa (Extended Data Fig. 9b). In our system, as space velocity is increased to extremely high values, we observe such a transition. The mathematical basis of this transition can be explained by the number of real solutions to a cubic equation (X3 – X + C = 0). These solutions are a function of control parameter C; as C is increased, two out of three real solutions disappear, indicating a bifurcation point.
Oscillatory systems can undergo bifurcations in which they transition from having a single stable steady state to having one unstable steady state with sustained oscillations. This transition, which is called a Poincare-Andronov-Hopf bifurcation, is apparent in our system. Hysteresis.
Hysteresis describes the failure of a system to return to its starting point when an external condition—or control parameter—is changed and then changed back (Extended Data Fig. 9b). The best-known example of hysteresis is the behavior of ferromagnetic materials in the magnetic field.
Bistable dynamic systems have two stable states (Extended Data Fig. 9b). Without a strong perturbation, the system does not transition from one state to the other. Consequently, when the control parameter is increased from zero, the system will remain at a first set of stable states (the lower solid line in Extended Data Fig. 9b) until a bifurcation point is reached; when the control parameter is again lowered, the system will remain at a second set of stable states (the upper line in Extended Data Fig. 9b) until a bifurcation point is reached. As a result, a hysteresis loop appears. Detailed chemical explanation.
Thioesters are carboxylic acid derivatives that can react with a variety of nucleophiles through nucleophilic substitution at their sp2 carbon. Such a reaction consists of two steps: (i) addition of the nucleophile; (ii) loss of the leaving group (e.g., thiolate). Our system contains
	 30	

AlaSEt as the thioester and three nucleophiles (hydroxyl, amines, and thiolates) with which it can react to form alanine (carboxylic acid), amides, and new thioesters (reaction with thiolate is also called thiol-thioester exchange). At pH 7.5-8, reactions with hydroxyls and amines are slow, and reactions with thiolate are fast,36 but reversible (the product is another thioester). Kent ligation relies on this difference in reaction rates. Thiol-thioester exchange goes quickly and is followed by a quick, entropically favorable intramolecular reaction with amine.
In our autocatalytic pathway, cysteamine (CSH) reacts with AlaSEt through Kent ligation. This ligation is followed by a fast thiolate-disulfide reaction, making the complete autocatalytic process fast,37 relative to other reactions in our system.
It is important to notice that two independent pathways result in alanine amides in our system: (i) the direct reaction of amines from CSSC, a noncatalytic reaction, and (ii) an autocatalytic pathway through CSH. The ratio between (i) and (ii) is 1/100 or more for our system. This ratio is responsible for enabling fine control over the autocatalytic process.
If we run an autocatalytic reaction in a CSTR, two new processes are added: (i) the continuous addition of thioester and disulfide through the inlet, and (ii) the removal of thiols and all other products through an outlet (removal of thiols slows amplification). The final concentration of thiols represents a balance between these two processes and the aforementioned reactions.
Maleimide is a conjugated carbonyl derivative, with which nucleophiles can readily react through conjugate addition (i.e., addition to the double bond next to the carbonyl group). Maleimide reacts with thiols much faster than with thioesters and disulfides involved in the autocatalytic process and, thus, gives rise to two steady states in a CSTR: (i) In the first state, which has very low concentrations of thiols, maleimide is supplied by flow faster than it is removed by reaction with thiols. (ii) In the second state, which has very high concentra-
	 31	

tions of thiols produced through the autocatalytic pathway, maleimide is consumed by thiols faster than it is supplied by flow.
The addition of acrylamide leads to oscillations. Acrylamide, similar to maleimide, is a conjugated carbonyl compound which reacts with the same nucleophiles as maleimide. Its reactivity, however, is several orders of magnitude lower than that of maleimide; as a result, acrylamide is not able to suppress the autocatalytic pathway, and is only able to remove thiols when the thioester is depleted. The mechanism by which the slow removal of an activator (i.e., the slow removal of thiols by acrylamide) converts a bistable system into an oscillatory system is described in detail by Epstein, I. R. & Pojman, J. A. An introduction to nonlinear chemical dynamics : Oscillations, waves, patterns, and chaos. (Oxford University Press, 1998), p 62 – 78. Briefly, acrylamide creates a mechanism of negative feedback—as the concentration of thiols is increased, their rate of acrylamide-mediated removal increases— which destabilize the high thiol steady state leaving only an unstable state at low flow rates. The system is forced to oscillate around this unstable state.
	
References
31.	 in	Handbook	of	Chemistry	and	Physics			(ed	W.	M.	Haynes)	(CRC,	2015).	 32.	 in	Chemical	handbook	Vol.	3		(ed	B.	P.	Nikolsky)	(Chimia,	Leningrad,	1971).	 33.	 Benesch,	R.	E.	&	Benesch,	R.	The	acid	strength	of	the	-SH	group	in	cysteine	and	
related	compounds.	J.	Am.	Chem.	Soc.	77,	5877-5881	(1955).	 34.	 Whitesides,	G.	M.,	Lilburn,	J.	E.	&	Szajewski,	R.	P.	Rates	of	thiol-disulfide	
interchange	reactions	between	mono-	and	dithiols	and	ellmans	reagent.	J.	Org.	 Chem.	42,	332-338	(1977).	 35.	 Hoops,	S.	et	al.	COPASI-	a	complex	pathway	simulator.	Bioinformatics	22,	30673074	(2006).	 36.	 Bracher,	P.	J.,	Snyder,	P.	W.,	Bohall,	B.	R.	&	Whitesides,	G.	M.	The	relative	rates	of	 thiol-thioester	exchange	and	hydrolysis	for	alkyl	and	aryl	thioalkanoates	in	 water.	Origins	Life	Evol.	B.	41,	399-412	(2011).	 37.	 Steinfeld,	J.	I.,	Francisco,	J.	S.	&	Hase,	W.	L.	Chemical	kinetics	and	dynamics.	2nd	 edn,		(Prentice	Hall,	1999).	
	 32	

Extended Data Figure 1. 1H NMR kinetic experiments showing the mechanism of autocatalysis in AlaSEt with CSSC reaction: (a) Scheme for the reaction between AlaSEt and CSH (b) Scheme for the reaction between AlaSEt and hexamethylendiamine (c) Kinetic plot for the between AlaSEt and CSH . Reaction conditions: D2O, pH 7.5, 500-mM phosphate buffer, 25 oC, AlaSEt (58 mM), CSH (88 mM); (d) Kinetic plot for the reaction between AlaSEt and hexamethylendiamine . Reaction conditions: D2O, pH 7.5, 500-mM phosphate buffer, 25 oC, AlaSEt (58 mM), hexamethylendiamine (88 mM); (e) Kinetic plot for “autocatalytic thiol network”. Reaction conditions: 500-mM phosphate buffer pH 7.5, 25 oC, D2O. Points represent experimental data, and solid lines represent the predictions of numerical simulations of the system. (We assume the following: all disulfide-thiol interchange reactions proceed with the same rate constant kSS = 0.65 s-1M-1, Kent ligation proceeds with rate constant kL = 0.41 s-1M-1, and hydrolysis proceeds with rate constant kW = 7.00.10-6 s-1). The poor fit of the experiment with [CSSC] = 23 mM is a result of a simplifying assumption that all thiolate-disulfide interchanges in the model occur with equal rate constants in the forward and reverse directions. (f) Time profiles for selected individual compounds in the “autocatalytic thiol network”. Reaction conditions: 500-mM phosphate buffer pH 7.5, 25 oC, D2O AlaSEt (46 mM); CSSC (46 mM).
	 33	

Extended Data Figure 2. Schematic representation of the thiol-disulfide interchange reactions, and thiol-thioester exchange reactions that we omitted from Fig. 1b of the main text of the paper.
	 34	

Extended Data Figure 3. Batch kinetic experiments. a-c. Kinetics plots for the reaction of AlaSEt with CSSC and maleimide. Reaction conditions: 500-mM phosphate buffer pH 7.5, 25 oC, D2O; AlaSEt (46 mM); CSSC (46 mM); maleimide: (a) 1.16 mM; (b) 2.31 mM; (c) 3.47 mM. d-e. Kinetics plots for batch reaction between AlaSEt, CSSC, maleimide and acrylamide. Reaction conditions: 1 M phosphate buffer pH 8, 25 oC, H2O; AlaSEt (46 mM); CSSC (92 mM); maleimide (10 mM); acrylamide: (d) 160 mM; (e) 320 mM.
	 35	

Extended Data Figure 4. Calibration curves for the determination of the total concentration of free thiols. (a) Calibration curve for the batch experiments [SH] (mM) = A*28.71; (b) Calibration curve for the detection system where microfluidic chip was coupled to the fiber optic spectrophotometer [SH] (mM) = (A - 0.01496)*21.32; (c) Trace of the absorbance from the calibration experiment for the detection system where glass flow cell was coupled to the Cary 60 UVVis spectrophotometer [SH] (mM) = (A - 0.03848)*7.16 (d) Calibration curve that was obtained from the data shown on c. ME stands for β-mercaptoethanol.
	 36	

Extended Data Figure 5. Bistability in the maleimide delayed autocatalytic thiol network. a-c. Modeling of the hysteresis in the maleimide delayed autocatalytic thiol network in CSTR with following parameters kSS = 0.65 s-1M-1; kL = 0.41 s-1M-1; kW = 9.26.10-6 s-1; kMal = 150 s-1M-1; [CSSC]0 = 100 mM; [AlaSEt]0 = 50 mM; [Maleimide]0 = (a) 1.16 mM; (b) 2.31 mM; (c) 3.47 mM. d-f. Hysteresis curves for bistability experiments with different maleimide concentrations. Reaction conditions: 0.5 M phosphate buffer pH 7.5, 25 oC; AlaSEt (47 mM); CSSC (92 mM); Maleimide: (d) 1.16 mM; (e) 2.31 mM; (f) 3.47 mM.
	 37	

Extended Data Figure 6. Identification of space velocities resulting in sustained oscillations. (a) Numerical simulations predicting sustained oscillations based on rate constants kSS = 0.444 s-1M-1; kL = 0.46 s-1M-1; kW = 6.64.10-6 s-1; kMal = 150 s-1M-1; kAAm = 0.014 s-1M-1 at: [CSSC]0 = 100 mM; [AlaSEt]0 = 50 mM; [Maleimide]0 = 10 mM; [Acrylamide]0 = 290 mM; SV = 0.0002 s-1 (b) Adjustment of the flow rate in the oscillatory experiments. Reaction conditions: 1 M phosphate buffer pH 8; AlaSEt (47 mM); CSSC (92 mM); Maleimide (10.3 mM); Acrylamide (321 mM).
	 38	

Extended Data Figure 7. Summary of all oscillatory experiments. Reaction conditions: 1 M phosphate buffer pH 8; AlaSEt (47 mM); CSSC. (92 mM); Maleimide (10.3 mM); Acrylamide (321 mM). The qualitatively different behavior of the system at SV = 2.8·10-4 s-1 is an indication of the bifurcation point. The apparent thiol concentration in the range 0.5 – 1 mM in experiments 2 and 3 of panel at SV = 2.8·10-4 s-1 is the result of a background reaction of AlaSEt with Ellman’s reagent that occurs with a rate constant that is about 4.104 times slower than the reported apparent second-order rate constant of the reaction between Ellman’s reagent and βmercaptoethanol at pH 7.36
	 39	

Extended Data Figure 8. Analysis of for the three variable model. (a) Results of the linear stability analysis of for the three variable model. The graphs show steady states, which were calculated from the model (S0 = 0.05 M, I0 = 0.01 M, k1 = 0.25 s-1M-1, k2 = 300 s-1M-1, k3 = 0.0035 s-1, k4 = 7.10-5 s-1), and their stability. The model is given by following differential questions: In this model, A corresponds to [RSH]; I to [Maleimide]; S to [AlaSEt]; k1-4 to rate constants; and k0 to space velocity. The right graph shows detailed analyses of the region selected by the dashed square. Details of the linear stability analysis can be found in the supplementary materials. (b) A phase plot computed from the three variable model. (c) Experimental hysteresis plot proving transition from oscillations to bistability with lowering of the feeding concentration of maleimide. Reaction conditions: 1 M phosphate buffer pH 8, 25 oC; AlaSEt (47 mM); CSSC (92 mM); Maleimide (4 mM); Acrylamide (321 mM). Dashed line shows correspondence between the phase plot and experimental data.
	 40	

Extended Data Figure 9. Details of the full numerical model and an illustration of bistability and hysteresis. (a) Schematic representation of the reactions that were considered in the full models of the bistable network and of the oscillatory network. (b) An example of a fold bifurcation. As the control parameter is increased, the system transitions from having two stable states (A and A’) to having just one (B). The transition from A to B is often called a critical transition. After the transition, lowering the control parameter will not return the system to A. This phenomenon is called hysteresis. When the control parameter is lowered sufficiently, the system will again transition from having two stable states (C and C’) to having just one (D). 	
	 41