Continuously tunable multistability in DNA replication networks

Domain: neuroscience
Published: 2026-02-05
Editorial outcome: major_revision

Abstract
Multistable systems, ideally, could stabilize at any desired, switchable state within a continuous spectrum. However, conventional systems, constrained by signal-mediated mutual activation or inhibition, are limited to a finite set of discrete steady states. Here, we propose a rational framework for achieving continuously tunable multistability through reversible displacement reaction-mediated competition between positive autoregulatory DNA polymerization/nicking modules. This framework harnesses the chemical energy of dNTP hydrolysis to suppress spontaneous interconversion between modules for stabilizing at any target state along a continuous compositional gradient. With unparalleled tunability, the framework enables continuous, orthogonal state transitions and concentration-adaptive molecular memory in response to transient stimuli. Moreover, the single-stranded DNAs generated by polymerization/nicking reactions can be customized with predefined structures and functions, enabling continuously multistable controlover downstream processes, e.g., biocatalysis and RNA transcription, while maintaining multistability. This framework establishes a versatile and robust platform for developing chemical and material systems with continuously tunable multistability.

Introduction
In cells, highly sophisticated molecular reaction networks1–4 orchestrate various biochemical activities to maintain cellular function and homeostasis. These networks operate via different mechanisms, including feedback5, feedforward6, intercommunication7, adaptation8, oscillation9, and multistability10, to achieve precise modulation of biological functions. Inspired by the complexity and fine-tuned modulation observed in natural systems, alternative strategies11–13 have been developed for bottom-up assembly of synthetic networks in systems chemistry14,15 and synthetic biology16,17 over the past decade. These efforts focus on replicating the properties of living systems18–20, such as environmental adaptation21–23 and complex information processing24–27. For instance, synthetic networks have been designed to govern temporal28–31, spatial32,33 or multistage34 self-assembly processes for emerging functions28,35, harmonize transitions between material steady states36–38, generate chemical or material spatial patterns through coupled reactions and diffusion39,40, or orchestrate spatiotemporal signaling and developmental patterning processes41–43. Multistability represents a fundamental form of functional self-organization, permitting systems to operate within multiple distinct regimes under identical external conditions. This phenomenon, characterized by inherent nonlinearity, has attracted growing interest across diverse scientific disciplines, including biology10,44–51, chemistry52–54, physics46, and material sciences55,56. In chemical and biological contexts, multistable molecular networks are crucial for maintaining molecularly distinct steady states, each associated with a specific functional role10,44–46. Additionally, these networks exhibit molecular memory47–51 whereby transient external signals induce sustained state changes. Such multistable behaviors are essential for dynamic systems performing diverse tasks without altering their structural setup, as seen in the processes of cell division49, apoptosis50, and differentiation47,51.

Enlightened by natural multistability, synthetic multistable networks have been constructed using small molecules37,52, peptides57,58, genelets59–62, PEN tools63,64, or CRISPRi systems65 as building blocks. These synthetic systems are paramount for advancing our understanding of multistability and establishing a foundation for controlling multistable processes. Such control has significant implications for therapeutic interventions and the creation of adaptive materials and synthetic systems with life-like properties. To further advance this field, three critical objectives must be addressed: (i) Expanding the number of accessible steady states to enhance the system’s capacity for complex decision-making and functionality; (ii) Enabling adaptive switching between steady states in response to dynamic environments; (iii) Regulating downstream processes through multistable networks. However, both natural and synthetic multistable networks to date predominantly rely on signal-mediated mutual activation or inhibition between network modules to establish bistability45. To create more intricate multistable systems, a common strategy involves the combination of multiple discrete bistable elements to form an overall multistable configuration. Therefore, a multistable network with N modules (A1 to AN) requires N×(N−1) inhibitory reactions to support mutual inhibition between the modules, as illustrated in Fig. 1a. For instance, a bistable network with two modules (A1 and A2) needs two inhibitory reactions (A1 inhibiting A2 and A2 inhibiting A1). For a tristable network with three modules (A1, A2, and A3), six inhibitory reactions where each module inhibits the other two (A1 inhibiting A2 and A3, A2 inhibiting A1 and A3, A3 inhibiting A1 and A2) are required. As the number of modules increases, the complexity of designing such networks grows quadratically, making it increasingly challenging to balance signaling activities and avoid unintended interactions. More importantly, these networks stabilize in discrete steady states, limiting their adaptability and flexibility in response to dynamic environmental cues. For example, bistable networks generally toggle between two distinct states (e.g., “ON” and “OFF”), precluding the stepless or fine-tuned control of functional states. These limitations raise a fundamental question: whether it is possible to engineer a multistable system that can stabilize in any desired, switchable state rather than in a limited set of discrete steady states. Addressing this challenge would represent a significant advancement in controlling multistability, enabling unprecedented control over dynamic behaviors and broadening the applicability of synthetic multistable networks across diverse domains66,67.

Although largely unexplored in molecular programming, continuous multistability has precedents in other fields. For example, in neuroscience, continuous attractor networks are well-established models for processing analog information68–70. The oculomotor integrator, a neural circuit that maintains arbitrary eye positions in response to transient inputs, provides a canonical example of persistent analog memory. Similarly, in mechanical metamaterials, arrays of bistable elements can adopt a continuum of structural configurations in response to global inputs, illustrating how discrete components can be organized to exhibit graded behavior71.

In this study, we present a rational design for constructing scalable and controllable continuously tunable multistability networks (CTMNs) through reversible displacement reaction-mediated competition among positive autoregulatory DNA polymerization/nicking modules, as illustrated in Fig. 1b. By dissipating free energy from a dNTP reservoir, the system operates under far-from-equilibrium conditions where competition among these modules suppresses their spontaneous transitions, thereby stabilizing the system at any desired compositional state, a phenomenon we term continuously tunable multistability. We validate this continuously tunable multistability by experimental data and mathematical modeling. The tunable interconversion between modules mediated by reversible displacement reactions enables continuous, orthogonal state transitions and concentration-adaptive molecular memory in response to transient stimuli. Moreover, reversible displacement reactions simplify the inherent complexity of conventional signaling pathways, offering a significant advantage in modular scalability. This approach allows for the construction of larger-scaled, more complex CTMNs by seamlessly incorporating additional DNA polymerization/nicking modules. A unique feature of these networks is their ability to achieve precise, state-specific DNA replication via polymerization/nicking reactions, enabling the on-demand synthesis of single-stranded DNAs (ssDNAs) that can encode diverse structural and functional information. These ssDNAs further extend the system’s versatility by regulating downstream processes, such as DNAzyme-catalyzed reactions72 and RNA transcription73, in a stepless and programmable manner, without compromising the system’s multistability.

Results
Design of a continuously tunable multistability network
We developed CTMNs based on reversible displacement reaction-mediated far-from-equilibrium competition between positive autoregulatory DNA polymerization/nicking modules, such as CTMN-1, illustrated in Fig. 2 and Supplementary Fig. 1. The network consists of four precisely engineered DNA duplexes: AA’, BB’, BA’, and AB’. Each duplex incorporates an exchangeable double-stranded domain (ZZ’), facilitating a reversible displacement reaction between competing modules AA’+BB’ and BA’+AB’74,75. This reversible displacement reaction offers two critical advantages over traditional designs. First, in contrast to traditional systems producing binary outcomes—one module is fully activated while the other is entirely suppressed, a reversible displacement reaction allows precise, proportional control over modules’ fractions by modulating their transitions, enabling the fine-tuning of network states. Second, employing a reversible displacement reaction to mediate competition between network modules could simplify the architectural design relative to traditional systems, which often rely on complex signaling pathways to establish mutual inhibition or activation between modules. Therefore, this competitive architecture could reduce assembly complexity and improve scalability and robustness (vide infra).

To engineer a far-from-equilibrium competition between these two modules, we coupled this reversible displacement reaction with a dNTP-fueled DNA polymerization/nicking reaction system, including DNA templates, DNA polymerase I, Klenow fragment (3′→5′ exo-), and nicking endonuclease Nt.BbvCI, and dNTPs. Specifically, A’ and B’ serve as DNA templates, with template domains P’1 and P’2 incorporated at their 5′ ends, respectively. Strand A interacts with A’ or B’ to form duplex AA’ or AB’, which initiate the dNTP-driven polymerization on the template domain ((DNA)n + dNTP → (DNA)n+1 + PPi), yielding AA’P1 or AB’P2. In contrast, strand B was modified with a 7-base tether at its 3′ end to inhibit the polymerization reaction of duplexes BA’ and BB’. In addition, polymerized products AA’P1 and AB’P2 were pre-engineered with a nicking site for Nt.BbvCI in their elongated domains. These nicking reactions generate AA’ and AB’, restoring the exchangeable domain ZZ’ and enabling these structures to re-enter both the strand displacement and polymerization/nicking cycles. Thus, this design establishes autonomous polymerization/nicking cycles on AA’ and AB’ driven by dNTP hydrolysis. Under appropriate conditions, the polymerization/nicking reactions could swiftly convert the majority of AA’ and AB’ into their polymerized forms, AA’P1 and AB’P2. The formation of AA’P1 isolates AA’ from the reversible displacement reaction, inhibiting the forward transition from AA’+BB’ to BA’+AB’ while maintaining the reverse conversion of BA’+AB’ back to AA’+BB’. This establishes a positive autoregulatory loop on AA’+BB’, which upregulates AA’+BB’ and downregulates BA’+AB’. To elucidate the thermodynamic and kinetic principles governing this regulatory architecture, a detailed mechanistic framework along with corresponding NUPACK-predicted results is provided in Supplementary Section 2 and Supplementary Figs. 2 and 3. Similarly, the generation of AB’P2 isolates AB’ from the reversible displacement reaction, inhibiting the reverse transition from BA’+AB’ to AA’+BB’ while retaining the forward conversion of AA’+BB’ into BA’+AB’. This yields an opposite, positive autoregulatory loop on BA’+AB’ that upregulates BA’+AB’ and downregulates AA’+BB’. These competitive dynamics result in a simultaneous and rapid inhibition of the forward and reverse transitions between AA’+BB’ and BA’+AB’, effectively “freezing” the reversible displacement reaction in a far-from-equilibrium manner. This mechanistic interplay between dNTP hydrolysis and the inhibition of strand exchange exemplifies how continuous energy consumption actively suppresses spontaneous, thermodynamically favored interconversion among network modules. Once the polymerization/nicking reactions reach a balanced state where AA’ and AB’ in the reversible displacement reaction are minimized, the reversible displacement reaction enters into an equilibrated state, allowing the network to reach a steady state. Therefore, the interplay between the polymerization/nicking reactions and reversible displacement reaction establishes a fine-tuned balance under any compositional state, endowing the network with continuously tunable multistability.

This design adheres to the fuel-driven, far-from-equilibrium principles observed in other multistable biochemical systems, such as genelet61,62 and PEN toolkits63,64, where sustained nucleotide hydrolysis fuels multistable behaviors. In these systems, RNA or DNA strands produced via polymerization serve as regulatory species, mediating mutual activation or inhibition to sculpt system dynamics. The concentrations of these species are tightly regulated through coordinated synthesis and degradation, allowing precise control of network behavior. Meanwhile, inert degradation byproducts gradually accumulate but remain functionally decoupled from the core regulatory circuitry. Our CTMN architecture is mechanistically distinct. In this case, the regulatory species are not ssDNA byproducts of the polymerization/nicking cycle (Px), but rather the double-stranded polymerized complexes AA’P1 and AB’P2. These complexes directly mediate inter-module competition and define the network’s multistability. Their concentrations are regulated through a reversible polymerization/nicking cycle that continuously regenerates AA’ and AB’, thereby preventing uncontrolled accumulation and supporting dynamic stability. Although byproducts Px accumulate over time, they are decoupled from the network’s regulatory architecture and do not impact its steady states. This mirrors design strategies in other synthetic biochemical systems, where non-functional degradation fragments or useless monomers accrue without influencing system dynamics.

Experimental validation and kinetic modeling
To monitor the network’s dynamic behavior, we modified strand A with quencher BHQ2 and labeled strands A’ and B’ with fluorophores FAM and ROX, respectively (details in Supplementary Table 1). To ensure well-defined initial conditions, each duplex component of the network modules was pre-annealed prior to its introduction into the reaction mixture. The time-dependent fluorescent changes associated with the network were recorded, and the measured data were subsequently converted to concentrations of the corresponding DNA duplexes using calibration curves that correlate fluorescence intensities with DNA duplex concentrations (Supplementary Fig. 4).

First, we investigated the reversible displacement reaction between modules AA’+BB’ and BA’+AB’ in the absence of DNA polymerase, nicking endonuclease, and dNTPs, thus excluding polymerization/nicking reactions. As depicted in Fig. 3a–c and Supplementary Fig. 5, the concentrations of the initially dominant modules gradually decline while those of the competing modules increase, eventually converging to a single thermodynamically equilibrated state under two different initial compositions: one with 100% AA’+BB’ and the other with 100% BA’+AB’. These dynamics were further validated by electrophoretic PAGE gel analysis (Supplementary Section 3 and Supplementary Fig. 6). We formulated a kinetic model to simulate these dynamics computationally (Supplementary Section 4). The fitted computational curves represented by dashed lines in Fig. 3b, c generally match the experimental data. The rate constants derived from this model are presented in Supplementary Table 2. Computational analysis reveals that the reaction rates of the initially dominant modules progressively decrease while those of the competing modules increase (Fig. 3d and for the explanation of asymmetry kinetics, see Supplementary Section 2). Once the forward and reverse rates are approximately equal, the system reaches a thermodynamically equilibrated state.

Next, we studied the impact of a single positive autoregulatory loop on the network dynamics by selectively modifying specific network components. To configure a network with an autoregulatory loop of BA’+AB’, we modified strands A’ by removing template domain P’1 (renaming it A’) while keeping template domain P’2 on strand B’. This modification selectively allows the polymerization/nicking reactions associated with AB’ in the presence of DNA polymerase, nicking endonuclease, and dNTPs, thereby establishing an autoregulatory loop on BA’+AB’ (Fig. 3e). For an alternative network configuration, we modified strand B’ by removing its template domain P’2 (designating it B’) while preserving template domain P’1 on strand A’. Under the same conditions, this arrangement selectively permits the polymerization/nicking reactions associated with AA’, thus generating an autoregulatory loop on AA’+BB’ (Fig. 3j). These structural adjustments allow us to investigate how variations in autoregulatory loops affect network dynamics, particularly regarding multistability.

We investigated the dynamics of these network configurations under two distinct initial compositions: one with 100% AA’+BB’ and the other with 100% BA’+AB’. For the network containing the BA’+AB’ autoregulatory loop, when initiated with 100% BA’+AB’, the network sustains its BA’+AB’-dominated composition over 12 h. Even with a different initial composition of 100% AA’+BB’, this network quickly shifts to an identical BA’+AB’-dominated state (Fig. 3f, g, and Supplementary Fig. 7). These results indicate that the BA’+AB’ autoregulatory loop endow the network with a single steady-state output, regardless of its initial composition. Conversely, the network with the AA’+BB’ autoregulatory loop exhibits distinct dynamics. When initiated with 100% AA’+BB’, the network sustains the initial composition over 12 h (Fig. 3k and Supplementary Fig. 8a). Interestingly, when beginning with 100% BA’+AB’, the network transitions to an intermediate steady state with a composition of 65% AA’+BB’ and 35% BA’+AB’ (Fig. 3l and Supplementary Fig. 8b). These results demonstrate that the AA’+BB’ autoregulatory configuration allows the network to achieve multiple steady states dependent on its initial composition.

To gain further insights into these dynamics, we developed a kinetic model (Supplementary Section 4) for each autoregulatory network configuration. Simulated curves generated by these models (dashed lines in Fig. 3f, g, k, l) are generally consistent with experimental data, and the derived rate constants are listed in Supplementary Table 2. Computational analysis indicates that removing the respective template domain alters the rate constants governing the reversible displacement reaction in both autoregulatory network configurations (Supplementary Figs. 9 and 10), as compared to the unmodified system. As shown in Fig. 3h, i, m, n, when the initially dominant module is in accord with that in the autoregulatory loop (e.g., the BA’+AB’ autoregulatory network initiated with 100% BA’+AB’, Fig. 3g or the AA’+BB’ autoregulatory network initiated with 100% AA’+BB’, Fig. 3k), the polymerization rate is initially high but declines rapidly, followed by a concurrent increase in the nicking rate of the polymerized duplex, implying that the polymerization/nicking reactions quickly convert the polymerizable duplex into its polymerized form (Fig. 3i, n). The produced polymerized duplex effectively isolates the polymerizable form in the dominant module from the reversible displacement reaction, thereby suppressing the transition from the initially dominant module to its competing module, as indicated by a sharp decrease in the transition rate and minimal production of the competing module (Fig. 3h, m). As the polymerization/nicking reaction cycle progresses, the nicking rate increases further, eventually balancing with the decreasing polymerization rate. This process drives the reversible displacement reaction to equilibrium swiftly, establishing a steady state that mirrors the initial composition.

In contrast, when the network’s initial configuration differs from that of the autoregulatory loop (e.g., the BA’+AB’ autoregulatory network initiated with 100% AA’+BB’, Fig. 3f, or the AA’+BB’ autoregulatory network initiated with 100% BA’+AB’, Fig. 3l), the activation of autoregulatory polymerization/nicking reactions necessitates a preceding conversion of the initially dominant module into the autoregulatory (polymerizable) module via the reversible displacement reaction. This conversion step modulates the kinetics of polymerization/nicking reactions, regulating the dynamic equilibrium between the polymerizable species and their polymerized counterparts. Consequently, this modulation impacts the isolation of the polymerizable species from the reversible displacement reaction, thereby influencing the inhibition dynamics on the transition from the autoregulatory module back to the initially dominant module. For the BA’+AB’ autoregulatory network initiated with 100% AA’+BB’ (Fig. 3f), the transition from AA’+BB’ to BA’+AB’ proceeds rapidly, leading to a swift accumulation of polymerizable duplex AB’. This rapid transition leads to an initial surge in the polymerization rate of AB’, followed by its subsequent decrease. Concurrently, the nicking rate associated with AB’P2 rises. Eventually, the polymerization/nicking rates (the interconversion between AB’ and its polymerized form) reach a balance at 17 nM h−1, identical to the balanced rate observed in the BA’+AB’ autoregulatory network initiated with 100% BA’+AB’ (Fig. 3i). This result indicates that the steady-state interconversion dynamics of AB’ and AB’P2 are invariant to the initial composition, stabilizing the network in the same BA’+AB’-dominant steady state irrespective of initial conditions.

Conversely, in the AA’+BB’ autoregulatory network initiated with 100% BA’+AB’ (Fig. 3l), the transition from BA’+AB’ to AA’+BB’ proceeds more slowly (Fig. 3n). As polymerizable duplex AA’ gradually accumulates, its polymerization rate rises rapidly at first but slows over time, accompanied by a proportional increase in the nicking rate of AA’P1. In this case, the balance between polymerization and nicking rate is achieved at 37 nM h−1, which is lower than 48 nM h−1 observed in the AA’+BB’ autoregulatory network initiated with 100% AA’+BB’ (Fig. 3n). This reduced rate reflects a smaller total content of AA’ (including AA’ and its polymerized form) and a greater residual concentration of BA’+AB’ in the network, resulting in an intermediate steady state with a composition of 65% AA’+BB’ and 35% BA’+AB’, which is distinct from the AA’+BB’-dominant steady state achieved when the network is initiated with 100% AA’+BB’. Thus, when the network’s initial configuration diverges from that of the autoregulatory loop, the transition rate of the initially dominant module to the autoregulatory (polymerizable) module significantly influences the subsequent polymerization kinetics associated with the autoregulatory module. These changes, in turn, alter the balanced rates of polymerization and nicking reactions, ultimately leading to distinct steady-state compositions. Such kinetic interdependencies are critical for the emergence of multistability in the AA’+BB’ autoregulatory network.

To further support our analysis, we employed this kinetic model to predict the network’s multistable behaviors by varying the initial composition, which was then validated experimentally. As depicted in Supplementary Fig. 11, the AA’+BB’ autoregulatory network demonstrates distinct dynamics depending on its initial composition. For an initial composition of 30% AA’+BB’ and 70% BA’+AB’, which is lower than the balanced fraction (65% AA’+BB’ and 35% BA’+AB’) observed in the network initiated with 100% BA’+AB’, the polymerization/nicking reaction remains limited by the transition rate between modules, and the network ultimately stabilizes at the same intermediate steady state as the one achieved with 100% BA’+AB’. Conversely, when the initial fraction of AA’+BB’ increases to 70%, exceeding the balanced fraction observed in the network initiated with 100% BA’+AB’, the network evolves towards a new intermediate steady state with a balanced polymerization/nicking rate of 41 nM h−1. These results demonstrate the network’s capacity for multistability, allowing the emergence of distinct steady states depending on its initial composition.

Enlightened by the observed multistability within a narrow compositional range in the network governed by a single autoregulatory loop of AA’+BB’, we sought to extend the multistable region by introducing two competitive autoregulatory loops, AA’+BB’ and BA’+AB’, in CTMN-1. In the presence of 30 U mL−1 DNA polymerase, 20 U mL−1 nicking endonuclease, and 0.1 mM dNTPs, we examined the dynamic behaviors of CTMN-1 across five distinct initial compositions, varying AA’+BB’ from 100%, 75%, 50%, 25% to 0%, with corresponding variations in BA’+AB’ from 0, 25, 50, 75 to 100%. This setup enabled us to assess the network’s dynamic behavior across the entire range of initial compositions. Supplementary Fig. 12 shows the time-dependent concentration profiles of AA’, AB’, BA’, and BA’ derived from the corresponding fluorescence data (Supplementary Fig. 13) using appropriate calibration curves (Supplementary Fig. 4). Across all initial conditions, CTMN-1 exhibits a two-phased dynamic behavior. The initial phase (~0.5 h) is characterized by a short relaxation, during which polymerization/nicking reactions exert an increasing, suppressive influence on the reversible displacement reaction. This is followed by a second, prolonged phase corresponding to the operational steady state. In this regime, polymerization/nicking reactions have fully established kinetic dominance, effectively halting further spontaneous interconversion between modules. However, this regime is not entirely static—compositional changes occur at a markedly slower rate.

Under all initial conditions, CTMN-1 reliably settles into composition-specific steady states (S1-1 to S1-5), each exhibiting small compositional drift over a 12-h window (varying from 0.55 to −0.15% h−1; Supplementary Fig. 12f). Comparable drift phenomena have been reported in conventional multistable biochemical systems58,61–64. To explore the source of this drift, we first characterized its time dependence at 0.1 mM dNTPs. Under these conditions, the steady-state composition remains nearly constant for ~36 h, after which the drift becomes pronounced, reaching ~0.7% h−1 during the final 12-h window (Supplementary Section 5 and Supplementary Fig. 14a, d). The delayed increase in drift suggests a progressive fuel depletion as the dominant contributor rather than an intrinsic instability of the network. To further evaluate the role of fuel availability, we systematically investigated the influence of the dNTP concentration on steady-state stability (Supplementary Section 5). Using S1-1 as a representative example, we found that at lower dNTP concentrations (0–0.02 mM), multistability is compromised, and compositional drift increases markedly (Supplementary Figs. 15 and 16). In contrast, raising the dNTP concentration to 0.2 mM reduced its drift rate to 0.16% h−1 by 24 h (Supplementary Fig. 14b, e). Further increasing dNTP concentration beyond 0.2 mM does not yield additional improvement, indicating that fuel-dependent stabilization saturates near this concentration (Supplementary Fig. 14c, f). To quantify robustness across the entire multistability regime, we next monitored the drift rates of five representative steady states (S1-1 to S1-5) at the optimal dNTP concentration of 0.2 mM. Under this condition, these drift rates were substantially suppressed (Fig. 3o–s and Supplementary Fig. 17), with rates by 24 h varying from 0.16 to −0.05% h−1 (Fig. 3t), corresponding to the mean drift rate of 0.04% h−1. Moreover, CTMN-1 maintains high compositional stability over at least 48 h under optimal fuel conditions. For S1-1, the drift rate during the final 24-h window remains as low as 0.07% h−1. These findings collectively underscore that the observed multistability is not a passive thermodynamic property but arises from the network’s sustained, energy-consuming, far-from-equilibrium dynamics. The polymerization/nicking cycles continuously consumed NTPs to “lock” the system into one of multiple steady states near its initial composition. Control experiments lacking the energy-consuming polymerization/nicking reactions further confirm this interpretation: all systems converge to a single equilibrium state regardless of initial conditions (Fig. 3b, c), reinforcing the necessity of ongoing energy consumption for sustaining multistability—a hallmark of many natural biochemical networks.

To evaluate the robustness of CTMN-1 under other biochemical parameters, we also examined the operational ranges of DNA polymerase and nicking endonuclease concentrations at 0.1 mM dNTPs (Supplementary Section 5 and Supplementary Figs. 18–23). Across a broad span of enzyme stoichiometries, CTMN-1 maintained continuously tunable multistability, indicating that the network architecture is not sensitive to moderate changes in enzymatic balance. At a higher fuel dNTP concentration of 0.2 mM, shifting the system toward a polymerase-rich, nicking-limited condition does not lead to noticeable changes in the drift rate (Supplementary Figs. 24 and 25). Together, these analyses show that CTMN-1 maintains reliable multistability across a broad range of enzymatic and fuel conditions.

The multistable behavior was further corroborated by PAGE gel electrophoresis analysis (Supplementary Section 3 and Supplementary Fig. 26), providing orthogonal confirmation of the composition-specific steady states. Notably, although CTMN networks were often initialized with stoichiometrically balanced DNA duplexes, symmetry in initial composition is not a prerequisite for multistability. As demonstrated in Supplementary Fig. 27, CTMN-1 reliably maintains a well-defined steady state even under an asymmetric compositional condition, further illustrating the network’s flexibility.

To gain deeper insights into the network’s continuously tunable multistability, we developed a kinetic model (Supplementary Section 4) to simulate the experimental results. As shown in Fig. 3d, the simulated curves are generally consistent with experimental results, and the derived rate constants are detailed in Supplementary Table 2. As shown in Fig. 3v–x, computational analysis reveals that, regardless of initial compositions, the polymerization reactions of AA’ and AB’ proceed rapidly with their initially high polymerization rates declining swiftly. This is followed by a concurrent increase in the nicking rates of AA’P1 and AB’P2. These polymerization reactions rapidly convert AA’ and AB’ into their respective polymerized forms, efficiently isolating these species from the reversible displacement reaction and minimizing their forward and reverse rates. This competitive interaction between two autoregulatory modules leads to the “freezing” of the dynamic transition between them. Under these conditions, the polymerization/nicking reactions and the forward/reverse displacement reaction swiftly reach equilibrium, resulting in five distinct steady states corresponding to their respective initial compositions. The predictive power of the kinetic model was further validated by simulating CTMN-1’s multistable behaviors under two additional initial compositions, which were experimentally validated (S1-6 and S1-7 in Fig. 3u and Supplementary Fig. 28). These results demonstrate CTMN-1’s capacity to stabilize at any compositional configuration, yielding linearly distributed, continuously tunable steady states and emphasizing the reliability and predictive reliability of the kinetic model.

To systematically characterize the continuously tunable multistability observed in CTMN-1, we have performed Linear Stability Analysis (LSA) based on the fitted kinetic model to map the network’s full accessible steady-state space. Drawing inspiration from computational neuroscience69,70, where both linear and nonlinear network models have been proposed to explain continuous attractors and graded memory states, we note that ideal continuous attractors, comprising an infinite number of exact steady states, are theoretically achievable in linear systems. In contrast, biologically realistic nonlinear systems more commonly give rise to approximate continuous attractors. Although they are not mathematically exact, they are typically robust and functionally meaningful within a biological context. Guided by this conceptual framework, we adopted a criterion for steady-state identification where the net concentration flux of each molecular species is constrained to be less than 1×10−10 nM h−1 in LSA, consistent with established methodologies in nonlinear dynamical systems modeling. In this case, we identified thousands of steady states (for detailed LSA simulations, see Supplementary Section 4). From these, 216 representative states were uniformly sampled across the full accessible compositional range and are depicted in Fig. 3u (circle dots). The distribution densely populates a linear manifold, providing strong theoretical support for the emergence of linear, continuous multistability in CTMN-1. This organization reflects the emergence of an approximate continuum of attractors, reinforcing our interpretation of CTMN-1 as a system that exhibits biologically plausible, continuously tunable multistability. As a control, we conducted an identical LSA in an enzyme-free variant of the system, using the same net concentration flux applied to the 216 representative steady states. Unlike the multistable behavior in CTMN-1, the enzyme-free system consistently yields a single fixed-point equilibrium. This contrast highlights the essential role of the competing polymerization/nicking reactions in generating and maintaining the multistability observed in CTMN-1.

Orthogonal and continuous state transitions and concentration-adaptive network memory
The ability of multistable networks to switch their steady states in response to transient stimuli is fundamental for adaptive behavior and molecular memory. In conventional multistable systems, state transitions typically arise from the interference of signaling pathways between network modules, which switch on/off their activity levels. In contrast, state transitions within the CTMN are governed by the modulation of the reversible displacement reaction. Owing to the reversible dynamic behavior of the reversible displacement reaction, the CTMN enables orthogonal and continuous state transitions, in which each module’s activity can be regulated in a gradient-like, stepless fashion.

To fulfill this, we designed two orthogonal transient triggering modules to induce orthogonal state transitions of CTMN-1 (Fig. 4a, c). The transient triggering module consists of an RNA triggering strand, rTAA’ or rTAB’, alongside the enzyme RNase H. In CTMN-1, the single-stranded segments of AA’ and AB’ serve as distinct target domains complementary to their respective RNA triggers, rTAA’ and rTAB’. The state transitions of the network proceed as follows: The hybridization of rTAA’ with AA’ stabilizes it as AA’-rTAA’, which promotes a reconfiguration of the network. This reconfiguration upregulates module AA’+BB’ and downregulates module BA’+AB’. The mechanistic rationale governing this reconfiguration is detailed in Supplementary Section 2 and Supplementary Figs. 29 and 30. Following the reconfiguration, RNase H catalyzes the cleavage of rTAA’ in the complex AA’-rTAA’, restoring the original AA’ while maintaining the reconfigured steady state. Similarly, the application of rTAB’ to CTMN-1 stabilizes AB’ via their hybridization, forming AB’-rTAB’. This induces another orthogonal reconfiguration, upregulating module BA’+AB’ while downregulating module AA’+BB’. The cleavage of rTAB’ in the AB’-rTAB’ complex by RNase H restores AB’ while retaining the configured steady state.

To investigate the orthogonal state transitions of CTMN-1, we established two steady states, S1-1 and S1-5, as before. CTMN-1 in states S1-1 and S1-5 was subsequently treated with rTAB’ or rTAA’, respectively, followed by RNase H digestion. As presented in Fig. 4b and Supplementary Fig. 31a, applying rTAB’ to CTMN-1 in steady state S1-1 stabilizes AB’, facilitating the forward transition from AA’+BB’ to BA’+AB’. Consequently, CTMN-1 reconfigures into a new steady state, in which the compositional fraction of BA’+AB’ increases while that of AA’+BB’ decreases. Subsequently, the catalytic degradation of rTAB’ in complex AB’-rTAB’ by RNase H restores AB’, allowing the system to maintain the newly established steady state. Similarly, introducing another orthogonal trigger rTAA’ to the network in steady state S1-5 induces a reverse transition from BA’+AB’ to AA’+BB’ (Fig. 4d and Supplementary Fig. 31b), reconfiguring the network into another steady state where AA’+BB’ upregulates while BA’+AB’ downregulates. The degradation of rTAA’ by RNase H recovers AA’ without altering the reconfigured network state. These results demonstrate that the transient stabilization of AB’ by rTAB’ or AA’ by rTAA’ induces orthogonal and sustained reconfiguration of CTMN-1, encoding transient signals into distinct steady states. These reconfigurations are driven by significant shifts in reversible displacement reaction dynamics, the forward rate increases from 46.6 to 172.0 nM h−1 upon rTAB’ application, while the reverse reaction accelerates from 7.8 to 200.0 nM h−1 after applying rTAB’ (Supplementary Fig. 32). To quantitatively understand the mechanism underlying these transitions, we developed a comprehensive kinetic model (Supplementary Section 4) to simulate the time-dependent compositional changes of CTMN-1 in response to rTAA’ or rTAB’. As illustrated in Supplementary Fig. 33, the simulated curves are generally consistent with experimental observations, validating the model’s reliability. Rate constants derived from model fitting were summarized in Supplementary Table 2. In control experiments without DNA polymerase and nicking endonuclease, the system, regardless of the initial conditions, converges to a thermodynamically equilibrated state (Supplementary Fig. 34). Upon treatment with the respective RNA trigger, the system undergoes orthogonal dynamic reconfiguration. However, once the RNA trigger is degraded by RNase H, the system reverts to its original thermodynamically equilibrated state. These outcomes suggest that transient triggers cannot exert sustained reconfiguration in a monostable system.

Next, we sought to explore the potential for continuous state transitions in CTMN-1. As depicted in Fig. 4e and Supplementary Fig. 35, three successive applications of rTAB’ to CTMN-1 in state S1-1 induce three sequential state transitions along the forward pathway from AA’+BB’ to BA’+AB’. Similarly, triplicate applications of rTAA’ to CTMN-1 in state S1-5 result in three consecutive state transitions in the reverse direction, reverting from BA’+AB’ back to AA’+BB’. These transition patterns were also simulated by our established kinetic model, as detailed in Supplementary Section 4. Simulated curves (Supplementary Fig. 36) are generally consistent with experimental observations, supporting the robustness of the model. These results demonstrate the orthogonal and continuous nature of state transitions in CTMN-1. It is essential to note that, following the addition of the trigger, a marked asymmetry arises in the concentration change of the DNA duplexes within each module. Specifically, the DNA duplex that directly interacts with the trigger undergoes a more pronounced concentration change compared to its non-interacting counterpart. This disparity arises from the interplay between the intrinsic duplex formation–dissociation equilibrium (e.g., AA’ ⇌ A+A’, Supplementary Fig. 3a) and the selective, sequence-specific hybridization between the trigger and its corresponding single-stranded components (e.g., rTAA’+A+A’ → rTAA’AA’). Consequently, the concentration of the directly involved in trigger binding is perturbed more significantly (Supplementary Fig. 30). Importantly, this initial asymmetry is further amplified over time by a molecular memory effect intrinsic to CTMN-1. While RNA triggers are degraded by RNase H, their reintroduction in subsequent cycles enables recursive activation, which reinforces specific molecular states through cumulative effects. This cumulative effect progressively biases the system toward the accumulation of particular molecular species. We also explored the system’s behavior in response to multiple RNA triggers introduced either concurrently or sequentially. When two competing RNA triggers with different ratios were added simultaneously, the system consistently stabilizes at distinct steady states in a tunable, concentration ratio-dependent way (Supplementary Fig. 37). This outcome demonstrates that CTMN-1 enables tunable and deterministic control over state selection through input stoichiometry. Moreover, the system exhibits full reversibility: after reaching a stabilized state with rTAA’, the subsequent introduction of rTAB’ successfully reverts the system to its original state, implying the reversibility of CTMN-1’s multistable dynamics (Supplementary Fig. 38). Lastly, CTMN-1 can adapt to varying concentrations of transient stimuli to adjust its steady states. As shown in Fig. 4f and Supplementary Fig. 39, progressively increasing the concentration of rTAB’ applied to CTMN-1 in state S1-1 shifts the system to steady states with BA’+AB’ increasing and AA’+BB’ decreasing. This concentration-adaptive behavior suggests that CTMN-1 can encode the concentration information of transient stimuli into distinct long-lasting steady states, a property termed concentration-adaptive network memory that has not been observed in conventional multistable systems.

Scaling the continuously tunable multistability network
In conventional multistable networks, integrating additional network modules for scalability requires the establishment of complex regulatory interactions between the introduced modules and each of the existing ones. As the number of introduced modules increases, the signaling complexity increases quadratically, posing a significant challenge to the scalability of these systems. In contrast, our design, based on the reversible displacement reaction-mediated competitive interaction between network modules, provides a more simplified approach to scalability. It enables the modular expansion of the CTMN by directly integrating multiple pre-designed DNA duplexes, each containing a shared exchangeable double-stranded domain and a distinct template domain. These duplexes naturally compete via reversible displacement reactions, thereby circumventing the need for intricate regulatory interactions typically required in conventional multistable systems.

We constructed a large-scale CTMN by incorporating additional pre-designed DNA duplexes into the basic two-module network that consists of four pre-designed DNA duplexes, as illustrated in Fig. 5a and Supplementary Fig. 40. Specifically, a three-module network, CTMN-2, is composed of six pre-designed DNA duplexes, CC’, CD’, CE’, DC’, DD’, and DE’. Each duplex features the common exchangeable double-stranded domain (ZZ’) as CTMN-1 and a distinct template domain, P’3 in C’, P’4 in D’, and P’5 in E’. Strand C was designed to activate the polymerization/nicking reactions on CC’, CD’, and CE’, whereas these reactions on DC’, DD’, and DE’ are inhibited owing to the 3′-tether modification on strand D. This structural configuration results in three positive autoregulatory modules: CC’+DD’+DE’, CD’+DC’+DE’, and DC’+DD’+CE’, each two of which engage in a reversible displacement reaction. To monitor the dynamic behavior of CTMN-2, we labeled strand C with quencher BHQ2, and strands C’, D’, and E’ with fluorophores FAM, ROX, and Cy5, respectively, as provided in Supplementary Table 3.

We initiated CTMN-2 with six representative compositions with the dominant module ranging from DD’+CC’+DE’ to CD’+DC’+DE’ and finally to DC’+DD’+CE’, and investigated the corresponding multistable behaviors. Figure 5b illustrates the time-dependent concentration changes of CC’, CD’, CE’, DC’, DD’, and DE’ in CTMN-2 under these conditions, which were derived from the corresponding fluorescence profiles (Supplementary Fig. 41) using appropriate calibration curves (Supplementary Figs. 42 and 43). Regardless of the initial compositions, the concentrations of all six DNA duplexes remain consistent over time, defining six distinct steady states within CTMN-2 (Fig. 5b, c). In contrast, control experiments in the absence of DNA polymerase and nicking endonuclease reveal that, regardless of the initial conditions, the network converges to a common thermodynamically equilibrated state within approximately 4 h (Supplementary Fig. 44), underscoring the key role of the autoregulatory loop in continuously tunable multistability. More interestingly, the steady states of CTMN-1 are linearly distributed, while those of CTMN-2 exhibit a more complex planar arrangement (Fig. 5d). To validate this observation, we conducted additional experiments aimed at more comprehensively mapping the steady-state landscape of CTMN-2 (Supplementary Fig. 45). These studies identify two distinct steady states residing within the interior of the triangular phase space—rather than along its edges—thus providing compelling evidence that CTMN-2 enables higher-order multistability. These results demonstrate the network architecture’s modular scalability. By incorporating additional components, CTMNs with higher-order continuously tunable multistability could be envisaged without the need to design complex regulatory interactions. This modular design broadens the range of accessible steady states and offers a versatile platform for exploring complex dynamic behaviors in synthetic systems.

State-dependent regulation of downstream processes
Next, we investigated the CTMN’s capability to generate state-specific outputs for downstream regulation. Traditional synthetic multistable systems typically depend on intricate signaling cascades to maintain stability, rendering their signaling components unsuitable for direct use in downstream modules without risking destabilization of the core network. As a result, prior designs have often required auxiliary elements—such as signal relays or converter domains—to mediate output functionality at the expense of added circuit complexity61. In contrast, our CTMN architecture supports direct, programmable production of distinct, state-resolved outputs without compromising system stability. This is fulfilled by a core design feature: the multistability of the CTMN is functionally decoupled from the identities and concentrations of the ssDNA species (Px) generated via the polymerization/nicking cycle. The decoupling allows for the ssDNA products to be freely designed with predefined sequences and structural motifs tailored for desired downstream functions. More importantly, state transitions within the CTMN directly modulate the concentrations of network modules that gate DNA replication, thereby exerting dynamic control over the identities and production rates of ssDNA outputs. The network further supports tunable regulation across a continuous spectrum, from full activation to complete repression, rather than enforcing binary ON/OFF responses, enhancing the versatility of the system for interfacing with complex downstream processes. Collectively, these features enable the CTMN to act as a programmable multistable controller that generates modular, state-dependent outputs for downstream reaction networks, without undermining its core dynamical properties.

To demonstrate state-dependent DNA replication, we took CTMN-1 as a case study (Fig. 6a). Transitions between the steady states of CTMN-1 result in significant shifts in the concentrations of modules AA’+BB’ and BA’+AB’. For example, from steady states S1-1 to S1-5, the network composition with AA’+BB’ dominated is gradually altered to that with a majority of BA’+AB’. These module concentration variations could modulate their DNA polymerization/nicking activities to directly control the production of ssDNAs P1 and P2. In state S1-1, the network with predominant AA’+BB’ and minimal BA’+AB’ leads to a large production of P1 but a negligible generation of P2. As the network’s state transitions from S1-1 to S1-5, the P1 production progressively decreases, while the P2 production increases. When the network reaches state S1-5, the P1 production becomes nearly imperceptible, while P2 emerges as the dominant product. This state-dependent modulation of the P1 and P2 production was validated by electrophoretic PAGE analysis (Fig. 6b). In state S1-1, the band corresponding to P1 is intense, whereas band P2 is barely visible, indicating the major production of P1 and the minimal yield of P2. As the network shifts from S1-1 to S1-5, the intensity of band P1 progressively decreases while band P2 intensifies, implying a gradual reduction in the P1 yield and an increasing yield of P2. By state S1-5, P1 is barely detectable, and P2 becomes the predominant product.

To explore the CTMN’s capacity for state-specific downstream regulation, we integrated CTMN-1 with functional modules, such as DNAzymes and T7 RNA transcription systems. Figure 6a illustrates the integration of CTMN-1 with a downstream biocatalytic module consisting of two switchable, Mg2+-ion-dependent DNAzymes. This module was designed with four DNAzyme subunits, M1, M2, N1, and N2. Subunits M1 and M2 interact with P1 generated by module AA’+BB’ to form a catalytically active DNAzyme, DNAzyme 1, while N1 and N2 interact with P2 produced by module BA’+AB’ to yield DNAzyme 2. Each DNAzyme cleaves a specific fluorophore/quencher-labeled substrate (Sub1 for DNAzyme 1 and Sub2 for DNAzyme 2) to produce distinct fluorescence signals.

We monitored the time-dependent fluorescence changes generated from DNAzyme-catalyzed cleavage reactions under different states of CTMN-1 (Supplementary Figs. 46−48). Using calibration curves that correlate the catalytic activity of each DNAzyme with its concentration (Supplementary Fig. 49), we quantified the state-dependent concentration changes of DNAzymes 1 and 2 across states S1-1 to S1-5 (Fig. 6c). In state S1-1, DNAzyme 1 exhibits maximal yield, while the yield of DNAzyme 2 remains minimal (Supplementary Fig. 47). This outcome is consistent with the network’s DNA replication dynamics, where the P1 production predominates while the P2 production is negligible. As the system transitions from S1-1 to S1-5, the yield of DNAzyme 1 decreases to a negligible level, whereas the yield of DNAzyme 2 increases correspondingly, paralleling the shift from the dominant P1 to P2 production within the network. By state S1-5, DNAzyme 1 achieves minimal yield, while the yield of DNAzyme 2 reaches its maximum (Supplementary Fig. 48). These results reveal the CTMN-1’s robust capacity to integrate with downstream DNAzyme modules, enabling fine-tuned, state-dependent regulation of biocatalytic processes.

Furthermore, we explored whether the sequences of ssDNAs produced by the network could be programmed with specific functionalities to regulate diverse biological processes in a state-dependent manner. As illustrated in Fig. 6d, we integrated CTMN-1 with T7 RNA polymerase (T7RNAP) transcriptional machinery73,76,77 as an exemplary system. The T7RNAP transcriptional machinery involves two switchable transcriptional units, F1 and F2, each containing a DNA duplex region encoding an RNA strand (F1 for rR1 and F2 for rR2) and an incomplete promoter domain for regulation. To couple CTMN-1 and the T7RNAP transcriptional machinery, template domains P’1 and P’2 in CTMN-1 were substituted by template domains P’6 and P’7, resulting in the production of two ssDNAs, P6 and P7, for activating transcriptional units F1 and F2, respectively. That is, in the absence of P6 and P7, transcriptional units F1 and F2 with the incomplete promoter domain could lead to an inefficient transcription of rR1 and rR2, respectively. However, as CTMN-1 functions in distinct steady states, variable concentrations of P6 and P7 can be produced to selectively activate transcriptional units F1 and F2 for state-dependent transcription of rR1 and rR2. To improve the signal-to-noise ratio of P6- and P7-stimulated activation of F1 and F2, threshold levels of single strands P’6 and P’7 were added to the integrated system.

For real-time monitoring of RNA transcription, we employed two fluorophore/quencher-labeled DNA duplexes, X1X’1 and X2X’2, as reporters (Fig. 6e, top panel), which can hybridize with rR1 and rR2 to generate fluorescence signals, respectively. We measured the time-dependent fluorescence changes corresponding to RNA transcription of the integrated system across five steady states, from S1-1 to S1-5 (Fig. 6e and Supplementary Fig. 50). In state S1-1, the fluorescence signal associated with the rR1 transcription reaches a maximum, whereas the signal from the rR2 transcription remains nearly undetectable. This result is consistent with the network in state S1-1, primarily producing P6, which activates transcriptional unit F1, and minimally generating P7, which activates F2. As the network transitions from S1-1 to S1-5, the fluorescence signal related to the rR1 transcription decreases, while the signal corresponding to the rR2 transcription increases. This agrees with the decrease in P6 production and the increase in P7 production in CTMN-1, leading to a decrease in F1 activation and an increase in F2 activation. By S1-5, the network primarily generates P7, resulting in an inefficient rR1 transcription and a maximal rR2 transcription. Importantly, integrating CTMN-1 with the RNA transcription module does not compromise the system’s multistability (Supplementary Fig. 51). These results highlight the programmability of ssDNAs produced by the network, enabling direct, state-dependent regulation of nucleic acid-responsive biochemical and biological processes. The ability to design customized DNA sequences to modulate downstream processes demonstrates the CTMN’s versatility in controlling complex biochemical systems.

Methods
Materials
Magnesium chloride hexahydrate (MgCl2·6H2O) and tris(hydroxymethyl)aminomethane hydrochloride (Tris-HCl) were purchased from Aladdin Chemical Reagent Co., Ltd. Hydrochloric acid (HCl) was sourced from Sinopharm Chemical Reagent Co., Ltd. Nt.BbvCI nicking endonuclease (10,000 U mL−1, 2.3 μM), RNase H (5000 U mL−1), 10× rCutSmart® buffer (1× buffer components: 50 mM potassium acetate, 20 mM Tris-acetate, 10 mM magnesium acetate, 100 μg mL−1 recombinant protein, pH 7.9@25°C), Klenow Fragment (3→5 exo-) (5000 U mL−1), deoxyribonucleotide triphosphates (dNTPs), and NEBuffer 2 (1× buffer components: 50 mM NaCl, 10 mM Tris-HCl, 10 mM MgCl2, 1 mM dithiothreitol, pH 7.9@25°C) were purchased from New England BioLabs Inc. Ammonium persulfate substitute (APS substitute), Tris-borate-EDTA (TBE, 5× premixed powder), Gel-Red (10,000×), acrylamide/bis-acrylamide (29:1, 30% solution), tetramethylethylenediamine (TEMED) substitute and DNA loading buffer were purchased from Beyotime Biotechnology. Co., Ltd. T7 High-yield RNA synthesis kit was purchased from Yeasen Biotechnology (Shanghai) Co., Ltd. Ultrapure water (electric resistance, 18.2 MΩ cm) obtained through a Millipore Milli-Q water purification system (Billerica, MA, USA) was used throughout the experiments. All DNA and RNA strands were purchased from Sangon Biotech Co., Ltd. (Shanghai, China).

The oligonucleic acid sequences used in the study are the following:
(1) A: 5′-AGCAATAACTAAC-BHQ2-ACAACTACTTTCC-3′
(2) A’: 5′-ATTGTTGTTGAGAGAAATGCTGAGGAAAGTAGTTGT-FAM-ATAAAGGTCGAATTG-3′
(3) B: 5′-ACAACTACTTTCCAATAATA-3′
(4) B’: 5′-TGTTGAATGATGAGAAATGCTGAGGAAAGTAGTTGT-ROX-AAGATGAAGTGGAA-3′
(5) M1: 5′-ATTGTTGTTGCTAATCACCCATGTTACCTG-3′
(6) M2: 5′-CTCTTCAGCGATATTAGGAGAAATGC-3′
(7) P1: 5′-TCAGCATTTCTCTCAACAACAAT-3′
(8) N1: 5′-TGTTGAATGATTTAATTTCACCCATGTTTCTGT-3′
(9) N2: 5′-GACCTCAGCGATAAATTAAGAGAAATG-3′
(10) P2: 5′-TCAGCATTTCTCATCATTCAACA-3′
(11) rTAA’: 5′-CAAUUCGACCUUUAUGUUAGUUAUUGCU-3′
(12) rTAB’: 5′-UUCCACUUCAUCUUGUUAGUUAUUGCU-3′
(13) A’ (without template domain P’1): 5′-GGAAAGTAGTTGT-FAM-ATAAAGGTCGAATTG-3′
(14) B’ (without template domain P’2): 5′-GGAAAGTAGTTGT-ROX-AAGATGAAGTGGAA-3′
(15) C: 5′-TGATCTATTATAC-BHQ2-ACAACTACTTTCC-3′
(16) C’: 5′-GTTGTTATTGTTGTTGAGAGAAATGCTGAGGAAAGTAGTTGT-FAM-ATAAAGGTCGAATTG-3′
(17) D: 5′-ACAACTACTTTCCAAAAATA-3′
(18) D’: 5′-GTAGTGTTAGTAAGATAGAGAAATGCTGAGGAAAGTAGTTGT-ROX-AAGATGAAGTGGAA-3′
(19) E’: 5′-TGAGTATTGGTTTAGTTGAGAAATGCTGAGGAAAGTAGTTGT-Cy5-ATTAAAGAGGAGAG-3′
(20) J1: 5′-GTTGTTATTGAGATTCACCCATGTTACCTG-3′
(21) J2: 5′-CTCTTCAGCGATAATCTTTGTTGAGAG-3′
(22) L1: 5′-TCAGCATTTCTCTCAACAACAATAACAAC-3′
(23) K1: 5′-GTAGTGTTAGATTAATTCACCCATGTTTCTGT-3′
(24) K2: 5′-GACCTCAGCGATAATTAATAAGATAGAGA-3′
(25) L2: 5′-TCAGCATTTCTCTATCTTACTAACACTAC-3′
(26) O1: 5′-TGAGTATTGGATCTTCACCCATGTTTCTGT-3′
(27) O2: 5′-CTTGTCAGCGATAAGATTTTAGTTGAG-3′
(28) L3: 5′-TCAGCATTTCTCAACTAAACCAATACTCA-3′
(29) Sub1: 5′-ROX-CAGGTATrAGGAAGAG-BHQ2-3′
(30) Sub2: 5′-FAM-ACAGAATrAGGAGGTC-BHQ1-3′
(31) Sub3: 5′-FAM-ACAGAATrAGGACAAG-BHQ1-3′
(32) A-noQ: 5′-AGCAATAACTAACTACAACTACTTTCC-3′
(33) A’(P’6): 5′-AGATGATATGTGTTAATACGATGCTGAGGAAAGTAGTTGT-TATAATGGTGGAAATG-3′
(34) A’(P’6)-FAM: 5′-AGATGATATGTGTTAATACGATGCTGAGGAAAGTAGTTGT-FAM-TATAATGGTGGAAATG-3′
(35) B’(P’7): 5′-TGTTGAATGATGATAATACGATGCTGAGGAAAGTAGTTGT-TATGATGATGTGGAA-3′
(36) B’(P’7)-ROX: 5′-TGTTGAATGATGATAATACGATGCTGAGGAAAGTAGTTGT-ROX-TATGATGATGTGGAA-3′
(37) Nt1: 5′-AGATGATATGTGTTAATACGACTCACTATAGGGAAGGTAAGTCTCGTATCAGCACTGACAAAGTCAGAAA-3′
(38) T1: 5′-TTTCTGACTTTGTCAGTGCTGATACGAGACTTACCTTCCCTATAGTGAGTAATTTAA-3′
(39) Nt2: 5′-TGTTGAATGATGATAATACGACTCACTATAGGGTGTTCTCTGTGTTGTTGTATCTGACAAAGTCAGAAA-3′
(40) T2: 5′-TTTCTGACTTTGTCAGATACAACAACACAGAGAACACCCTATAGTGAGTAATTTAA-3′
(41) X1: 5′-Cy5-ATTGCTGATACGAGACTTACCTTC-3′
(42) X’1: 5′-GAAGGTAAGATCAGCAAT-BHQ2-3′
(43) X2: 5′-FAM-CAATACAACAACACAGAGAACAACC-3′
(44) X’2: 5′-GGTTGTTCTCGTTGTATTG-BHQ1-3′
(45) P’6: 5′-AGATGATATGTGTTAATACGATGCTGA-3′
(46) P’7: 5′-TGTTGAATGATGATAATACGATGCTGA-3′
(47) A-tether: 5′-CATAGATTTGACATAGAAGAGCAATAACTAACTACAACTACTTTCC-3′
(48) A’-tether: 5′-GGTGTACTGAATATTGTTGTTGAGAGAAATGCTGAGGAAAGTAGTTGTTATAAAGGTCGAATTG-3′
(49) B’-no F: 5′-TGTTGAATGATGAGAAATGCTGAGGAAAGTAGTTGT-AAGATGAAGTGGAA-3′
(50) P-tether1: 5′-TCAGCATTTCTCTCAACAACAATATTCAGTACACC-3′

The ribonucleobase cleavage site, rA, in the substrates of the different Mg2+-ion-dependent DNAzymes is indicated in bold, the respective Mg2+-ion-dependent DNAzyme sequences are underlined, and the sequence domains recognized by Nt.BbvCI is presented in italics and in bold.

Instruments and measurements
Fluorescence spectra were recorded with a Cary Eclipse Fluorometer. FAM, ROX, and Cy5 excitations were performed at 488, 580, and 635 nm, respectively. FAM, ROX, and Cy5 emissions were recorded at 520, 608, and 668 nm, respectively. PAGE Gel electrophoresis was performed on a Hoefer SE600 electrophoresis unit.

Preparation of duplexes AA’, AB’, BA’, BB’, CC’, CD’, CE’, DC’, DD’, DE’
Taking AA’ as an example, a mixture of A and A’ (500 nM each) in 1× rCutSmart buffer was annealed at 40°C for 10 min, cooled down to 25°C at a rate of 0.3°C min−1, and allowed to equilibrate at 25°C for 2 h to yield duplex AA’.

Probing multistable behaviors in CTMN-1
To demonstrate multistability in the network, seven steady states (S1-1 – S1-7) were established by adjusting the initial composition of the network. All the assays were prepared in 100 μL of 1× rCutsmart buffer. For the establishment of steady state S1-1, a mixture of AA’ (100 nM) and BB’ (100 nM) was treated with DNA polymerase (30 U mL−1), Nt.BbvCI (20 U mL−1) and dNTPs (0.1 mM), that is, the network was initiated with a composition of 100% AA’+BB’ and 0% BA’+AB’. The fluorescence changes of the resulting mixtures were monitored in real-time utilizing a fluorometer at 37°C. Different initial compositions were used to establish other steady states. For S1-2, the network was initiated with a mixture of AA’ (75 nM), BB’ (75 nM), AB’ (25 nM), and BA’ (25 nM), a composition of 75% AA’+BB’ and 25% BA’+AB’; For S1-3, the network was initiated with a mixture of AA’ (50 nM), BB’ (50 nM), AB’ (50 nM), and BA’ (50 nM), a composition of 50% AA’+BB’ and 50% BA’+AB’; For S1-4, the network was initiated with a mixture of AA’ (25 nM), BB’ (25 nM), AB’ (75 nM), and BA’ (75 nM), a composition of 25% AA’+BB’ and 75% BA’+AB’; For S1-5, the network was initiated with a mixture of AB’ (100 nM) and BA’ (100 nM), a composition of 0% AA’+BB’ and 100% BA’+AB’; For S1-6, the network was initiated with a mixture of AA’ (65 nM), BB’ (65 nM), AB’ (35 nM), and BA’ (35 nM), a composition of 65% AA’+BB’ and 35% BA’+AB’; For S1-7, the network was initiated with a mixture of AA’ (35 nM), BB’ (35 nM), AB’ (65 nM), and BA’ (65 nM), a composition of 35% AA’+BB’ and 65% BA’+AB’. It should be noted that, unless otherwise specified, the concentrations of DNA polymerase, Nt.BbvCI and dNTPs were fixed at 30 U mL−1, 20 U mL−1, and 0.1 mM, respectively.

Probing orthogonal and continuous state transitions of CTMN-1
For orthogonal state transitions, samples of CTMN-1 in steady states S1-1 and S1-5 were prepared as before and incubated at 37°C for 4 h. Subsequently, RNA triggers rTAB’ (200 nM) and rTAA’ (200 nM) were introduced into the prepared samples, respectively. After 4 h, each of the resulting mixtures was treated with RNase H (5.00 U mL−1). The fluorescence changes of the resulting mixtures were recorded in real-time at 37°C.

For continuous state transitions, the prepared samples of CTMN-1 in steady states S1-1 and S1-5 were subjected to rTAB’ (1 μM)/RNase H (3.75 U mL−1) and rTAA’ (1 μM)/RNase H (3.75 U mL−1), respectively. Every 2 h, a second and a third round of rTAB’ (1 μM) and rTAA’ (1 μM) were added, respectively. The fluorescence changes of the resulting mixtures were monitored in real-time at 37°C.

Probing multistable behaviors in CTMN-2
To demonstrate multistability in CTMN-2, six steady states (S2-1 – S2-6) were established in network CTMN-2 by adjusting the initial compositions of the network. All the assays were prepared in 100 μL of 1× rCutsmart buffer. For steady state S2-1, a mixture of CC’ (100 nM), DD’ (100 nM), and DE’ (100 nM) was treated with DNA polymerase (100 U mL−1), Nt.BbvCI (50 U mL−1) and dNTPs (0.1 mM), that is, the network was initiated with a composition of 100% DD’+CC’+DE’, 0% CD’+DC’+DE’, and 0% DC’+DD’+CE’. The fluorescence changes of the resulting mixtures were monitored in real-time using a fluorometer at 38°C. Different initial compositions were used to establish the remaining steady states. For S2-2, the network was initiated with a mixture of CD’ (100 nM), DC’ (100 nM), and DE’ (100 nM), a composition of 0% DD’+CC’+DE’, 100% CD’+DC’+DE’, and 0% DC’+DD’+CE’; For S2-3, the network was initiated with a mixture of DC’ (100 nM), DD’ (100 nM), and CE’ (100 nM), a composition of 0% DD’+CC’+DE’, 0% CD’+DC’+DE’, and 100% DC’+DD’+CE’; For S2-4, the network was initiated with a mixture of CC’ (50 nM), CD’ (50 nM), DE’ (100 nM), DC’ (50 nM), and DD’ (50 nM), a composition of 50% DD’+CC’+DE’, 50% CD’+DC’+DE’, and 0% DC’+DD’+CE’; For S2-5, the network was initiated with a mixture of CC’ (50 nM), CE’ (50 nM), DD’ (100 nM), DC’ (50 nM) and DE’ (50 nM), a composition of 50% DD’+CC’+DE’, 0% CD’+DC’+DE’, and 50% DC’+DD’+CE’; For S2-6, the network was initiated with a mixture of CD’ (50 nM), CE’ (50 nM), DC’ (100 nM), DD’ (50 nM) and DE’ (50 nM), a composition of 0% DD’+CC’+DE’, 50% CD’+DC’+DE’, and 50% DC’+DD’+CE’; For S2-7, the network was initiated with a mixture of CC’ (40 nM), CD’ (30 nM), CE’ (30 nM), DC’ (60 nM), DD’ (70 nM) and DE’ (70 nM), a composition of 40% DD’+CC’+DE’, 30% CD’+DC’+DE’, and 30% DC’+DD’+CE’; For S2-8, the network was initiated with a mixture of CC’ (20 nM), CD’ (60 nM), CE’ (20 nM), DC’ (80 nM), DD’ (40 nM) and DE’ (80 nM), a composition of 20% DD’+CC’+DE’, 60% CD’+DC’+DE’, and 20% DC’+DD’+CE’.

Probing state-dependent biocatalytic reactions of DNAzymes mediated by CTMN-1
To assess state-dependent biocatalytic activities of DNAzymes 1 and 2 dictated by CTMN-1, samples of CTMN-1 in five steady states (S1-1 – S1-5) were prepared as before and incubated at 37°C for 12 h. For the state-dependent biocatalytic reaction of DNAzyme 1, 88 μL of the prepared samples were then treated with a mixture of M1 (7.5 μL of 10 μM), M2 (7.5 μL of 10 μM), and MgCl2·6H2O (1.51 μL of 1 M) with final concentrations to be 0.5 μM, 0.5 μM, and 20 mM, respectively. After the incubation at 25°C for 2 h, Sub1 (3 μL of 100 μM) was added to these samples. Subsequently, the time-dependent fluorescence changes generated by DNAzyme 1-catalyzed cleavage of Sub1 were monitored at 25°C. For DNAzyme 2, N1 (7.5 μL of 10 μM), N2 (7.5 μL of 10 μM), and Sub2 (3 μL of 100 μM) were used instead of M1 (7.5 μL of 10 μM), M2 (7.5 μL of 10 μM) and Sub1 (3 μL of 100 μM), respectively.

Preparation of switchable T7 RNAP RNA transcription module
The switchable T7 RNAP RNA transcription module consists of incomplete transcriptional units F1 and F2 and reported units X1X’1 and X2X’2. For the preparation of incomplete transcriptional units F1 (Nt1T1) and F2 (Nt2T2), a mixture of Nt, T, Nt, and T, each of 5 µM, in 20 mM Tris-HCl buffer, pH=7.9, including 10 mM MgCl2·6H2O, was annealed at 85°C, gradually cooled to 25°C at a rate of 1.0°C min−1, and equilibrated at 25°C for 2 h. For the preparation of reported units X1X’1 (for rA transcription) and X2X’2 (for rA transcription), a mixture of X, X’, X and X’, each of 40 µM, in 20 mM Tris-HCl buffer, pH=7.9, including 10 mM MgCl2·6H2O was annealed at 85°C, gradually cooled to 25°C at a rate of 1.0°C min−1, and equilibrated at 25°C for 2 h.

Probing state-dependent RNA transcription mediated by modified CTMN-1
Samples of modified CTMN-1 in five steady states (S1-1 – S1-5) were prepared by adjusting the initial concentrations of duplexes as before, in which A-noFQ, A’(P’6), and B’(P’7) were used instead of A, A’, and B’, respectively. The resulting samples were treated with the prepared T7 RNAP RNA transcription module, forming an integrated system capable of state-dependent RNA transcription. As an example, the integrated system in steady state S1-1 was prepared by mixing A-noFQ/A’(P’6) (100 nM), BB’(P’) (100 nM), DNA polymerase (35 U mL−1), Nt.BbvCI (30 U mL−1), dNTPs (0.4 mM), F1 (100 nM), F2 (100 nM), X1X’1 (400 nM), X2X’2 (400 nM), P’6 (30 nM), P’7 (80 nM), NTPs (1.43 mM) and T7 RNAP (72 U mL−1) in 1× rCutSmart buffer. To evaluate state-dependent RNA transcription, fluorescence signals generated from the integrated system across various steady states were monitored in real-time using a fluorometer at 37°C.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this work, the authors present a reaction network that uses competing displacement reactions and a DNA polymerization/nicking module to create an \"arbitrarily\" multi-stable system. They build up their system step-by-step:\n1. Starting with a simple equilibrium between two competing displacement reactions, balancing two states (AA' + BB' and AB' + BA').\n2. Adding a nicking/polymerase sub-system to one of the two states, creating a far-from-equilibrium bi-stable system.\n3. Implementing the nicking/polymerase sub-system for both states, leading to the desired multi-stable system. The system self-regulates to create a steady state at the initially present concentrations of the two states, whatever this may be.\n4. Adding a submodule that uses complementary RNA strands to stabilize either state, which can be removed using RNase H, effectively allowing the shift of the relative amount of each state.\nThe authors also show a more complex system with the same sub-systems that has three competing states instead of two. Furthermore, the authors implement two functional downstream reactions that do not influence the multi-stable system itself.\nThe work is novel, showing a new path to multi-stability which does not require increasingly complex reaction networks to add additional states, and instead uses a more straightforward and modular approach to achieving multi-stability through competing displacement reactions and DNA polymerization/nicking modules.\nHowever, I am somewhat concerned about the fundamental thesis of this manuscript: is \"arbitrary multi-stability\" truly useful? Many of the interesting properties of bi- and multi-stability are caused by the discrete nature of the stable states: switching behavior, hysteresis, oscillations, stabilization under perturbations, memory, etc. More stable points are interesting, sure, but when transitions are continuous, like in this system, have we not lost much of the interesting behavior of these systems? What exactly is the added benefit of this system compared to mixing two DNA strands directly without using a complex dissipative system? In my opinion, this point must be properly addressed in the manuscript.\nOn top of that I have some other major points:\nI would like to see at least one steady state for the three-module network that is not at the edge of the triangle in Fig. 5d.\nIn all cases I see here, the modules (let's say AA' + BB') start out balanced: there is an equal amount of AA' and BB'. Is this required for the system to work? What happens when we create an imbalanced stoichiometry, keeping the total stoichiometry correct by applying an inverse perturbation on the AB'/BA' module?\nPossibly related to that: the data shows that the two components in a module do not always have the same concentrations, sometimes showing quite a big difference, see for instance Fig. 4e. What is the origin of this?\nMany plots showing the development of the steady-states do not seem to be truly steady and show drift of a few percent per hour even after 12 hours (Fig. 3d S1-1 for instance). Is this a question of waiting longer for true stability, are these effects from consumption of dNTPs, aging reagents, or is something else going on?\nThe authors use a kinetic model to predict the behavior of the system, and mention it \"aligns closely\" with the experimental system. However, in several plots it is clear that the model deviates and does not fit even qualitatively (Fig. 3c-v, Fig. 3d-S1-5, Sup. Fig. 9c&d, Sup. Fig. 8d, Sup. Fig. 12a, Sup. Fig. 17f). I realize the deviations are not too big, but can the authors formulate a reason why this is the case?\nIs there a reason why the rTaa/rTab & RNase H reactions are not investigated with the kinetic model? The relevant rate constants seem to be determined in Sup. Fig. 20.\nI wonder about the robustness of the system. The authors convincingly show that a steady state is found at or around the initial composition of the system, but how robust is this when concentrations of the various catalytic components, dNTPs, etc., are varied/consumed? A small discussion in the SI would be appreciated. I realize some data is included (Sup. Fig. 13-18), but some written discussion would be helpful.\nFig. 3d, bottom row, second from the left as well as Fig. 5d (and their equivalent plots in Fig. 1a&b): These plots do not add anything and suggest a higher dimensional system than is actually presented. For Fig. 3, it is really not clear what the message of the plot is supposed to be anyway. In Fig. 5, a more appropriate way of showing this would be a ternary plot, which correctly shows 2 orthogonal dimensions instead of 3 - which the authors do actually use in Fig. 1b.\nMiscellaneous minor points:\n- In my version of the SI, I see black lines around figure edges.\n- Supplementary Fig. 29b has hours displayed far away from the plots they are supposed to represent.\n- If 3D plots are used, one should label their axes; most plots are missing the initial compositions axis (e.g., something like \"fraction AA' + BB'\" in Fig. 3).\n- In Fig. 3, the usage of both a, b, c and i, ii, iii as figure labels independent of each other at the same time is confusing. In general, it is very hard to refer to specific panels: Fig. 3d has 10 (!) unlabeled sub-panels; these should be numbered or labeled somehow. Figure 3 is the main offender, but it is worth going over all figures again.", "reviewer_2": "In this paper, the authors present some observations about a network based on polymerization nicking reaction on two templates competing for a single trigger (where the stoichiometry is balanced by introduction of a non-extendable pseudo-trigger). The polymerization/nickase cycles slow down the trigger exchange reaction, but does not abolishes it on longtime scales, a phenomenon that the author wrongly interpret as a new form of multistability. However, continuous multistablity does not exist mathematically. Slow relaxation does. In addition it is inappropriate to discuss steady states (and hence their stability) in a system where some of the variables (here compounds) can only increase, which is the case here for strand P, which is continuously produced and accumulates. Furthermore mathematical support for multistability may not come from kinetic ODE fitting but from Linear Stability Analysis (LSA), which I encourage the author to try (once they have corrected their kinetic model). If the authors prefer an experimental test, they can also attempt to adjust the balance between polymerization and nicking activities, upon which they will probably find that the system does relax to a single AA’/AB’ ratio, which a time scale that decrease with pol/nick ratio, but this does not mean the system is multistable (the situation presented in Fig 3a being the limit pol/nick->0).", "reviewer_3": "The manuscript by Zhong et al. introduces a novel strategy for constructing multistable biochemical networks using reversible DNA strand displacement combined with polymerization/nicking modules. While the authors propose an ambitious framework for achieving arbitrarily programmable multistability, several critical aspects require clarification to fully support the claims. Notably, the extent and scalability of this \"arbitrarily programmable\" capability—specifically, how many distinct steady states can be robustly engineered and maintained—remain unclear. Regrettably, the current evidence does not sufficiently validate the authors' central thesis. Below are my specific concerns:\n1. The current schematic is challenging to interpret. To improve accessibility, the authors should explicitly label the 3′ and 5′ ends of DNA strands to indicate polymerase extension directions. Additional molecular details in the diagram would help readers visualize the design principles.\n2. The domain-level labels on DNA complexes are insufficiently clear. Revising the figures to highlight domain names with distinct colors or fonts would aid reader comprehension.\n3. The description leaves uncertainty about whether the four duplexes (AA', BB', BA', AB') are pre-annealed or formed spontaneously from single strands. Clarifying this hierarchy in the text is essential for reconstructing the system’s operational logic.\n4. The distinction between AA'-P1/AB'-P2 and AA'/AB' remains unclear. If these complexes are functionally equivalent, consistent naming would reduce confusion. If not, a side-by-side comparison of their structural differences should be provided.\n5. The manuscript does not sufficiently explain the driving forces behind the reversible transitions (AA'+BB' -- BA'+AB'). Specifically: What thermodynamic or kinetic factors enable bidirectional conversion? How does AA'P1 formation selectively inhibit the forward reaction? A mechanistic diagram or kinetic model would help justify this design.\n6. The assertion that dNTP hydrolysis suppresses spontaneous interconversion lacks experimental or theoretical support in the \"Design\" section. The authors should explicitly describe how dNTP consumption is coupled to module stabilization.\n7. The roles of strands A (BHQ2-labeled), A' (FAM), and B' (ROX) are ambiguously described. A dedicated subsection or table mapping strands to their modifications would prevent misinterpretation of the monitoring strategy.\n8. The significant difference in equilibration kinetics between the two initial conditions (100% AA'+BB' vs. BA'+AB') appears inconsistent with the proposed thermodynamic framework. The authors should address whether this asymmetry arises from kinetic traps or unintended design biases. (Fig. 3a)\n9. The claim that rTAA' hybridization upregulates AA'+BB' requires mechanistic elaboration. How does AA'-rTAA' formation preferentially stabilize specific modules? A diagram illustrating the energy landscape or binding competition would strengthen this argument.\n10. The manuscript would benefit from exploring scenarios involving simultaneous or sequential RNA inputs: How would competing RNA triggers at varying concentrations influence state selection? Is the system fully reversible if a second RNA is introduced after rTAA'-induced stabilization?\n11. The observed convergence near initial concentrations raises questions about true multistability. To substantiate claims, the authors should: Define quantitative criteria for distinguishing steady states. Compare experimental results with theoretical predictions of bistable/multistable systems\n12. The gel lanes for DNAzyme 2 products lack sufficient contrast (Fig. 6b). Repeating the experiment with adjusted exposure times or alternative visualization methods (e.g., SYBR Gold staining) would improve credibility."}, "author_response": "Re: Manuscript ID: NCOMMS-25-07099-T\nTitle: \"Arbitrary Multistability in DNA Replication Networks\"\nThe following point-by-point listed corrections/explanations were introduced into\nthe paper:\nReviewer #1:\nComment 1: “In this work, the authors present a reaction network that uses competing\ndisplacement reactions and a DNA polymerization/nicking module to create an\n\"arbitrarily\" multi-stable system. They build up their system step-by-step:\n1. Starting with a simple equilibrium between two competing displacement reactions,\nbalancing two states (AA' + BB' and AB' + BA').\n2. Adding a nicking/polymerase sub-system to one of the two states, creating a far-\nfrom-equilibrium bi-stable system.\n3. Implementing the nicking/polymerase sub-system for both states, leading to the\ndesired multi-stable system. The system self-regulates to create a steady state at the\ninitially present concentrations of the two states, whatever this may be.\n4. Adding a submodule that uses complementary RNA strands to stabilize either state,\nwhich can be removed using RNase H, effectively allowing the shift of the relative\namount of each state.\nThe authors also show a more complex system with the same sub-systems that has three\ncompeting states instead of two. Furthermore, the authors implement two functional\ndownstream reactions that do not influence the multi-stable system itself.\nThe work is novel, showing a new path to multi-stability which does not require\nincreasingly complex reaction networks to add additional states, and instead uses a\nmore straightforward and modular approach to achieving multi-stability through\ncompeting displacement reactions and DNA polymerization/nicking modules.”\nResponse: We sincerely thank the reviewer for their thorough and insightful summary\nof our work, which highlights the conceptual framework, core design principles, and\nstep-by-step experimental development of our system. We are especially grateful for\nthe reviewer’s recognition that our work represents a novel and scalable approach to\nengineering multistable biochemical networks. We believe that the design principles\ndemonstrated here offer a robust foundation for developing synthetic molecular systems\nwith programmable and tunable multistable behavior.\nComment 2: “However, I am somewhat concerned about the fundamental thesis of this\nmanuscript: is \"arbitrary multi-stability\" truly useful? Many of the interesting\nproperties of bi- and multi-stability are caused by the discrete nature of the stable states:\nswitching behavior, hysteresis, oscillations, stabilization under perturbations, memory,\netc. More stable points are interesting, sure, but when transitions are continuous, like\nin this system, have we not lost much of the interesting behavior of these systems? What\nexactly is the added benefit of this system compared to mixing two DNA strands directly\nwithout using a complex dissipative system? In my opinion, this point must be properly\naddressed in the manuscript.”\nResponse: We sincerely thank the reviewer for raising this fundamental and thought-\nprovoking question regarding the conceptual significance and practical utility of\narbitrary multistability. We fully agree that classical bistable and multistable systems—\ncharacterized by discrete attractor states—play a central role in both natural and\nsynthetic biochemical circuits. These systems enable a wide range of behaviors,\nincluding hysteresis, memory, switch-like transitions, and robustness to perturbations.\nSuch properties have been instrumental in enabling binary decisions and nonlinear\ninformation processing in diverse molecular and cellular contexts (e.g., genetic toggle\nswitches and cell fate determination, ref. 50—Proc. Natl Acad. Sci. USA 106, 6638–\n6643 (2009) and ref. 51—Science 375, eabg9765 (2022)).\nIn contrast, our study explores a fundamentally different and complementary\nparadigm: continuous arbitrary multistability, where the system stabilizes not to a\nfinite set of discrete states, but to any compositional states along a continuous spectrum.\nWhile discrete multistability is optimal for binary or categorical functions, many\nbiological processes require the maintenance or encoding of graded, analog information.\nA salient example is the neural control of eye position, where the oculomotor integrator\nsustains arbitrary eye angles across a continuum (ref. 68—Proc. Natl. Acad. Sci. USA\n93, 13339–13344 (1996), ref. 69—Nat. Neurosci. 5, 775–782 (2002) and ref. 70—Curr.\nOpin. Neurobiol. 13, 204–211 (2003)). This function is achieved via neural circuits with\ncontinuous attractor dynamics, enabling persistent memory of analog inputs and precise\nmotor control.\nOur AMN aims to mimic this mode of operation at the molecular level. Specifically,\nit allows the network to stabilize at any desired compositional state within a continuous\nrange through programmable chemical interactions, allowing reversible and continuous\ntransitions between steady states. Unlike traditional multistable systems with fixed\nswitching thresholds, the AMN architecture supports smooth, analog modulation and\nmemory, thereby endowing it with new functional modalities such as graded signal\nprocessing, proportional memory, and biomimetic analog computation—\ncapabilities not accessible with traditional discrete-state architectures.\nNotably, although the transitions in our system are continuous, the AMN retains\nseveral key dynamic features often associated with bistability. For example, hysteresis\n(refers to the failure of a system to return to its starting point when an external condition\nor control parameter is changed and then changed back, ref. 52—Nature 537, 656–660\n(2016) and ref. 70—Curr. Opin. Neurobiol. 13, 204–211 (2003)) is observed upon the\nmodulation of the input species (e.g., rT ): once the system transitions to a new state\nAA'\nupon rT addition, the subsequent degradation of rT (e.g., by RNase H) does not\nAA' AA'\nrestore the original state (Fig. 4d). This path-dependent irreversibility reflects\nmolecular memory, wherein the current state encodes the history of prior stimuli, a\nhallmark behavior of multistable systems.\nWith respect to perturbation handling, classical discrete multistable systems tend to\nreject small perturbations by remaining within their attractor basins. In contrast, our\nAMN offers graded adaptability, dynamically adjusting its state in proportion to the\nmagnitude of the input. As shown in Fig. 4f, the system can resolve a range of input\nconcentrations into distinct steady states. We refer to this behavior as concentration-\nadaptive molecular memory, wherein the system encodes and preserves the\nquantitative features of its environment.\nAnother key strength of the AMN architecture lies in its modularity and its ability to\nconnect downstream reactions without interference. In many conventional multistable\nnetworks, interfacing with downstream processes risks perturbing the core dynamics,\nas key intermediates must be rerouted. In contrast, our AMN system operates via a\ndNTP-fueled polymerization/nicking reaction cycle, which consumes dNTPs to\ngenerate single-stranded DNA waste products (e.g., P , P in Fig. 6a). These ssDNA\n1 2\noutputs are functionally decoupled from the multistable core and can thus be\nindependently harnessed to drive downstream reaction modules without perturbing the\nsystem’s multistability. This waste-to-signal conversion strategy improves the\nmodularity and functionality of synthetic molecular networks.\nFinally, we emphasize that the observed multistability is fundamentally dependent\non the dissipative, far-from-equilibrium nature of the AMN. A thermodynamically\npassive mixture of DNA duplexes, lacking the energy-driven polymerase/nicking cycle,\nwould invariably relax to a single thermodynamically equilibrated state, regardless of\nthe initial composition (Fig. 3b and 3c and new Supplementary Fig. 5). In contrast, the\ncontinuous energy input in our system sustains multiple, programmable steady states\n(new Supplementary Figs. 13−15, see Comment 6) — an essential feature for\nmimicking the far-from-equilibrium behaviors observed in biological counterparts.\nIn response to this insightful comment, we have clarified these conceptual and\nfunctional distinctions between discrete and continuous multistability in the revised\nmanuscript (pages 4, 21, and 22). Also, we have included a discussion of the dNTP-\ndriven operation of the system in the revised manuscript (pages 5, 6, and 11) and in the\nnew Supplementary Section 5. Furthermore, we have expanded the discussion to clarify\nthe unique advantages of the AMN framework in the revised manuscript (page 18).\nSupplementary Fig. 5. Dynamic behaviors of network AMN-1 in the absence of DNA\npolymerase, nicking endonuclease, and dNTPs. a−e, Time-dependent fluorescence changes\nassociated with network AMN-1 with five different initial compositions: 100% AA'+BB' and\n0% BA'+AB' (a); 75% AA'+BB' and 25% BA'+AB' (b); 50% AA'+BB' and 50% BA'+AB' (c);\n25% AA'+BB' and 75% BA'+AB' (d); and 0% AA'+BB' and 100% BA'+AB' (e). f−h, Time-\ndependent concentration changes of AA', AB', BA', and BB' derived from fluorescence data\nshown in panels b−d by applying the corresponding calibration curves (Supplementary Fig. 4).\ni, Compositional states of the systems after 12 h, represented as the fraction of module BA'+AB'\nunder the five initial conditions described in panels a–e. Derived time-dependent concentration\nchanges of AA', AB', BA', and BB' from panels a and e are depicted in Fig. 3b and 3c.\nComment 3: “I would like to see at least one steady state for the three-module network\nthat is not at the edge of the triangle in Fig. 5d.”\nResponse: We appreciate the reviewer's constructive suggestion. In response, we have\nconducted additional experiments to better characterize the steady-state landscape of\nAMN-2, as shown in the new Fig. 5d and the new Supplementary Fig. 41.\nAs expected, these new investigations revealed two distinct steady states within the\ninterior of the triangular phase space, rather than being confined to its edges as\npreviously shown. The results further support the existence of a planar continuum\n(triangular phase space) of steady states in AMN-2. A discussion of these findings has\nbeen included in the revised manuscript (page 16).\nFig. 5｜Design of a three-module network, AMN-2, and its planarly distributed steady states.\na, Schematic representation of AMN-2, consisting of three autoregulatory DNA\npolymerization/nicking reaction modules. b, Time-dependent concentration changes of CC', CD',\nCE', DC', DD', and DE' in AMN-2 initiated by six different initial compositions, derived from data\nin Supplementary Fig. 37 using the corresponding calibration curves (Supplementary Fig. 39). c,\nCompositional fraction variations in modules DD'+CC'+DE', CD'+DC'+DE', and DC'+DD'+CE'\nat\nt = 0 and 12 h, derived from data shown in b. d, Planarly distributed steady states. Experimentally\ncharacterized steady states marked with pentagrams (S to S ) after 12 h, represented as the\n2-1 2-8\nrelative fractions of the three modules. The distribution exhibits a planar manifold, providing strong\nsupport for the emergence of high-order multistability in AMN-2.\nSupplementary Fig. 41. Dynamic behaviors of network AMN-2 in steady states S and\n2-7\nS . a,b, Time-dependent fluorescence changes of network AMN-2 in different steady states:\n2-8\nS , initiated with a composition of 40% DD'+CC'+DE', 30% CD'+DC'+DE', and 30%\n2-7\nDC'+DD'+CE' (a) and S , initiated with a composition of 20% DD'+CC'+DE', 60%\n2-8\nCD'+DC'+DE', and 20% DC'+DD'+CE' (b). c,d, Time-dependent concentration changes of CC',\nCD', CE', DC', DD', and DE' derived from fluorescence data shown in a and b by applying the\ncorresponding calibration curves (Supplementary Fig. 39).\nComment 4: “In all cases I see here, the modules (let's say AA' + BB') start out\nbalanced: there is an equal amount of AA' and BB'. Is this required for the system to\nwork? What happens when we create an imbalanced stoichiometry, keeping the total\nstoichiometry correct by applying an inverse perturbation on the AB'/BA' module?”\nResponse: We thank the reviewer for this insightful and technically important question.\nIndeed, in all cases presented in the manuscript, modules (e.g., AA'+BB') were initially\nprepared with equimolar concentrations of the respective DNA duplexes. However, this\nstoichiometric symmetry is not a strict requirement for the system's operation.\nTo evaluate the impact of imbalanced initial stoichiometries, we investigated a\nrepresentative case with 70% AA' and 30% BB'. The reviewer raises an important point\nregarding whether such an imbalance could be compensated by applying an inverse\nperturbation to a complementary module, such as AB'+BA', in an effort to preserve the\noverall component stoichiometry (i.e., total A, A', B, and B'). However, in practice,\nsuch a compensation is not feasible: increasing the concentrations of AB' and BA' to\nrestore B and B' levels inherently introduces excess A and A', thereby disturbing the\nglobal balance of the network. Consequently, a perfect restoration of overall\nstoichiometry via inverse perturbation is not achievable.\nTo examine the impact of this initial asymmetry, we performed additional\nexperiments on AMN-1 initialized with 70% AA' and 30% BB'. The system was\nanalyzed both in the absence and presence of DNA polymerase, nicking endonuclease,\nand dNTPs. The results have been included in the new Supplementary Fig. 17.\nThe results indicate that in the absence of DNA polymerase, nicking endonuclease,\nand dNTPs, the system gradually evolves toward a thermodynamic equilibrium. In\ncontrast, when the polymerization/nicking reactions are active in the presence of DNA\npolymerase, nicking endonuclease, and dNTPs, the network stabilizes near its initial\nimbalanced state. These results demonstrate that the emergence and maintenance of\nmultistability in our DNA replication network do not critically depend on perfect\nstoichiometric symmetry at initialization, underscoring the robustness of the system to\ncompositional perturbations.\nWe have discussed this stoichiometric symmetry issue in the revised manuscript\n(page 11).\nSupplementary Fig. 17. Dynamic behaviors of network AMN-1 with an initial composition\nof 70% AA'+30% BB' in the absence/presence of DNA polymerase, nicking endonuclease,\nand dNTPs. a,b, Time-dependent fluorescence changes of the network in the absence (a) and\npresence (b) of DNA polymerase, nicking endonuclease, and dNTPs. c,d, Time-dependent\nconcentration changes of AA', AB', BA', and BB' derived from fluorescence data shown in a\nand b by applying the corresponding calibration curves (Supplementary Fig. 4).\nComment 5: “Possibly related to that: the data shows that the two components in a\nmodule do not always have the same concentrations, sometimes showing quite a big\ndifference, see for instance Fig. 4e. What is the origin of this?”\nResponse: We thank the reviewer for this insightful observation. As noted in Fig. 4e,\nthe concentrations of the two DNA duplexes within a given module (e.g., AA' and BB'\nin module AA'+BB') are initially comparable. However, upon addition of an RNA\ntrigger, a marked divergence in their concentrations emerges. Below, we provide a\nmechanistic explanation, supported by thermodynamic analysis and NUPACK\nsimulations, to elucidate the origin of this asymmetry.\nIn AMN-1, each module consists of two DNA duplexes (e.g., AA' and BB'), which\ncoexist in dynamic equilibrium with their corresponding single-stranded DNA\ncomponents (A/A' and B/B'). We have used NUPACK to simulate the impact of RNA\ntriggers on the network composition in the absence of DNA polymerase, nicking\nendonuclease, and dNTPs. NUPACK simulation (new Supplementary Fig. 3a) shows\nthat, in the absence of RNA triggers, the network exists in a formation-dissociation\nequilibrium. Specifically, the concentrations of the DNA duplexes are approximately\nequal, and the concentrations of their associated single strands are approximately equal.\nUpon introduction of an RNA trigger (e.g., rT , new Supplementary Fig. 26b), two\nAA'\ninterrelated processes are selectively perturbed: (i) the reversible strand displacement\nbetween modules (e.g., AA'+BB' ⇌ BA'+AB'), and (ii) the intramodular equilibria\nbetween duplexes and their single-stranded components. RNA trigger rT binds\nAA'\npreferentially to its complementary sequences (either single strands A and A' or duplex\nAA'), forming a more stable hybrid complex (rT AA'). This preferential binding\nAA'\neffectively sequesters A and A', thereby shifting the module equilibrium to favor the\nformation of rT AA' through the dissociation of BA' and AB'. Concurrently, the\nAA'\ndissociation of BA' and AB' releases additional B and B' for the formation of BB'.\nThermodynamic simulations (new Supplementary Fig. 25a) reveal that the free\nenergy of rT AA' (ΔG ≈ −40.88 kcal·mol−1) is more negative (i.e., more favorable)\nAA'\nthan that of the BB' duplex (ΔG ≈ −16.41 kcal·mol−1). This difference leads to\npreferential stabilization and accumulation of rT AA', while BB' continues to\nAA'\nparticipate in dynamic equilibria with its corresponding single strands. As a result, a\nconcentration asymmetry between the duplexes gradually emerges.\nThis asymmetry is further amplified over time due to a memory effect inherent to\nAMN-1: although the RNA trigger is degraded by RNase H, it can be re-added across\nmultiple triggering cycles. This recursive reinforcement progressively biases the system\ntoward accumulation of specific species (AA'), which is clearly observed in Fig. 4e (\nleft\npanel). Importantly, this phenomenon is not unique to the rT /AA'+BB' system.\nAA'\nSimilar asymmetric behaviors are observed in response to other RNA triggers, such as\nrT , as shown in Fig. 4d and Fig. 4e (right panel).\nAB'\nTo clarify this phenomenon, we have revised the manuscript accordingly (page 14).\nSupplementary Fig. 3. NUPACK simulations of thermodynamically equilibrated\ncompositions in AMN-1 and AA'P -induced network reconfiguration, in the absence of\n1, confirming the coexistence of hybridized species (AA', BB', AB' and BA') with identical\nconcentrations, and unpaired single strands, A, A', B and B'. b, Compositional simulation of\nAA'P -induced network reconfiguration, where AA'P +BB' becomes dominant owing to the\n1 1\ninhibition of the forward reaction (AA'+BB' → BA'+AB').\nSupplementary Fig. 26. NUPACK simulations of trigger-induced state transitions in\nAMN-1 in the absence of DNA polymerase, nicking endonuclease, and dNTPs. a,\nCompositional simulation of rT -induced network reconfiguration, where BA'+rT AB'\nAB' AB'\nbecomes dominant owing to the inhibition of the inverse reaction (BA'+AB' → AA'+BB'). a,\nCompositional simulation of rT -induced network reconfiguration, where rT AA'+BB' are\nAA' AA'\ndominated due to the inhibition of the forward reaction (AA'+BB' → BA'+AB'). Together, these\nsimulations demonstrate that specific RNA triggers can selectively stabilize distinct low-energy\nnetwork configurations, thereby enabling programmable transitions between steady states.\nImportantly, the introduction of either trigger induces emergent asymmetries in the\nconcentration changes among DNA duplexes within each module. Compared to the unperturbed\nbaseline (Supplementary Fig. 3a), equilibrium simulations reveal the coexistence of both\nhybridized duplexes (AA', BB', AB', and BA') and their corresponding single-stranded\ncomponents (A, A', B, and B'). Upon trigger binding, single-stranded DNA species\ncomplementary to the RNA trigger engage in sequence-specific hybridization (e.g.,\nrT +A+A'→ rT AA'), resulting in the formation of more thermodynamically stable\nAA' AA'\ncomplexes. This interaction perturbs the native dissociation–association equilibrium of the\nduplexes (e.g., AA' ⇌ A+A'), leading to disproportionate compositional shifts among DNA\nspecies directly involved in trigger hybridization relative to those that are not.\nSupplementary Fig. 25. Trigger-induced state transition of network AMN-1 in the absence\nof DNA polymerase, nicking endonuclease, and dNTPs. a, Free energy landscape of AMN-\n1 network reconfiguration before and after treatment with rT , based on NUPACK simulations.\nAA'\nBefore triggering, the ΔG values of modules BA'+AB' and AA'+BB' are equal, suggesting the\nreversible module interconversion. However, the treatment with rT results in the formation\nAA'\nof rT AA' with more negative free energy, thus suppressing the forward reaction (AA'+BB'\nAA'\n→ BA'+AB'). b, Schematic mechanism of rT -triggered reconfiguration of the network under\nAA'\nnon-enzymatic conditions. Prior to triggering, the network resides in an equilibrated state\ncomprising modules AA'+BB' and AB'+BA' with comparable free energy, bidirectionally\nconverting through both antagonistic and agonistic regulatory interactions. The introduction of\nrT leads to the formation of rT AA' with more negative free energy, thus driving the\nAA' AA'\nreconfiguration of the network where AA'+BB' are dominated. Antagonistic and agonistic\nconnectivities are denoted by arrows and T-shaped lines, respectively.\nComment 6: “Many plots showing the development of the steady-states do not seem to\nbe truly steady and show drift of a few percent per hour even after 12 hours (Fig. 3d\nS1-1 for instance). Is this a question of waiting longer for true stability, are these effects\nfrom consumption of dNTPs, aging reagents, or is something else going on?”\nResponse: We thank the reviewer for the thoughtful observation regarding the drift in\nsome of the steady-state trajectories, particularly in Fig. 3d (trajectory S ). This raises\n1-1\nan important point about the robustness of steady states in AMN-1.\n(i) Quantitative characterization of steady-state drifts\nTo address this concern, we have conducted a detailed analysis of the system's long-\nterm behavior. Our results suggest that the network dynamics can be broadly divided\ninto two regimes. The initial phase (~0.5 h) is characterized by a short relaxation, during\nwhich polymerization/nicking reactions exert an increasing, suppressive influence on\nthe reversible displacement reaction. This is followed by a second prolonged phase,\nwhich corresponds to the operational steady state. In this regime,\npolymerization/nicking reactions have fully established kinetic dominance, effectively\nhalting further spontaneous interconversion between network modules. However, this\nregime is not entirely static—compositional changes occur at a markedly slower rate.\nTo quantitatively assess this drift, we calculated the drift rate as the change in the\nconcentration of each module per hour, normalized by the total concentration of all\nmodules. This calculation provides a practical metric to evaluate the long-term stability\nof the system. Under all initial conditions, AMN-1 reliably settles into composition-\nspecific steady states (S to S ) with small compositional drifts (revised Fig. 3t, see\n1-1 1-5\nComment 6). The largest drift was observed in S , with a maximum rate of ~0.5%·h−1\n1-1\nover a 12-hour period. It should be noted that comparable drift phenomena have been\nreported in conventional multistable biochemical systems (ref. 58—Nat. Commun. 10,\n4636 (2019), ref. 61—Nat. Chem. 11, 829–838 (2019), ref. 62—Nat. Chem. 14, 1224–\n1232 (2022), ref. 63—Proc. Natl. Acad. Sci. USA 109, E3212–E3220 (2012) and ref.\n64—Nat. Commun. 7, 13474 (2016)).\n(ii) Experimentally investigating the origin of the drift and strategies to improve\nlong-term robustness\nTo determine the origin of this drift, we have conducted a series of long-duration\nexperiments, focusing on several potential factors: enzyme inactivation (e.g., DNA\npolymerases and nicking enzymes), nonspecific side reactions, and consumption of\ndNTPs. Among these, our data strongly suggest that dNTP depletion is the primary\ndriver of the observed drift. As dNTPs serve as the energy source for polymerization,\ntheir gradual consumption naturally limits the system’s ability to maintain a far-from-\nequilibrium state over time.\nAt the original dNTP concentration of 0.1 mM, the drift becomes pronounced after\napproximately 36 h, reaching a drift rate of ~0.7%·h−1 during the final 12 hour-window\n(new Supplementary Fig. 13). In contrast, increasing dNTP concentration to 0.2 mM\nsignificantly mitigates the drift (~0.07%·h−1 by 48 h). These results indicates that fuel\ndepletion is the primary cause of the observed drift. Under conditions of adequate dNTP\nsupply, the AMN-1 network is capable of maintaining a stable steady-state composition\nwith minimal drift.\nTo further assess the energy dependence of steady-state stability, we systematically\nvaried the dNTP concentration within the range of 0–0.02 mM (new Supplementary\nFigs. 14 and 15). At suboptimal dNTP levels, multistability is compromised, and\ncompositional drift increases markedly. Conversely, higher dNTP concentrations\nenhance both the persistence and fidelity of the steady states.\nOverall, these results support the conclusion that gradual dNTP depletion is the\nprimary driver of drift in the system. Notably, increasing the dNTP concentration offers\na simple and effective strategy for improving the long-term stability and robustness of\nthe network.\nThe new experiments and analyses have been summarized and included in the\nrevised manuscript (page 11) and in the new Supplementary Section 5.\nFig. 3｜Arbitrary multistability of AMN-1. a, Schematic representation of AMN-1 without\nautoregulatory loop, in the absence of DNA polymerase, nicking endonuclease, and dNTPs. b,c,\nTime-dependent concentration changes of AA', BB', BA', and AB' in the network initiated by two\ndistinct compositions: 100% AA'+BB' (b) and 100% BA'+AB' (c). d, Computationally derived\nreversible displacement rates. e, Schematic representation of a modified AMN-1, where A' was\nreplaced by A', a variant without template domain P ', leading to a single autoregulatory loop of\ntwo initial compositions: 100% AA'+BB' (f) and 100% BA'+AB' (g). h,i,\nCorresponding\ncomputationally derived reversible displacement rates (h) and polymerization/nicking rates\nassociated with AB' (i). j, Schematic representation of another modified AMN-1 with a single\nautoregulatory loop of AA'+BB'. Here, B' was replaced by B', a variant without template domain\nP '; k,l, Time-dependent concentration profiles of AA', BB', BA', and AB' in the network initiated\n2\nfrom two initial compositions: 100% AA'+BB' (k) and 100% BA'+AB' (l). m,n, Computationally\nderived reversible displacement rates (m) and polymerization/nicking rates associated with AA' (n).\no−s, Time-dependent concentration changes of AA', BB', BA', and AB' in the AMN-1 with two\ncompeting autoregulatory loops under five initial compositions, where the proportion of AA'+BB'\nvaries from 100% (o), 75% (p), 50% (q), and 25% (r) to 0% (s), with corresponding variations in\nBA'+AB' from 0%, 25%, 50%, and 75% to 100%; The network converges to five different steady\nstates (S to S ). t, Drift rates of BA'+AB' by 12 h under these steady states, where the drift rate\n1-1 1-5\nis defined as the change in the concentration of each module per hour, normalized by the total\nconcentration of all modules. u, Linearly distributed steady states. Uniformly sampled theoretically\npredicted steady states of AMN-1 (216 circle dots) represented in terms of the fractions of AA'+BB'\nversus BA'+AB'. The distribution densely populates a linear manifold, providing strong theoretical\nsupport for the emergence of linear, continuous multistability in AMN-1. Pentagram markers (S\n1-1\nto S ) correspond to the experimentally characterized steady states after 12 h. v−x,\n1-7\nComputationally derived reversible displacement rates (v) and polymerization/nicking rates of AB'\n(w) and AA' (x). Time-dependent concentration profiles in b/c, f/g, k/l, and o−s were derived from\nSupplementary Figs. 5, 7, 8, and 12, respectively, using the respective calibration curves\n(Supplementary Fig. 4). Solid curves represent experimentally determined data, and dashed curves\nindicate computationally simulated results.\nSupplementary Fig. 13. Dynamic behaviors of network AMN-1 in steady state S at the\n1-1\ndNTP concentrations of 0.1 and 0.2 mM in 48 h. a,b, Time-dependent fluorescence changes\nof AMN-1 in S at different dNTP concentrations: 0.1 mM (a) and 0.2 mM (b). c and d, Time-\n1-1\ndependent concentration changes of AA', AB', BA', and BB' derived from the fluorescence data\nin a and b using corresponding calibration curves (Supplementary Fig. 4).\nSupplementary Fig. 14. Fluorescence kinetics of AMN-1 in steady state S at the dNTP\n1-1\nconcentrations of 0 − 0.020 mM. a−e, Time-dependent fluorescence changes related to AMN-\n1 in S at different dNTP concentrations: 0 mM (a); 0.001 mM (b); 0.005 mM (c); 0.010 mM\n1-1\n(d); 0.020 mM (e).\nSupplementary Fig. 15. Time-dependent composition changes of AMN-1 in steady state\nS at the dNTP concentrations of 0 − 0.020 mM. a−e, Time-dependent concentration\n1-1\nchanges of AA', AB', BA', and BB' derived from fluorescence data shown in Supplementary\nFig. 14 by applying the respective calibration curves (Supplementary Fig. 4): AMN-1 in S at\n1-1\ndifferent dNTP concentrations: 0 mM (a); 0.001 mM (b); 0.005 mM (c); 0.010 mM (d); 0.020\nmM (e).\nComment 7: “The authors use a kinetic model to predict the behavior of the system,\nand mention it \"aligns closely\" with the experimental system. However, in several plots\nit is clear that the model deviates and does not fit even qualitatively (Fig. 3c-v, Fig. 3d-\nS1-5, Sup. Fig. 9c&d, Sup. Fig. 8d, Sup. Fig. 12a, Sup. Fig. 17f). I realize the deviations\nare not too big, but can the authors formulate a reason why this is the case?”\nResponse: We thank the reviewer for pointing out specific instances where the kinetic\nmodel deviates from the experimental data. As correctly noted, although the model\nmostly captures the qualitative trends and key dynamics of the system, we acknowledge\nthat certain discrepancies remain. We have elaborated on the possible sources of these\ndeviations:\nTheoretical considerations and sources of model–experiment deviations\nIt is worth noting that, in certain experiments, observable deviations exist between\ntheoretical simulations and experimental results. These discrepancies primarily arise\nfrom intentional model simplifications that aim to strike a balance between biochemical\nrealism and predictive tractability. Below, we outline the major contributing factors and\nanalyze their potential impact in detail:\n(a) Simplification of the reversible displacement reaction kinetics: To reduce the\nmodel’s complexity and limit the number of free parameters, we assumed that the rate\nconstants, k and k , governing the reversible displacement reaction between modules\n1 2\nAA'+BB' and BA'+AB' (Supplementary Eq. 5), are identical to those used in other\nrelated reactions (Supplementary Eqs. 6−8). While the interacting domains in these\nreactions are identical in sequence, they differ in other structural domains (e.g., AA'-P\n1\nvs. AA'), which could affect hybridization dynamics and displacement efficiency.\nFurthermore, the model does not account for the formation-dissociation equilibria\nbetween duplexes and their constituent strands, an omission that might influence the\neffective reaction rates. These simplifications were necessary to preserve\ngeneralizability and analytical tractability but inevitably introduce sources of error that\ncontribute to the observed model-data deviations.\n(b) Assumption of constant enzymatic activity and linear kinetics: In the kinetic\nmodel, we assumed DNA polymerase and nicking enzyme activity as constant and\napproximated polymerization and nicking as first-order reactions, depending solely on\nenzyme concentration. In actual experimental conditions, enzyme activity may decline\nover time due to degradation, inactivation, or substrate depletion. Moreover, enzyme\nkinetics may exhibit nonlinear behaviors, which our model does not capture. While\nincorporating these effects would increase quantitative accuracy, it would also\nsignificantly increase model complexity.\n(c) Exclusion of dNTP consumption dynamics: In the kinetic model, we did not\ninclude dNTP depletion, but assumed that the polymerization process was unaffected\nby the consumption of dNTPs. However, as demonstrated in the new Supplementary\nFigs. 13−15 (see Comment 6), the concentration of dNTPs can significantly influence\nthe polymerization kinetics, thus impacting overall network dynamics. Incorporating\nthis factor could improve the model's accuracy but would also substantially increase its\ncomplexity and require additional experimental quantification.\nOverall, the discrepancies between the model and experimental data primarily stem\nfrom deliberate simplifications intended to balance biochemical realism with analytical\nand conceptual clarity. Although the current model does not capture every detail, it\nsuccessfully reproduces the essential multistable behavior and provides a robust\nmechanistic framework for understanding and engineering synthetic DNA replication\nnetworks.”\nThese clarifications have been included in Supplementary Section 4. In light of the\nreviewer’s suggestion, we have replaced the phrase “aligns closely” with “is generally\nconsistent with” throughout the text to more accurately reflect the level of agreement\nbetween the model and the experimental data.\nComment 8: “Is there a reason why the rTaa/rTab & RNase H reactions are not\ninvestigated with the kinetic model? The relevant rate constants seem to be determined\nin Sup. Fig. 20.”\nResponse: We sincerely thank the reviewer for this thoughtful and constructive\ncomment. In our original kinetic modeling, which comprised 10 reactions, we\nintentionally focused on capturing the core design principles and fundamental\nmechanisms underlying multistability in the DNA replication network. This simplified\nmodel enabled us to systematically analyze the network’s intrinsic dynamics without\novercomplicating the parameter space.\nIncorporating trigger-induced transitions, such as those involving rT and rT ,\nAA' AB'\nwould have required a substantial expansion of the model—adding approximately 13\nadditional reactions to account for the complex interactions between the trigger species\nand the network components. For instance, modeling rT -mediated activation requires\nAA'\nintroducing species, such as rT AA', rT AA'P , and their downstream intermediates\nAA' AA' 1\nto represent the reversible strand displacement, polymerization/nicking cycles, and\nRNase H-mediated degradation. To maintain tractability, we chose to exclude these\ntrigger dynamics from the initial model.\nIn response to the reviewer’s suggestion, we have now extended our kinetic model\nto explicitly include the trigger-induced transitions between steady states. Specifically,\nwe have incorporated the relevant reactions for both rT and rT , as detailed in\nAB' AA'\nSupplementary Eqs. 5−14 and 41−53 for rT and Supplementary Eqs. 5−14 and 68−80\nAB'\nfor rT . Simulations were performed using MATLAB R2019a with rate constants (k\nAA' 1\n− k , k , k , k and k ) listed in Supplementary Table 2.\n8 13 14 18 19\nThe extended model was used to simulate the system’s transient response to trigger\nstrand addition, following the experimental trajectories shown in Fig. 4b, 4d, and 4e.\nThe fitted kinetic profiles are presented in the newly added Supplementary Figs. 29 and\n32. Furthermore, the newly derived rate constants (k – k and k – k ), specific to\n15 17 20 22\nthe trigger-modulated reactions, were extracted through model fitting and are now\nreported in Supplementary Table 2.\nThe expanded model and computational simulation procedures have been included\ninclued\nin Supplementary Section 4 and mentioned in the revised manuscript (pages 13 and 14).\nSupplementary Fig. 29. Computationally simulated orthogonal state transitions of AMN-\n1 corresponding to experimental results in Fig. 4b and 4d. a,b, Computationally simulated\ntime-dependent concentration profiles of AA', BB', AB', and BA' during sequential treatment\nwith different RNA triggers and RNase H: transition from state S upon the sequential\n1-1\ntreatment with RNA trigger rT and RNase H (a) and transition from state S following the\nAB' 1-5\nsequential addition of RNA trigger rT and RNase H (b). These simulations are generally\nAA'\nconsistent with the experimental results, demonstrating the orthogonality and programmability\nof state transitions in AMN-1.\nSupplementary Fig. 32. Computationally simulated continuous state transitions of AMN-\n1 corresponding to the experimental results in Fig. 4e. a,b, Computationally simulated, time-\ndependent concentration profiles of AA', BB', AB', and BA' during a multi-step protocol\ninvolving initial treatment with RNase H and a specific RNA trigger, followed by sequential\naddition of the designated RNA triggers: transitions initiated from state S upon sequential\n1-1\naddition of rT (a); transitions initiated from state S upon sequential addition of rT (b).\nAB' 1-5 AA'\nThese simulations highlight the capacity of AMN-1 to undergo continuous, trigger-guided\ntransitions, demonstrating the continuity and programmability of state control within the DNA\nreplication network.\nComment 9: “I wonder about the robustness of the system. The authors convincingly\nshow that a steady state is found at or around the initial composition of the system, but\nhow robust is this when concentrations of the various catalytic components, dNTPs,\netc., are varied/consumed? A small discussion in the SI would be appreciated. I realize\nsome data is included (Sup. Fig. 13-18), but some written discussion would be helpful.”\nResponse: We thank the reviewer for this thoughtful and important comment\nregarding the robustness of the system under varying biochemical conditions. In\nresponse, we have expanded our analysis and discussion in the revised Supplementary\nInformation. Specifically, we investigated the effects of varying concentrations of the\ncatalytic components, DNA polymerase and nicking endonucleases, as presented in the\noriginal Supplementary Figs. 13–18 (Supplementary Figs. 18–23 in the revised\nSupplementary Information). Additionally, we examined the influence of the fuel\nsupply, particularly dNTP concentrations, on system stability and multistable outcomes.\nThe changes are summarized as follows:\nRobustness of multistability under varying reaction conditions\nThe multistable behavior of AMN-1 demonstrates notable robustness across a range\nof enzymatic and fuel (dNTP) concentrations. The key findings regarding the system’s\ntolerance to variations in reaction conditions are summarized below:\n(i) Robustness to enzymatic concentration variations: To evaluate the system's\ntolerance to fluctuations in enzymatic input, we independently varied the\nconcentrations of DNA polymerase and nicking enzyme while keeping all other\nreaction parameters constant. When the nicking enzyme concentration was fixed at 20\nU·mL−1, AMN-1 maintains multistability over a concentration range of 20–35 U·mL−1\nfor the DNA polymerase (Supplementary Figs. 18−20). Conversely, with the DNA\npolymerase concentration kept constant at 30 U·mL−1, the system retains multistable\ndynamics across a concentration range of 10–25 U·mL−1 for the nicking enzyme\n(Supplementary Figs. 21−23). These results indicate that AMN-1 exhibits a robust\nmultistable response to moderate fluctuations in catalytic activity, highlighting its\ntolerance to enzymatic variability.\n(ii) Robustness to dNTP concentration and depletion: We next investigated how\nvariations in dNTP concentration and consumption dynamics influence the long-term\nstability of the system. Experiments conducted with different initial dNTP\nconcentrations (new Supplementary Figs. 13–15, see Comment 6) reveal that dNTP\ndepletion plays a central role in driving compositional drifts, ultimately compromising\nmultistability over time.\nUnder original operating conditions (0.1 mM dNTPs), AMN-1 reliably settles into\ncomposition-specific steady states (S to S ), each determined by its initial\n1-1 1-5\nconfiguration, with minimal compositional drifts (revised Fig. 3t, see Comment 6). The\nlargest drift was observed in S , with a maximum rate of ~0.5%·h−1 over 12 h. To\n1-1\nprobe the source of this drift and evaluate stability over longer timescales, we conducted\nextended experiments under S while varying dNTP concentrations (new\n1-1\nSupplementary Fig. 13, see Comment 6). At a limiting dNTP concentration of 0.1 mM,\nthe drift becomes pronounced after approximately 36 h, reaching a drift rate of\n~0.7%·h−1 during the final 12-hour window. In contrast, increasing dNTP concentration\nto 0.2 mM significantly mitigates the drift, reducing the rate to ~0.07%·h−1 over 48 h.\nThese results indicate that dNTP depletion is the principal factor limiting long-term\ncompositional stability.\nTo further assess the energy dependence of steady-state stability, we tested a broader\nrange of suboptimal dNTP concentrations (0–0.02 mM, new Supplementary Figs. 14\nand 15). At suboptimal dNTP levels, multistability is compromised, and compositional\ndrift increases markedly. Conversely, higher dNTP concentrations enhance both the\npersistence and fidelity of the steady states. These findings underscore the essential role\nof energy input in sustaining far-from-equilibrium dynamics and maintaining steady\nstates.\nOverall, AMN-1 displays a high degree of robustness to variations in both enzymatic\nand fuel concentrations, maintaining multistable behavior across a broad operational\nwindow. Nevertheless, the system's stability is ultimately constrained by resource\ndepletion, particularly under energy-limited conditions.\nTo address the reviewer’s request for additional discussion, we have added a new\nSupplementary Section 5 (\"Robustness of Multistability under Varying Reaction\nConditions\") in the Supplementary Information and have summarized the key\nconclusions of these findings in the main text (page 11).\nComment 10: “Fig. 3d, bottom row, second from the left as well as Fig. 5d (and their\nequivalent plots in Fig. 1a&b): These plots do not add anything and suggest a higher\ndimensional system than is actually presented. For Fig. 3, it is really not clear what the\nmessage of the plot is supposed to be anyway. In Fig. 5, a more appropriate way of\nshowing this would be a ternary plot, which correctly shows 2 orthogonal dimensions\ninstead of 3 - which the authors do actually use in Fig. 1b.”\nResponse: We sincerely thank the reviewer for the insightful comment regarding the\nclarity and dimensional interpretation in Figs. 1a, 1b, 3d (bottom row, second from the\nleft in the original manuscript), and 5d. We fully agree that the original 3D\nrepresentations may have inadvertently conveyed an impression of higher\ndimensionality than is present in the system, potentially leading to confusion about the\nnature of the compositional space.\nIn response, we have implemented several substantial revisions to the relevant\nfigures, captions, and corresponding text to enhance clarity, aimed at enhancing clarity,\ninterpretability, and dimensional accuracy. The changes are summarized as follows:\n(i) Clarification and revision of Figs. 1 and 3 (new Fig. 3u)\nFig. 1: We have revised Fig. 1 and its caption to more clearly distinguish between\ndiscrete and continuous multistability. Specifically, the right panel of Fig. 1a now more\nexplicitly represents the discrete attractor states of conventional multistable networks,\nwhile Fig. 1b illustrates the ability of AMNs to support a continuous distribution of\nsteady states along linear or planar manifolds in compositional space. This visual\ncontrast is intended to reinforce the conceptual distinction between conventional and\narbitrary multistability.\nFig. 3u (previously Fig. 3d, bottom row, second from the left): This panel has\nbeen significantly revised to more accurately convey the continuous nature of the\nsteady-state distributions in AMN-1. In the updated figure (See Comment 6), we plot\nthe steady-state compositions of AA'+BB' versus BA'+AB' after 12 hours of reaction,\nas determined experimentally. The resulting data form a linear distribution, consistent\nwith the emergence of continuous multistability.\nTo further characterize this behavior, we have performed Linear Stability Analysis\n(LSA) based on the fitted kinetic model to map the full accessible steady-state space of\nthe network. Drawing inspiration from computational neuroscience (ref. 69—Nat.\nNeurosci. 5, 775–782 (2002) and ref. 70—Curr. Opin. Neurobiol. 13, 204–211 (2003)),\nwhere both linear and nonlinear network models have been proposed to explain\ncontinuous attractors and graded memory states, we note that ideal continuous attractors,\ncomprising an infinite number of exact steady states, are theoretically achievable in\nlinear systems. In contrast, biologically realistic nonlinear systems more commonly\ngive rise to approximate continuous attractors. Although they are not mathematically\nexact, they are typically robust and functionally meaningful within a biological context.\nGuided by this conceptual framework, we adopted a criterion for steady state\nidentification where the net concentration flux of each molecular species is constrained\nto be less than 1×10−10 nM·h−1 in LSA, consistent with established methodologies in\nnonlinear dynamical systems modeling. From this set, 216 representative steady states\nwere uniformly sampled and are plotted in the new Fig. 3u (as circular markers, see\nComment 6). The resulting distribution forms a densely populated linear manifold,\nproviding strong theoretical support for the presence of continuous multistability along\na one-dimensional subspace.\n(ii) Replacement of 3D schematic plots in Figs. 1b and 5d\nFollowing the reviewer’s suggestion, we have replaced the original 3D schematic\ndiagrams in Figs. 1b and 5d with 2D ternary plots (for the updated Fig. 5d, see Comment\n3). These revised plots more accurately reflect the intrinsic two-dimensionality of the\ncompositional space in AMN-2 and eliminate any misleading suggestion of a three-\ndimensional distribution of steady states. The revised figure captions now clearly state\nthat the steady states lie on a two-dimensional manifold within the ternary\ncompositional space, thereby reinforcing the correct interpretation of the system’s\ndimensionality.\nA detailed description of the LSA has been provided in Supplementary Section 4 and\nsummarized in the main text (page 12).\nFig. 1｜Comparison of reported multistability networks and synthetic arbitrary multistability\nnetworks developed in this work. a, Schematic of previously reported bistable and tristable\nnetworks with discrete steady states. Typically, in a bistable network, modules A and A mutually\n1 2\ninhibit one another through the generated inhibitory signals, leading to two discrete, point-like\nsteady states (I and II) in the compositional space. Similarly, in a tristable network, signal-mediated\nmutual inhibition among modules A , A , and A results in three discrete steady states (III, IV, and\n3 4 5\nV). b, Schematic representation of AMNs demonstrating arbitrary multistability. In AMN-1, far-\nfrom-equilibrium competition between two autoregulatory modules leads to stabilization at any\npoint along a continuous linear manifold in the compositional space, rather than discrete states. In\nAMN-2, a three-module design enables stabilization across a continuous two-dimensional planar\nmanifold in the compositional space.\nComment 11: “Miscellaneous minor points:\n- In my version of the SI, I see black lines around figure edges.\n- Supplementary Fig. 29b has hours displayed far away from the plots they are\nsupposed to represent.\n- If 3D plots are used, one should label their axes; most plots are missing the initial\ncompositions axis (e.g., something like \"fraction AA' + BB'\" in Fig. 3).\n- In Fig. 3, the usage of both a, b, c and i, ii, iii as figure labels independent of each\nother at the same time is confusing. In general, it is very hard to refer to specific panels:\nFig. 3d has 10 (!) unlabeled sub-panels; these should be numbered or labeled somehow.\nFigure 3 is the main offender, but it is worth going over all figures again.”\nResponse: We thank the reviewer for their careful and detailed assessment of both the\nmain manuscript and Supplementary Information. We appreciate the specific feedback\nregarding figure formatting, clarity, and labeling, and have now addressed all points\nraised as follows:\n(i) Black lines around figures in the Supplementary Information: All supplementary\nfigures have now been re-prepared with standardized settings to eliminate unintended\nblack borders and ensure uniform presentation across all panels.\n(ii) Misalignment of time labels in Supplementary Fig. 29b: We have corrected the\nplacement of the hour markers in the original Supplementary Fig. 29b to align them\naccurately with the corresponding data curves (see updated Supplementary Fig. 43). In\naddition, we identified a similar issue in Supplementary Fig. 30a, which has also been\nrevised accordingly (see updated Supplementary Fig. 44).\n(iii) Missing axis labels in 3D plots: We have added the missing axes, labeled as\ndifferent systems, to all relevant 3D plots (see revised Fig. 3d, 3h, 3i, 3m, 3n, 3v-3x and\nSupplementary Fig. 11e and 11f).\n(iv) Inconsistent and confusing panel labeling in Fig. 3: We have revised the labeling\nscheme of all subpanels in Fig. 3. All panels now follow a consistent alphabetical order\n(a, b, c, ...) throughout the figure, avoiding the previous mixture of Roman numerals\nand letters.\n(v) General figure review for consistency: In light of the reviewer’s comment, we\nhave thoroughly reviewed all figures in both the main text and Supplementary\nInformation. We have standardized labeling conventions, ensured that all axes and units\nare clearly defined, and refined figure layouts to enhance visual clarity and uniformity\nacross the manuscript.\nWe are grateful to the reviewer for identifying these important issues. We believe the\nimplemented revisions significantly improve the readability and professional\npresentation of our work.\nReviewer #2:\nComment 12: “In this paper, the authors present some observations about a network\nbased on polymerization nicking reaction on two templates competing for a single\ntrigger (where the stoichiometry is balanced by introduction of a non-extendable\npseudo-trigger). The polymerization/nickase cycles slow down the trigger exchange\nreaction, but does not abolishes it on longtime scales, a phenomenon that the author\nwrongly interpret as a new form of multistability. However, continuous multistablity\ndoes not exist mathematically. Slow relaxation does.”\nResponse: We thank the reviewer for raising these insightful and fundamental points\nregarding the term “arbitrary multistability within a continuous spectrum” from both\ntheoretical and long-timescale observational perspectives. We address the comment in\ntwo parts, concerning: (i) the conceptual legitimacy of “continuous multistability”\nand\n(ii) the distinction between true multistability and slow relaxation dynamics.\n(i) On the conceptual validity of “continuous multistability”\nWe agree that “continuous multistability” is unconventional in the context of\nsynthetic biochemical networks, where multistability is typically defined as a system\npossessing a finite number of well-separated, discrete steady states (e.g., bistability or\ntristability, ref. 58—Nat. Commun. 10, 4636 (2019), ref. 61—Nat. Chem. 11, 829–838\n(2019) and ref. 62—Nat. Chem. 14, 1224–1232 (2022)). Nevertheless, closely\nanalogous phenomena exhibiting effectively continuous ranges of steady states have\nbeen rigorously described in other disciplines.\nSoft-robotics and mechanical metamaterials: As stated in ref. 71—Proc. Natl.\nAcad. Sci. USA 117, 5217–5221 (2020), mechanical systems comprising arrays of\nbistable mechanical elements (e.g., clamped-clamped curved beams) can produce\nmultistable behaviors where a single control input (e.g., pressure) yields a functionally\ncontinuous range of output configurations. Although the number of possible stable\nstates is mathematically finite and bounded by the number of bistable elements, the\nsystem can be approximated as a continuous manifold of configurations. As the number\nof coupled units increases, the number of distinct, accessible configurations grows\ncombinatorially, approaching a quasi-continuum in practice.\nNeuroscience and continuous attractor networks: Neural systems, such as the\noculomotor integrator, serve as another well-established example. These systems\nmaintain eye position across a continuous angular range and are modeled as continuous\nattractor networks, capable of storing analog signals as a continuum of persistent,\nsteady states (ref. 68—Proc. Natl. Acad. Sci. USA 93, 13339–13344 (1996), ref. 69—\nNat. Neurosci. 5, 775–782 (2002) and ref. 70—Curr. Opin. Neurobiol. 13, 204–211 (2003)).\nNotably, theoretical studies based on these systems strongly support the\nexistence of the linear, continuous attractors. Specifically, both linear and nonlinear\nnetwork models have been proposed to explain continuous attractors and graded\nmemory states. In linear models, a theoretically infinite number of steady states can be\nachieved. In contrast, biologically realistic nonlinear systems more commonly give rise\nto approximate continuous attractors. Although they are not mathematically exact, they\nare typically robust and functionally meaningful within a biological context. Therefore,\nthese models offer a theoretical foundation for continuous multistability in biological\nsystems.\nIn both cases, the systems operate far from thermodynamic equilibrium and exhibit\nbehaviors consistent with continuous multistability, thereby reinforcing the conceptual\nlegitimacy of this phenomenon.\nIn analogy, in our DNA-based architectures, multiple modules are simultaneously\nstabilized through far-from-equilibrium competition, constrained by global mass\nbalance, and maintained out of equilibrium via enzymatic fueling. This dynamic\nconfiguration gives rise to what we describe as arbitrary multistability within a\ncontinuous spectrum—a regime in which the system robustly stabilizes a broad,\ntunable range of compositionally distinct steady states. Each of these states corresponds\nto a unique distribution of module concentrations, and together they form an effectively\ncontinuous landscape of accessible configurations. Notably, this behavior is not merely\ntheoretical; it has been experimentally realized through the far-from-equilibrium\ndynamics of the coupled DNA replication network.\nWe have revised the manuscript (Introduction, page 4) to clarify the conceptual\nfoundation of this terminology, referencing established frameworks from neuroscience\nand materials science.\n(ii) On distinguishing slow relaxation from true multistability\nWe appreciate the reviewer’s concern that the observed persistence of system states\nmight arise from slow kinetic relaxation rather than genuine multistability. To address\nthis, we performed long-timescale experiments under varying fuel concentrations to\ndetermine whether the observed behavior arises from slow drift or active stabilization.\nUnder all initial conditions, AMN-1 reliably settles into composition-specific steady\nstates (S to S ) with small compositional drift (revised Fig. 3t, see Comment 6). The\n1-1 1-5\nlargest drift is observed in S , with a maximum rate of ~0.5%·h−1 over a 12-hour\n1-1\nperiod. Comparable drift phenomena have been reported in conventional multistable\nbiochemical systems (ref. 58—Nat. Commun. 10, 4636 (2019), ref. 61—Nat. Chem. 11,\n829–838 (2019), ref. 62—Nat. Chem. 14, 1224–1232 (2022), ref. 63—Proc. Natl. Acad.\nSci. USA 109, E3212–E3220 (2012) and ref. 64—Nat. Commun. 7, 13474 (2016)). To\ninvestigate the source of this drift and assess long-term stability, we conducted extended\nexperiments under S while varying dNTP concentrations (new Supplementary Fig.\n1-1\n13, see Comment 6). At the original dNTP concentration of 0.1 mM, the drift becomes\npronounced after approximately 36 h, reaching a drift rate of ~0.7%·h−1 during the final\n12-hour window. In contrast, increasing dNTP concentration to 0.2 mM significantly\nmitigates the drift (~0.07%·h−1 by 48 h), indicating that fuel depletion is the primary\ncause of this drift. Under conditions of adequate dNTP supply, the AMN-1 network is\ncapable of maintaining a stable, steady-state composition with minimal drift, implying\nthe robustness of AMN-1’s multistability. To further assess the energy dependence of\nsteady-state stability, we tested a broader dNTP range (0–0.02 mM, new Supplementary\nFigs. 14 and 15, see Comment 6). At suboptimal dNTP levels, multistability is\ncompromised, and compositional drift increases markedly. Conversely, higher dNTP\nconcentrations enhance both the persistence and fidelity of the steady states.\nThese findings collectively indicate that the multistable behavior is an emergent\nproperty of the system’s nonequilibrium dynamics and is actively maintained by\nenzymatic fuel consumption. In the presence of sufficient dNTPs, interconversion\nbetween module states is strongly suppressed, and the system remains kinetically\n“locked” in its selected state.\nIn summary, we have clarified the conceptual basis of arbitrary multistability within\na continuous spectrum by drawing upon established theoretical frameworks from\nneuroscience and mechanical systems. Furthermore, we have presented new\nexperimental data demonstrating that AMN-1's long-term state stability is an active,\nfuel-dependent process and not the result of slow passive relaxation.\nThese discussions have been included in the revised manuscript (pages 4 and 11).\nThe experimental results have also been included in the new Supplementary Section 5.\nComment 13: “In addition it is inappropriate to discuss steady states (and hence their\nstability) in a system where some of the variables (here compounds) can only increase,\nwhich is the case here for strand P, which is continuously produced and accumulates.”\nResponse: We thank the reviewer for raising this important issue, which emphasizes a\nfundamental consideration about the inappropriateness of discussing steady states due\nto the continuous production and accumulation of the single-stranded DNA species P .\nx\nTo address this concern, we have clarified the role of P within the AMN architecture\nx\nand elaborated on its functional decoupling from the regulatory loop that governs\nmultistability. This distinction is critical in understanding how our system maintains\nwell-defined steady states despite the continuous accumulation of P .\nx\nBiochemical multistable systems are inherently open and operate far from\nthermodynamic equilibrium. They require continuous energy input to sustain their\ndynamics and maintain distinct attractor states. In conventional multistable\narchitectures, such as networks based on genelet or PEN\n(polymerase/exonuclease/nickase) toolkits (ref. 61—Nat. Chem. 11, 829–838 (2019),\nref. 62—Nat. Chem. 14, 1224–1232 (2022), ref. 63—Proc. Natl. Acad. Sci. USA 109,\nE3212–E3220 (2012) and ref. 64—Nat. Commun. 7, 13474 (2016)), this energy input\nis provided by nucleotide hydrolysis, which drives the production of RNA or DNA\nspecies. These strands serve as regulatory species, mediating mutual activation or\ninhibition to shape the system's multistability. Importantly, their concentrations are\nmanaged through synthesis and degradation cycles, creating a multistable far-from-\nequilibrium system. During degradation, inert byproducts (RNA or DNA fragments)\nare continuously generated and accumulate, but they remain functionally decoupled\nfrom the core regulatory loops and thus do not interfere with the system’s multistable\nbehavior.\nOur AMN architecture adheres to the same core principles of synthesis-degradation,\nbut it functions through a distinct mechanism. Specifically, it also relies on dNTP\nhydrolysis for energy input. However, the regulatory species in the AMN system are\nthe double-stranded polymerized complexes AA'P and AB'P , rather than the single-\n1 2\nstranded DNA species P . These polymerized complexes mediate competitive\nx\ninteractions between network modules, thereby generating and maintaining\nmultistability. Their concentrations are dynamically regulated through a reversible\ncycle of template-guided polymerization and enzymatic nicking, which regenerates the\nprecursor duplexes (AA' and AB') and releases P as a byproduct. In other words, it is\nx\nthe dynamics of AA'P₁ and AB'P₂—not the accumulation of P —that govern the\nx\nmultistability of the system. As a result, the network can settle into well-defined steady\nstates, even as P continues to accumulate over time.\nx\nFurthermore, unlike typical inert byproducts, P maintains a well-defined nucleotide\nx\nsequence and can be intentionally programmed as a signaling intermediate. This feature\noffers a key design advantage, allowing P to transmit information to downstream\nx\nmodules without disrupting the internal regulatory dynamics and robustness of the\nmultistable core. In contrast, traditional architectures often rely on co-opting regulatory\nstrands for signal export, which could destabilize the multistable behavior or\ncompromise network modularity.\nWe have emphasized the distinction between regulatory and non-regulatory species,\nand between inert byproducts and programmable signal carriers in the revised main text\n(page 6). We have also elaborated on the unique role of P as a functionally isolated,\nx\nexportable signal carrier that enables the scalable and hierarchical design of synthetic\nmolecular circuits without compromising core multistability in the revised main text\n(page 18).\nComment 14: “Furthermore mathematical support for multistability may not come\nfrom kinetic ODE fitting but from Linear Stability Analysis (LSA), which I encourage\nthe author to try (once they have corrected their kinetic model). If the authors prefer\nan experimental test, they can also attempt to adjust the balance between\npolymerization and nicking activities, upon which they will probably find that the\nsystem does relax to a single AA’/AB’ ratio, which a time scale that decrease with\npol/nick ratio, but this does not mean the system is multistable (the situation presented\nin Fig 3a being the limit pol/nick->0).”\nResponse: We sincerely thank the reviewer for this insightful and constructive\ncomment. We fully agree that Linear Stability Analysis (LSA) provides a powerful\ntheoretical approach to assess the stability of steady states. In response, we have\nclarified the relationship between kinetic ODE fitting and LSA, applied LSA to map\nthe full accessible steady-state space of the AMN network and explained how the\nbalance between polymerization and nicking activities influences the system’s\ndynamics, as follows:\n(i) Complementarity of kinetic ODE fitting and LSA:\nOur modeling framework follows a standard two-stage process. First, kinetic ODE\nmodeling is employed to capture the nonlinear dynamics of the system and to fit\nexperimentally determined kinetics, thereby inferring unknown rate constants. Once the\nmodel parameters are reliably constrained, LSA is then performed to assess the stability\nof candidate steady states.\nIn this sense, ODE fitting and LSA are not alternatives but complementary steps: the\nformer provides a quantitatively grounded and mechanistically faithful model, while\nthe latter enables the characterization of the system’s stability landscape. In our study,\nthe kinetic model was first calibrated using experimental data to ensure fidelity to\nobserved system dynamics. We then applied LSA based on this model to investigate\nwhether the network could support multiple steady states.\n(ii) LSA reveals approximate continuous multistability\nIn response to the reviewer’s suggestion, we have performed LSA based on the fitted\nmodel to explore the full accessible steady-state space of AMN-1. Drawing inspiration\nfrom computational neuroscience (ref. 69—Nat. Neurosci. 5, 775–782 (2002) and ref.\n70—Curr. Opin. Neurobiol. 13, 204–211 (2003)), where both linear and nonlinear\nnetwork models have been proposed to explain continuous attractors and graded\nmemory states, we note that ideal continuous attractors, comprising an infinite number\nof exact steady states, are theoretically achievable in linear systems. In contrast,\nbiologically realistic nonlinear systems more commonly give rise to approximate\ncontinuous attractors. Although they are not mathematically exact, they are typically\nrobust and functionally meaningful within a biological context. Guided by this\nconceptual framework, we adopted a criterion for steady state identification where the\nnet concentration flux of each molecular species is constrained to be less than 1×10−10\nnM·h−1 in LSA, consistent with established methodologies in nonlinear dynamical\nsystems modeling. In this case, we identified thousands of stable steady states. From\nthese, 216 representative states were uniformly sampled across the full accessible\ncompositional range and are depicted in the new Fig. 3u (circle dots, see Comment 6).\nThe distribution densely populates a linear manifold, providing strong theoretical\nsupport for the emergence of linear, continuous multistability in AMN-1. This\norganization reflects the emergence of an approximate continuum of attractors,\nreinforcing our interpretation of AMN-1 as a system that exhibits biologically plausible\narbitrary multistability. As a control, we conducted an identical LSA in an enzyme-free\nvariant of the system, using the same net concentration flux applied to the 216\nrepresentative steady states. Unlike the multistable behavior exhibited by AMN-1, the\nenzyme-free system consistently yields a single fixed-point equilibrium. This contrast\nhighlights the essential role of the competing polymerization/nicking reactions in\ngenerating and maintaining the multistability observed in AMN-1.\n(iii) Experimental validation through enzyme ratio modulation:\nRegarding the reviewer’s suggestion for experimental validation by varying the\npolymerase-to-nickase (pol/nick) ratio, we agree this is a valuable means to probe the\nsystem’s dynamical regimes. As the reviewer notes, in the limit of vanishing\npolymerase activity (i.e., pol/nick → 0), the network is expected to relax to a unique\nAA'/AB' ratio, consistent with the monostable regime as shown in Fig. 3b and 3c.\nHowever, this does not contradict the presence of multistability under more balanced\nenzymatic conditions. In fact, it reflects a well-established property of multistable\nbiochemical systems: multistability typically arises only within a finite, bounded region\nof parameter space. When operating outside this region—for example, under extreme\nreaction parameters—networks often revert to monostable behaviors (ref. 58—Nat.\nCommun. 10, 4636 (2019) and ref. 63—Proc. Natl. Acad. Sci. USA 109, E3212–E3220\n(2012)). Our original study included a series of experiments in which the pol/nick\nratio was systematically varied. As shown in Supplementary Figs. 18–23, AMN-1\npreserves its multistable behavior over a broad operational window, confirming the\nsystem’s robustness.\nWe have incorporated a new section in the Supplementary Information (Section 5)\nto detail how enzymatic concentrations affect system behavior. The corresponding\nLSA-based analysis and experimental results have been included in the revised main\ntext (page 12) and Supplementary Section 4.\nReviewer #3:\nComment 15: “The manuscript by Zhong et al. introduces a novel strategy for\nconstructing multistable biochemical networks using reversible DNA strand\ndisplacement combined with polymerization/nicking modules. While the authors\npropose an ambitious framework for achieving arbitrarily programmable multistability,\nseveral critical aspects require clarification to fully support the claims. Notably, the\nextent and scalability of this \"arbitrarily programmable\" capability—specifically, how\nmany distinct steady states can be robustly engineered and maintained—remain\nunclear. Regrettably, the current evidence does not sufficiently validate the authors'\ncentral thesis.”\nResponse: We sincerely thank the reviewer for raising a key issue concerning the extent\nand scalability of the system’s “arbitrarily programmable” nature. In response, we have\nelucidated this issue from theoretical foundations and experimental limitations using\nAMN-1 as a representative case study:\n(i) Theoretical foundations\nAs a proof-of-concept, our work demonstrates that arbitrary multistability can be\nrealized through the reversible displacement reaction-mediated competition between\npositively autoregulatory DNA polymerization/nicking modules. The relative fraction\nbetween the two modules can, in principle, be continuously tuned to support a broad\nand densely populated range of steady states.\nIn our framework, the term “arbitrarily programmable” encompasses two key\ndimensions: Positional arbitrariness: The system can stabilize at any fractional\ncomposition between the two competing modules, thereby forming a continuous and\ntunable multistable manifold. Numerical arbitrariness: In theory, each unique\ncompositional fraction corresponds to a distinct steady state, enabling the potential for\nan arbitrarily large number of distinguishable states.\nTo validate these claims, we performed a Linear Stability Analysis (LSA) based on\nour kinetic model. This analysis revealed the existence of thousands of steady states,\neach characterized by net fluxes below 1×10−10 nM·h−1 for all molecular species—\nqualifying them as distinct attractors in the system’s compositional space. From this set,\nwe uniformly sampled 216 representative steady states and plotted them in the new Fig.\n3u (see Comment 6), which shows a dense distribution along a linear manifold. This\nprovides strong theoretical support for a continuum of programmable steady states in\nAMN-1.\n(ii) Experimental limitations\nWhile the theoretical capacity for multistability is extensive, we acknowledge that\nexperimental realizability is bounded by factors such as resolution limits, biochemical\nconstraints, and the robustness of individual steady states. In particular, the ability to\nexperimentally distinguish a greater number of steady states depends critically on\nminimizing steady-state drift and enhancing system stability over time.\nTo this end, we conducted a series of experimental optimizations focusing on key\nparameters such as enzyme ratios and nucleotide availability: (i) Polymerase-to-\nnickase ratio perturbations: As presented in Supplementary Figs. 18–23, AMN-1\nmaintains its steady-state behavior under variations in polymerase-to-nickase (pol/nick)\nratios, indicating good steady-state retention despite enzymatic fluctuations. (ii) dNTP\nconcentration effects: In the revised manuscript, we present new data evaluating the\nimpact of dNTP concentrations on steady-state stability (new Supplementary Figs. 13–\n15, see Comment 6). These results indicate that dNTP depletion is a major factor\nlimiting long-term multistability. By increasing the dNTP concentration to 0.2 mM,\nAMN-1 is capable of maintaining a steady-state composition with minimal deviation\n(~0.07%·h−1 by 48 h, new Supplementary Fig. 13).\nEmpirically, we observed that AMN-1 accesses a compositional range of ~5% to 92%\nfor the BA'+AB' module (Fig. 3o and 3s). Assuming a conservative distinguishability\nthreshold based on drift rate of 0.07 %·h−1 and a multistability duration of 48 h,\nwe\nestimate the number of resolvable steady states as follows: (92%－5%)/0.07 %·h−1 ×\n48 h ≈ 26. This estimate reflects a lower bound under current experimental conditions,\nand we anticipate that the number of experimentally distinguishable steady states could\nbe significantly increased through further biochemical stabilization and improved\ndetection resolution.\nWe have revised the manuscript to provide a more rigorous and precise definition of\n“arbitrarily programmable” in the context of our system (page 11). Also, a detailed\ndiscussion of both the theoretical framework and experimental considerations\nunderpinning the scalability of programmable multistability is now included in the main\ntext (pages 11, 12 and 21).\nComment 16: “The current schematic is challenging to interpret. To improve\naccessibility, the authors should explicitly label the 3′ and 5′ ends of DNA strands to\nindicate the directions of polymerase extension. Additional molecular details in the\ndiagram would help readers visualize the design principles.”\nResponse: We thank the reviewer for this valuable suggestion regarding the clarity of\nour schematics. In response, we have revised the network diagrams in the updated Figs.\n2 and 5a (see Comment 3) to explicitly label the 3' and 5' ends of all DNA strands,\nthereby clarifying the directionality of polymerase-mediated strand extension.\nTo further enhance the accessibility of our design, we have included new\nSupplementary Figs. 1 and 36, which present detailed sequence-level diagrams of the\nsystem’s core components. These additions provide a more comprehensive view of the\nmolecular architecture and better illustrate the design principles underlying the\nnetwork’s function.\nFig. 2｜Schematic design of a basic two-module network, AMN-1. AMN-1 was constructed\nbased on reversible displacement reaction-mediated competition between two autoregulatory\nmodules, AA'+BB' and BA'+AB'. Panel I—Schematic showing dNTP-fueled DNA\npolymerization/nicking reactions. Panel II—Schematic representation of AMN-1.\nSupplementary Fig. 1. Detailed sequences of all DNA strands used in AMN-1. Each strand\nis annotated with color-coded and numbered domains representing distinct functional elements,\nincluding toeholds, branch migration domains, and other structural elements.\nSupplementary Fig. 36. Detailed sequences of all DNA strands used in AMN-2. Each strand\nis annotated with color-coded and numbered domains representing distinct functional elements,\nincluding toeholds, branch migration domains, and other structural elements.\nComment 17: “The domain-level labels on DNA complexes are insufficiently clear.\nRevising the figures to highlight domain names with distinct colors or fonts would aid\nreader comprehension.”\nResponse: We thank the reviewer for this helpful and constructive suggestion. In\nresponse, we have revised Figs. 2 and 5a (see Comments 16 and 3) to improve the\nclarity of domain-level annotations in the DNA complexes. Specifically, distinct\ndomains are now labeled using a consistent numbering system. To further enhance\nvisual comprehension, we have applied corresponding color schemes and numerical\nidentifiers in the associated sequence-level representations. These enhanced\nannotations are presented in the updated Supplementary Figs. 1 and 36 (see Comment\n16). These revisions substantially improve the clarity and readability of the figures.\nComment 18: “The description leaves uncertainty about whether the four duplexes\n(AA', BB', BA', AB') are pre-annealed or formed spontaneously from single strands.\nClarifying this hierarchy in the text is essential for reconstructing the system’s\noperational logic.”\nResponse: We appreciate the reviewer's valuable suggestion. Indeed, this clarification\nis essential for accurately conveying the assembly hierarchy and operational logic of\nthe network. In response, we have revised the relevant description in the revised\nmanuscript (page 6) to clarify that the four duplexes (AA', BB', BA', and AB') are pre-\nannealed prior to their incorporation into the system.\nComment 19: “The distinction between AA'-P1/AB'-P2 and AA'/AB' remains unclear.\nIf these complexes are functionally equivalent, consistent naming would reduce\nconfusion. If not, a side-by-side comparison of their structural differences should be\nprovided.”\nResponse: We appreciate the reviewer's insightful comment. In the context of the AMN\ndesign, the complexes called AA'-P /AB'-P and AA'/AB' are functionally equivalent.\n1 2\nWe recognize that inconsistent naming might cause unnecessary confusion. In response,\nwe have revised the manuscript to adopt consistent terminology throughout, replacing\n“AA'-P ” and “AB'-P ” with “AA'” and “AB'”, respectively. This change has been\n1 2\nimplemented uniformly across all figures, main text, and supplementary materials to\nensure clarity and reduce potential ambiguity.\nComment 20: “The manuscript does not sufficiently explain the driving forces behind\nthe reversible transitions (AA'+BB' -- BA'+AB'). Specifically: What thermodynamic or\nkinetic factors enable bidirectional conversion? How does AA'P1 formation selectively\ninhibit the forward reaction? A mechanistic diagram or kinetic model would help justify\nthis design.”\nResponse: We sincerely thank the reviewer for this insightful and essential comment\nregarding the mechanistic underpinnings of the reversible transitions between modules\nAA'+BB' and AB'+BA', as well as the role of AA'P formation in modulating the\n1\ndirectionality of this reaction. We fully agree that a more comprehensive explanation\nof the relevant thermodynamic and kinetic factors, accompanied by mechanistic\ndiagrams, would enhance the clarity and rigor of our design. We have elucidated this\nissue from the following aspects:\n(i) Thermodynamic basis for the reversible displacement reaction\nThe interconversion between AA'+BB' and AB'+BA' (AA'+BB' ⇌ AB'+BA')\nproceeds via unbiased branch migration rather than conventional toehold-mediated\nstrand displacement. This mechanism is inherently reversible and governed by\nthermodynamic equilibrium, with the rate and extent of interconversion determined by\nthe relative duplex stabilities and reaction conditions (e.g., temperature, ionic strength).\nTo ensure the thermodynamic reversibility of the system, we employed NUPACK\nsimulations to guide the design of the four duplex modules. As shown in the new\nSupplementary Fig. 2, the calculated Gibbs free energies (ΔG) of the four engineered\nduplexes are approximately −16.33 kcal·mol−1, indicating near isoenergetic states. Thus,\nΔG of the reversible displacement reaction was calculated to be approximately 0,\nsuggesting that the reaction is a bidirectional conversion. This is supported by the\nexperimental results, where the system with different initial compositions converges\ninto a single equilibrated state under non-enzymatic conditions (new Supplementary\nFig. 5, see Comment 2).\n(ii) AA'P -induced kinetic modulation\npathway can be achieved through kinetic modulation. Specifically, we leveraged the\npolymerization of AA' into a more stable complex, AA'P , to introduce such kinetic\n1\nasymmetry. As demonstrated by NUPACK simulation (new Supplementary Fig. 2), the\nresulting AA'P exhibits markedly increased thermodynamic stability (ΔG = −51.70\n1\nkcal·mol−1) and does not participate in further reversible displacement reaction with\nBB'. The formation of AA'P effectively sequesters AA' from the reaction pool, thereby\n1\nsuppressing the forward pathway (AA'+BB' → AB'+BA') and favoring the reverse\nprocess. This process drives the conversion of AB'+BA' into AA'P +BB'. Over time,\n1\nthe reaction network is predominantly populated by module AA'P +BB', as further\n1\nsupported by NUPACK analysis in the new Supplementary Fig. 3b (see Comment 5).\nBy selectively depleting AA' through its irreversible incorporation into a highly\nstable complex, we establish a kinetically asymmetric landscape that allows controlled\nmodulation of network states in response to external perturbations.\nWe have included these discussions in the new Supplementary Section 2 and\nmentioned them in the revised main text (page 5).\nSupplementary Fig. 2. Free energy landscape of AMN-1 network reconfiguration before\nand after DNA polymerization, based on NUPACK simulations. Before DNA\npolymerization, the ΔG values of modules BA'+AB' and AA'+BB' are equal, suggesting the\nreversible module interconversion. However, after DNA polymerization, AA' was extended to\ngenerate elongated product AA'P with more negative free energy, thus suppressing the forward\n1\nreaction (AA'+BB' → BA'+AB').\nComment 21: “The assertion that dNTP hydrolysis suppresses spontaneous\ninterconversion lacks experimental or theoretical support in the \"Design\" section. The\nauthors should explicitly describe how dNTP consumption is coupled to module\nstabilization.”\nResponse: We thank the reviewer for this insightful comment, which highlights the\nneed to clarify the mechanistic coupling between dNTP hydrolysis and the suppression\nof spontaneous interconversion among network modules. In response, we have revised\nthe manuscript to address this point from both a theoretical and experimental\nperspective, as outlined below:\n(i) Theoretical framework\ndNTP hydrolysis provides the driving force for template-directed DNA\npolymerization, catalyzed by DNA polymerase according to the reaction: (DNA) +\nn\ndNTP → (DNA) + PP. In our system, duplexes AA' and AB' serve as templates for\nn+1 i\npolymerase-mediated extension, yielding the polymerized products AA'P and AB'P ,\n1 2\nrespectively. These extended duplexes exhibit enhanced thermodynamic stability and\nare rendered inert to the reversible displacement reaction. As a result, AA' and AB' are\nselectively and irreversibly sequestered from the dynamic pool, thereby suppressing the\nforward (AA'+BB' → AB'+BA') or reverse (AB'+BA' → AA'+BB') interconversion\npathway. This mechanism enables the network to establish a dNTP-driven, non-\nequilibrium steady state, thereby preventing relaxation toward thermodynamic\nequilibrium.\n(ii) Experimental evidence\nTo experimentally validate the role of dNTP hydrolysis in modulating\ninterconversion dynamics, we examined the behavior of AMN-1 initialized in the\nAA'+BB' state under varying concentrations of dNTPs. In the absence of dNTPs\n(Supplementary Fig. 15a), the system evolves toward equilibrium, consistent with\nspontaneous interconversion between modules. However, as dNTP concentrations are\nprogressively increased, the rate of interconversion diminishes, and the system\nstabilizes near its initial composition. At sufficiently high dNTP levels, interconversion\nis effectively suppressed (new Supplementary Figs. 13−15, see Comment 6), indicating\nthat dNTP-driven polymerization sequesters the reactive duplexes required for\ndisplacement-based interconversion, thereby stabilizing the network composition.\nTogether, these theoretical and experimental findings provide strong support for the\nconclusion that dNTP hydrolysis, through its role in driving DNA polymerization,\nactively suppresses module interconversion by forming stable, polymerized species\n(AA'P and AB'P ) and depleting the dynamic intermediates required for the reversible\n1 2\ndisplacement reaction. These clarifications have now been incorporated into the revised\nmain text (pages 5, 6 and 11) to enhance mechanistic transparency and reinforce the\ncoupling between dNTP consumption and network stabilization.\nComment 22: “The roles of strands A (BHQ2-labeled), A' (FAM), and B' (ROX) are\nambiguously described. A dedicated subsection or table mapping strands to their\nmodifications would prevent misinterpretation of the monitoring strategy.”\nResponse: We sincerely thank the reviewer for pointing out the ambiguity in our\noriginal description of the strands A (BHQ2-labeled), A' (FAM-labeled), and B' (ROX-\nlabeled), particularly about their roles within the fluorescence monitoring strategy. To\naddress this concern and improve clarity, we have introduced the new Supplementary\nTables 1 and 3, which provide a clear mapping of each strand to its corresponding\nfluorophore or quencher modification, along with its functional role in the network. In\naddition, we have revised the paragraph titled \"Quantitative assessment of the\ncompositions in networks AMN-1 and AMN-2 through fluorescence measurement\" in\nthe Supplementary Information (pages 9 and 10). These updates have been mentioned\nin the main text (pages 6 and 16).\nSupplementary Table 1. Sequence identities, chemical modifications, duplex\nconfigurations, and corresponding fluorescence outputs of DNA duplexes used in AMN-1\nStrand A Strand B\n(Internal BHQ2 label) (No modification)\nStrand A' Duplex AA' Duplex BA'\n(Internal FAM label) (FAM-BHQ2, quenched) (FAM)\nStrand B' Duplex AB' Duplex BB'\n(Internal ROX label) (ROX-BHQ2, quenched) (ROX)\nSupplementary Table 3. Sequence identities, chemical modifications, duplex\nconfigurations, and corresponding fluorescence outputs of DNA duplexes used in AMN-2\nStrand C Strand D\n(Internal BHQ2 label) (No modification)\nStrand C' Duplex CC' Duplex DC'\n(Internal FAM label) (FAM-BHQ2, quenched) (FAM)\nStrand D' Duplex CD' Duplex DD'\n(Internal ROX label) (ROX-BHQ2, quenched) (ROX)\nStrand E' Duplex CE' Duplex DE'\n(Internal Cy5 label) (Cy5-BHQ2, quenched) (Cy5)\nComment 23: “The significant difference in equilibration kinetics between the two\ninitial conditions (100% AA'+BB' vs. BA'+AB') appears inconsistent with the proposed\nthermodynamic framework. The authors should address whether this asymmetry arises\nfrom kinetic traps or unintended design biases. (Fig. 3a)”\nResponse: We thank the reviewer for this perceptive observation. The difference in\nequilibration kinetics between the two initial conditions (100% AA'+BB' vs 100%\nBA'+AB') was indeed unintentional and appears to stem from subtle design biases\nrather than kinetic traps:\nOur design strategy employed NUPACK software to design sequences with nearly\nperfect thermodynamic symmetry, followed by experimental validation and iterative\nrefinement. According to the NUPACK’s thermodynamic predictions (new\nSupplementary Fig. 3a, see Comment 5), the system should yield an equilibrium state,\nwith equivalent concentrations of modules AA'+BB' and BA'+AB' (theoretical\nequilibrium constant K ≈ 1). This implies that the forward and reverse reaction rate\nconstants are expected to be approximately equal.\nHowever, as shown in Fig. 3b and 3c (original Fig. 3a), the experimentally measured\nequilibrium deviates from this prediction: the system equilibrates to a state in which the\nconcentration of BA'+AB' exceeds that of AA'+BB', corresponding to a measured\nequilibrium constant of K ≈ 6. This indicates a significant asymmetry in the effective\nforward and reverse reaction rates, inconsistent with the thermodynamic symmetry\npredicted by the design.\nWhile all four species share the same core binding domain (1/1'), they differ in their\nauxiliary sequences. We hypothesize that this kinetic asymmetry arises from sequence-\ndependent features not fully captured by the current thermodynamic models used in\nNUPACK. Potential contributing factors include local variations in GC content,\nsecondary structure formation in single-stranded regions, or positional effects at domain\njunctions—each of which may differentially influence strand displacement kinetics and\nhybridization efficiency.\nThese findings and hypotheses have been included in the new Supplementary Section\n2 and mentioned in the revised main text (page 7).\nComment 24: “The claim that rTAA' hybridization upregulates AA'+BB' requires\nmechanistic elaboration. How does AA'-rTAA' formation preferentially stabilize\nspecific modules? A diagram illustrating the energy landscape or binding competition\nwould strengthen this argument.”\nResponse: We thank the reviewer for this insightful and constructive comment. In\nresponse, we have substantially revised the Supplementary Information to include a\nmore detailed mechanistic explanation (new Supplementary Section 2), along with a\nschematic illustration of the associated energy landscape (new Supplementary Fig. 25,\nsee Comment 5).\nSpecifically, we show that the addition of rT leads to the formation of a\nAA'\nthermodynamically stable duplex (rT AA') with AA' (ΔG ≈ −40.88 kcal·mol−1,\nAA'\ncalculated using NUPACK). This hybridization effectively sequesters AA' from the\ndynamic pool, thereby reducing its participation in the reversible strand displacement\nreaction (AA'+BB' ⇌ AB'+BA'). As a result, the equilibrium shifts toward the reactant\nmodule (AA'+BB'), leading to its preferential stabilization.\nThese clarifications have been mentioned in the main text (page 13).\nComment 25: “The manuscript would benefit from exploring scenarios involving\nsimultaneous or sequential RNA inputs: How would competing RNA triggers at varying\nconcentrations influence state selection? Is the system fully reversible if a second RNA\nis introduced after rTAA'-induced stabilization?”\nResponse: We thank the reviewer for this insightful and important suggestion regarding\nthe impact of simultaneous or sequential RNA triggers on the state-transition dynamics\nof our system. In response, we conducted a series of additional experiments to\nsystematically investigate how the introduction of multiple RNA inputs—either\nconcurrently or in sequence—affects state selection and reversibility.\nUnder simultaneous input conditions, we investigated the behavior of AMN-1in\nresponse to varying molar ratios of two competing RNA triggers, rT and rT (new\nAA' AB'\nSupplementary Fig. 33). These experiments revealed that the steady state of the system\nis highly sensitive to the relative concentrations of the competing RNA species,\nindicating a concentration-dependent selection mechanism. These findings highlight\nthe tunability and programmability of the network, wherein input stoichiometry can be\nused to bias the system toward specific steady states.\nUnder sequential input conditions, we first exposed the system to rT , which\nAA'\nrobustly induced the corresponding steady state. Subsequent introduction of RNase H\neffectively removed rT . Notably, the system retained the rT -induced state despite\nAA' AA'\nthe degradation of the initial trigger, indicating a form of memory. However, upon the\naddition of rT , the system transitioned back to the original steady state (new\nAB'\nSupplementary Fig. 34). These observations confirm that the system exhibits both\nmemory and reversibility, depending on the sequence of competing RNA triggers.\nWe believe these additional experiments significantly enhance the manuscript by\ndemonstrating the system’s capacity for dynamic adaptation, input prioritization, and\nprogrammable multistability. The new data and associated discussions have been\nincluded in the revised manuscript (page 14).\nSupplementary Fig. 33. Dynamic behavior of AMN-1 in response to competing RNA\ntriggers at varying concentration ratios. a,b, Time-dependent fluorescence changes of\nAMN-1 in S after treatment with RNase H (3.00 U·mL−1) and a mixture of RNA triggers\n1-1\nrT and rT at different concentrations: 0.9 μM rT and 0.1 μM rT (a), and 0.7 μM rT\nAB' AA' AB' AA' AB'\nand 0.3 μM rT (b). c,d, Time-dependent concentration changes of AA', AB', BA', and BB',\nAA'\nderived from fluorescence data in panels a and b using respective calibration curves\n(Supplementary Fig. 4).\nSupplementary Fig. 34. Reversible state transitions of AMN-1 in steady state S . a, Time-\n1-3\ndependent fluorescence changes associated with AMN-1 in steady state S upon sequential\n1-3\naddition of RNA trigger rT (0.2 μM), enzyme RNase H (3.00 U·mL−1), and the second RNA\nAA'\ntrigger, rT (1.7 μM). b, Time-dependent concentration changes of AA', AB', BA', and BB'\nAB'\nderived from fluorescence data presented in a using the corresponding calibration curves\n(Supplementary Fig. 4).\nComment 26: “The observed convergence near initial concentrations raises questions\nabout true multistability. To substantiate claims, the authors should: Define\nquantitative criteria for distinguishing steady states. Compare experimental results\nwith theoretical predictions of bistable/multistable systems”\nResponse: We sincerely thank the reviewer for this thoughtful and constructive\ncomment. We fully agree that a rigorous definition of steady states, alongside a\nsystematic comparison between experimental observations and theoretical models, is\nessential to validate the presence of true multistability. In response, we have addressed\nthese points based on the following clarifications and supplied new results:\n(i) Clarification of convergence behavior and justification of true multistability\nTo address this concern, we have carried out a detailed analysis of the compositional\ndrift. Across all initial conditions, AMN-1 exhibits a two-phase dynamic behavior. The\ninitial phase (~0.5 h) is characterized by a short relaxation, during which\npolymerization/nicking reactions exert an increasing, suppressive influence on the\nreversible displacement reaction. This is followed by a second prolonged phase,\nwhich corresponds to the operational steady state. In this regime,\npolymerization/nicking reactions have fully established kinetic dominance, effectively\nhalting further spontaneous interconversion between network modules. However, this\nregime is not entirely static—compositional changes occur at a markedly slower rate.\nTo quantify this drift within the steady-state regime, we calculated the drift rate as\nthe change in the concentration of each module per hour, normalized by the total\nconcentration of all modules. This calculation provides a practical metric to evaluate\nthe long-term stability of the system.\nUnder all initial conditions, AMN-1 reliably settles into composition-specific steady\nstates (S to S ) with slight compositional drift (revised Fig. 3t, see Comment 6). The\n1-1 1-5\nmost considerable drift is observed in S , with a maximum rate of ~0.5%·h−1 over a\n1-1\n12-hour period. Comparable drift phenomena have been reported in conventional\nmultistable biochemical systems (ref. 58—Nat. Commun. 10, 4636 (2019), ref. 61—\nNat. Chem. 11, 829–838 (2019), ref. 62—Nat. Chem. 14, 1224–1232 (2022), ref. 63—\nProc. Natl. Acad. Sci. USA 109, E3212–E3220 (2012) and ref. 64—Nat. Commun. 7,\n13474 (2016)).\nTo investigate the source of this drift and assess long-term stability, we conducted\nextended experiments under S while varying dNTP concentrations (new\n1-1\nSupplementary Fig. 13, see Comment 6). At the original dNTP concentration of 0.1\nmM, the drift becomes pronounced after approximately 36 h, reaching a drift rate of\n~0.7%·h−1 during the final 12-hour window. In contrast, increasing dNTP concentration\nto 0.2 mM significantly mitigates the drift (~0.07%·h−1 by 48 h), indicating that fuel\ndepletion is the primary cause of this drift. Under conditions of adequate dNTP supply,\nthe AMN-1 network is capable of maintaining a stable, steady-state composition with\nminimal drift, implying the robustness of AMN-1’s multistability.\nTo further assess the energy dependence of steady-state stability, we tested a broader\ndNTP range (0–0.02 mM, new Supplementary Figs. 14 and 15, see Comment 6). At\nsuboptimal dNTP levels, multistability is compromised, and compositional drift\nincreases markedly. Conversely, higher dNTP concentrations enhance both the\npersistence and fidelity of the steady states.\nThis analysis and results have been included in the revised manuscript (page 11) and\ndetailed in the new Supplementary Section 5.\n(ii) Quantitative criteria for steady-state identification\nTo objectively identify and distinguish steady states, we have introduced a set of\nquantitative criteria based on compositional stability. A system state is classified as a\nsteady state if it satisfies the following conditions:\n It undergoes an initial transient phase marked by significant compositional\nfluctuations;\n Subsequently, it enters a second, long-lasting phase during which the compositional\nprofile remains consistent within a tolerance of 0.5%·h−1.\nThese criteria enable reproducible classification of steady states across all\nexperimental conditions and help clarify the distinctions between true multistability and\ntransient convergence.\nA detailed description of these criteria has been included in the Supplementary\nInformation (page 10).\n(iii) Comparison of experimental results with theoretical predictions\nEmpirically, we observed that AMN-1 accesses a compositional range of ~5% to 92%\nfor the BA'+AB' module (Fig. 3o and 3s). Assuming a conservative distinguishability\nthreshold based on drift rate of 0.07%·h−1 and a multistability duration of 48 h, we\nestimate the number of resolvable steady states as follows: (92% − 5%)/0.07 %·h−1 ×\n48 h ≈ 26. This estimate provides a lower bound on the number of experimentally\ndistinguishable steady states under current operating conditions.\nTo compare with theoretical predictions, we performed a Linear Stability Analysis\n(LSA) based on our kinetic model, which revealed a large number of stable steady states.\nApplying a conservative threshold for compositional flux (<1×10−10 nM·h−1 per\nspecies), we identified thousands of theoretically stable states, each representing a\ndistinct attractor in the system’s compositional space. From these, we uniformly\nsampled 216 representative steady states. We plotted their distributions in the new Fig.\n3u (see Comment 6). These states densely populate a linear manifold within the\nexperimentally accessible range, providing strong theoretical support for the emergence\nof approximate continuous multistability in AMN-1.\nCollectively, these experimental and theoretical results consistently point to the\npresence of an approximate continuum of compositional attractors, reinforcing our\ninterpretation of AMN-1 as a system capable of biologically relevant arbitrary\nmultistability. While experimental resolution currently limits the full realization of this\ntheoretical landscape, we anticipate that the number of resolvable steady states could\nbe substantially increased through further optimization.\nThe experimental estimation and theoretical simulation of the system’s possible\nsteady states have been discussed in the main text (page 12).\nComment 27: “The gel lanes for DNAzyme 2 products lack sufficient contrast (Fig.\n6b). Repeating the experiment with adjusted exposure times or alternative visualization\nmethods (e.g., SYBR Gold staining) would improve credibility.”\nResponse: We appreciate the reviewer’s comment regarding the low contrast in the gel\nlanes for the products presented in Fig. 6b. In response, we have repeated the gel\nelectrophoresis using optimized imaging conditions. Specifically, we adjusted the\nexposure time and the concentration of the staining reagent to enhance the sensitivity\nand contrast of the DNA bands. These changes have greatly enhanced the clarity and\nresolution of the gel image. The updated gel image has been included in the revised\nmanuscript as new Fig. 6b.\nRe: Manuscript ID: NCOMMS-25-07099-A\nOriginal title: \"Arbitrary Multistability in DNA Replication Networks\"\nRevised title: \"Continuously Tunable Multistability in DNA Replication Networks\"\nThe following point-by-point listed corrections/explanations were introduced into\nthe paper:\nReviewer #1:\nComment 1: “The authors have done convincing and in-depth revisions. I also\nappreciate clarifications on conceptual novelty! For me, the paper is ready for\npublication.”\nResponse: We sincerely thank the reviewer for the positive evaluation and for\nrecognizing the conceptual clarifications and revised analyses. We greatly appreciate\nthe recommendation for publication.\nReviewer #3:\nComment 2: “The authors have addressed most of my concerns, and thanks for the\ntheoretical clarifications regarding scalability (Linear Stability Analysis) and\nexperimental optimizations (dNTP titration, enzyme ratios). However, the central claim\nof “arbitrarily programmable” multistability lacks experimental substantiation, as\ndetailed below:\n1. Theoretical claims ≠ experimental innovation in DNA programming: Breakthroughs\nin DNA-programmed systems must prioritize empirical validation. While LSA suggests\n\"thousands of steady states\" (Response §i), only 26 resolvable states (Response §ii)\nwere experimentally achieved—a marginal figure that contradicts \"arbitrary\"\nprogrammability. This gap between theory and practice undermines the innovation’s\ncore promise.\n2. \"Arbitrary\" requires experimental scalability: Programmable biological systems are\njudged by engineered outcomes, not theoretical potential. The reported 5–92%\ncompositional range (Response §ii) with high drift (0.07%/h) confirms severe\nbiochemical constraints, not arbitrary control. DNA computing should transcend\nlimitations (e.g., via error correction), not rationalize them.\n3. Terminological rigor concerns: Redefining \"arbitrary\" as \"context-dependent\"\n(Revised Manuscript p. 11) dilutes scientific precision. Equating \"arbitrary\nprogrammability\" with \"26 states under narrow conditions\" (Response §ii)\nmisrepresents scalability and sets a misleading precedent for the field.”\nResponse: We sincerely thank the reviewer for the thorough and constructive\nevaluation. The comments have prompted us to refine our conceptual framing, enhance\nterminological rigor, and conduct additional experiments to more precisely define the\noperational limits of the system. A summary of these revisions and our detailed\nresponses is provided below.\n1. Terminological precision and conceptual alignment\nWe fully agree with the reviewer that the previous use of the term “arbitrary\nmultistability” lacked the required rigor and risked implying an experimental capability\nbeyond what our current data supports. In response to this concern, we have replaced\n“arbitrary multistability” with the more precise term “continuously tunable\nmultistability” throughout the revised manuscript and Supplementary Information.\nRevising this terminology serves two purposes:\nFirst, the revised terminology more accurately reflects our theoretical and\nexperimental results. LSA results indicate that our network architecture can, in principle,\ngenerate a continuum of multistable steady states with continuously tunable\ncompositions. Consistent with this prediction, our experiments realize a continuously\ntunable range of steady-state compositions. Together, theoretical simulations and\nexperimental observations converge on the presence of continuously tunable\nmultistability—although the number of experimentally resolvable steady states is\nsmaller than theoretical expectation. This discrepancy arises from practical constraints\ninherent to biochemical reaction systems—including dNTP depletion, finite enzyme\nactivities, and limited measurement resolution—which collectively limit the\nexperimental accessibility of the whole theoretical state space.\nSecond, the revised terminology preserves the mechanistic novelty and conceptual\nsignificance of our network architecture without overstating its current experimental\nrealizability. The core mechanistic advance—a far-from-equilibrium, competition-driven\ndesign that enables continuous tuning of steady-state compositions—remains\nfully intact. This design principle fundamentally distinguishes our architecture from\nconventional multistable networks, which typically support only a small number of\ndiscrete fixed points—by contrast, our framework enables a continuous range of\ntunable steady-state compositions, representing a distinct and previously unexploited\nmode of multistability.\nThus, this terminological refinement serves not only to preserve the conceptual and\nmechanistic significance of the design but also to avoid overstating the experimental\naccessibility steady states under current laboratory conditions. All relevant sections of\nthe manuscript have been updated to ensure terminological consistency, conceptual\ncoherence, and alignment between theoretical claims and experimental results.\n2. Experimental scalability\nTo address the reviewer’s concerns regarding experimental scalability, we conducted\na series of additional optimization experiments. The results clarify the biochemical\nconstraints that define the system’s operational boundaries and provide a more\nquantitative, transparent assessment of the number of experimentally resolvable steady\nstates.\ni) dNTP concentration optimization\nIn the original manuscript, increasing the dNTP concentration from 0.10 mM to 0.20\nmM significantly reduced the drift in steady state S . To further assess whether\n1-1\nadditional increases in substrate availability could suppress this drift, we performed an\nadditional experiment on steady state S (100% AA' + BB') at a higher dNTP\n1-1\nconcentration of 0.30 mM (new Supplementary Fig. 14c and 14f). The drift observed\nat 0.30 mM dNTPs showed no appreciable improvement relative to that at 0.20 mM.\nThis indicates that the effect of dNTP supplementation saturates near 0.20 mM under\nour experimental conditions. Thus, 0.20 mM was chosen for the optimal dNTP\nconcentration for subsequent experiments. The corresponding discussion has been\nincorporated into the revised manuscript (p. 11).\nSupplementary Fig. 14. Dynamic behaviors of network CTMN-1 in steady state S at the\n1-1\ndNTP concentrations of 0.1, 0.2, and 0.3 mM. a−c, Time-dependent fluorescence changes of\nCTMN-1 in S at different dNTP concentrations: 0.1 mM (a), 0.2 mM (b), and 0.3 mM (c). d−f,\n1-1\nTime-dependent concentration changes of AA', AB', BA', and BB' derived from the fluorescence\ndata in a−c using corresponding calibration curves (Supplementary Fig. 4).\nii) Enzyme-ratio optimization\nUsing the optimal dNTP concentration (0.2 mM), we next investigated the influence\nof varying enzyme ratios on the drift observed at steady state S (100% AA'+BB'). To\n1-1\nthis end, we increased the DNA polymerase concentration from 30.0 U·mL−1 to 35.0\nU·mL−1 while decreasing the nicking enzyme concentration from 20.0 U·mL−1 to 15.0\nU·mL−1 (new Supplementary Figs. 24 and 25), thereby shifting the system toward a\npolymerase-rich, nicking-limited regime. Under these conditions, no meaningful\nreduction in drift was observed relative to the original enzyme ratio. These results\nsuggest that simple stoichiometric adjustments of enzymes do not overcome the drift\nconstraints of the current architecture. The corresponding discussion has been\nincorporated into the revised manuscript (p. 12).\nSupplementary Fig. 25. Time-dependent compositional dynamics of CTMN-1 in steady state\nS under fixed dNTP concentration at 0.2 mM while varying DNA polymerase and nicking\n1-1\nendonuclease concentrations. a−c, Experimentally determined (solid curves) time-dependent\nconcentration profiles of AA', AB', BA', and BB' derived from fluorescence data shown in\nSupplementary Fig. 24 and converted using the corresponding calibration curves (Supplementary\nFig. 4). Experimental conditions: CTMN-1 in steady state S with 0.2 mM dNTP and either 30\n1-1\nU·mL−1 DNA polymerase and 15 U·mL−1 nicking endonuclease (a); 35 U·mL−1 DNA polymerase\nand 15 U·mL−1 nicking endonuclease (b); or 35 U·mL−1 DNA polymerase and 20 U·mL−1 nicking\nendonuclease (c). Dashed curves represent the corresponding computationally simulated results. (d)\nDrift rates of the CTMN-1 network as a function of polymerase-to-nickase ratio, generated by fixing\ndNTP concentration at 0.2 mM while varying DNA polymerase and nicking endonuclease levels.\niii) System-wide drift quantification\nIn the original manuscript, the experimental scalability of steady states was estimated\nsolely from the network's drift rate at steady state S\n1-1\n. This single-state metric does not\nadequately capture stability variations across the entire multistable landscape. To\nestablish a global stability metric spanning the entire multistability range, we instead\nquantified the average drift rate across five representative steady states covering the\nexperimentally accessible multistable interval. This system-wide metric offers a more\nappropriate proxy for estimating the number of experimentally resolvable steady states\nthan the previous single-state measurement.\nAccordingly, under optimized biochemical conditions (0.20 mM dNTP, 30.0 U·mL−1\nDNA polymerase, 20.0 U·mL−1 nicking enzyme), we measured the drift rates at five\nrepresentative steady states (S\n1-1\nto S\n1-5\n). The resulting drift rates—0.16 %·h−1,\n0.13 %·h−1, 0.04 %·h−1, −0.06 %·h−1, and −0.05 %·h−1 for S\n1-1\nto S\n1-5\n, respectively—\nare all substantially lower than those observed at 0.10 mM dNTP. The mean drift rate\nacross these five states (0.04 %·h−1) was used as the global stability parameter for\nsubsequent scalability estimation. These results are shown in new Supplementary Fig.\n17 and new Fig. 3o–3u. The original data have been moved to Supplementary Fig. 12,\nand the implications are discussed in the revised manuscript (p. 11).\nSupplementary Fig. 17. Fluorescence kinetics of CTMN-1 in steady states S − S at the\n1-1 1-5\ndNTP concentration of 0.2 mM. a−e, Time-dependent fluorescence changes of CTMN-1 in steady\nstates S − S .\n1-1 1-5\nFig. 3｜Continuously tunable multistability of CTMN-1. a, Schematic representation of CTMN-\n1 without autoregulatory loop, in the absence of DNA polymerase, nicking endonuclease, and\ndNTPs. b,c, Time-dependent concentration changes of AA', BB', BA', and AB' in the network\ninitiated by two distinct compositions: 100% AA'+BB' (b) and 100% BA'+AB' (c). d,\nComputationally derived reversible displacement rates. e, Schematic representation of a modified\nCTMN-1, where A' was replaced by A', a variant without template domain P ', leading to a single\n1\nautoregulatory loop of BA'+AB'. f,g, Time-dependent concentration changes of AA', BB', BA', and\nAB' in the network with two initial compositions: 100% AA'+BB' (f) and 100% BA'+AB' (g). h,i,\nCorresponding\ncomputationally derived reversible displacement rates (h) and\npolymerization/nicking rates associated with AB' (i). j, Schematic representation of another\nmodified CTMN-1 with a single autoregulatory loop of AA'+BB'. Here, B' was replaced by B', a\nvariant without template domain P '; k,l, Time-dependent concentration profiles of AA', BB', BA',\n2\nand AB' in the network initiated from two initial compositions: 100% AA'+BB' (k) and 100%\nBA'+AB' (l). m,n, Computationally derived reversible displacement rates (m) and\npolymerization/nicking rates associated with AA' (n). o−s, Time-dependent concentration changes\nof AA', BB', BA', and AB' in the CTMN-1 with two competing autoregulatory loops under five\ninitial compositions, where the proportion of AA'+BB' varies from 100% (o), 75% (p), 50% (q),\nand 25% (r) to 0% (s), with corresponding variations in BA'+AB' from 0%, 25%, 50%, and 75% to\n100%, under 0.2 mM dNTPs; The network converges to five different steady states (S to S ). t,\n1-1 1-5\nDrift rates of BA'+AB' by 24 h under these steady states, where the drift rate is defined as the change\nin the concentration of each module per hour, normalized by the total concentration of all modules.\nu, Linearly distributed steady states. Uniformly sampled theoretically predicted steady states of\nCTMN-1 (216 circle dots) represented in terms of the fractions of AA'+BB' versus BA'+AB'. The\ndistribution densely populates a linear manifold, providing strong theoretical support for the\nemergence of linear, continuous multistability in CTMN-1. Pentagram markers (S to S )\n1-1 1-7\ncorrespond to the experimentally characterized steady states. v−x, Computationally derived\nreversible displacement rates (v) and polymerization/nicking rates of AB' (w) and AA' (x). Time-\ndependent concentration profiles in b/c, f/g, k/l, and o−s were derived from Supplementary Figs. 5,\n7, 8, and 17, respectively, using the respective calibration curves (Supplementary Fig. 4). Solid\ncurves represent experimentally determined data, and dashed curves indicate computationally\nsimulated results.\niv) Reporting resolvable states as a function of observation window\nThe number of experimentally resolvable steady states depends on (i) the intrinsic\nstability of the system (as captured by the average drift rate under defined biochemical\nconditions and measurement protocols) and (ii) the observation window duration.\nTo quantify this dependence, we estimated the number of resolvable states for\nobservation windows of 24 h and 48 h, using the mean drift rate (0.04 %·h−1) measured\nat 0.2 mM dNTP as a global stability parameter. Under these conditions, our system\nsupports 92 resolvable steady states over a 24 h observation window and 46 resolvable\nstates over 48 h (26 resolvable states over 48 h in the original manuscript).\nA detailed description of how drift rate is mapped to state separability, together with\nthe underlying assumptions and estimation criteria, is provided in the Supplementary\nInformation (p. S11):\n“Quantifying resolvable steady states from drift measurements. To estimate the number of\nexperimentally distinguishable steady states over a given observation window, we followed\nfour steps:\n(i) Mean drift rate\nAs a global metric of steady-state stability, the mean drift rate was estimated from five\nrepresentative steady states (S to S ) spanning the accessible multistable regime under the\n1-1 1-5\noptimal condition of 0.20 mM dNTPs, yielding a value of 0.04 %·h−1.\n(ii) Observation window\nAs the number of experimentally resolvable steady states depends on the observation window\nduration, 24-h and 48-h windows were selected to quantify these states.\n(iii) Accessible compositional range\nThe experimentally accessible compositional range was empirically defined as the AB'+BA'\ncomposition interval spanned by extremal steady state S and S . The lower and upper bounds\n1-1 1-5\nwere measured as 7% and 96%, respectively, yielding a total accessible compositional dynamic\nrange of: 96% − 7% = 89%.\n(iv) Number of resolvable steady states\nThe maximum number of resolvable steady states was calculated by dividing the accessible\ncompositional range by the product of the mean drift rate and the observation window duration.\nFor a 24-h window: 89% ÷ (0.04 %·h−1 × 24 h) = 92 states; for a 48-h window: 89% ÷\n(0.04 %·h−1 × 48 h) = 46 states.”\nv) Contextualizing our results within the field\nAlthough the optimization experiments yield only moderate improvements in the\nnumber of resolvable steady states due to practical constraints imposed by biochemical\ndrift and measurement resolution, our system nonetheless achieves, to the best of our\nknowledge, the highest experimentally validated number of distinct steady states\nreported for any biochemical multistable networks to date (ref. 46: Multistability in\nphysical and living systems: characterization and applications (eds. A.N. Pisarchik &\nA.E. Hramov, Springer International Publishing, Cham; 2022); ref. 51: Science 375,\neabg9765 (2022); and ref. 64: Nat. Commun. 7, 13474 (2016)).\nCrucially, this performance does not originate from incremental parameter tuning.\nInstead, it arises from the mechanistic innovation introduced in this work: a far-from-\nequilibrium, competition-driven architecture that enables a continuous range of tunable\nsteady-state compositions in DNA reaction networks. This architectural principle yields\na new mode of multistable behavior, which we term continuously tunable multistability,\ndistinct from previous designs that typically limited to a small number of discrete fixed\npoints.\nThus, even when constrained by the practical limits of current biochemical\nimplementations, our framework provides a fundamentally new design mechanism for\nrealizing and systematically modulating multistable behavior in DNA-based reaction\nnetworks. We have expanded the discussion in the revised manuscript (p. 21) to\nhighlight the implications for the design of biochemical computation and memory.\nvi) Conceptual strategies to extend the number of experimentally realizable\nsteady states\nIn response to the reviewer’s suggestion regarding potential routes to further enhance\nexperimental scalability, we have included a brief conceptual discussion of possible\nextensions in the revised manuscript (p. 22). These strategies lie beyond the scope of\nthe present study but illustrate how the network architecture could be integrated with\nadditional design features:\nContinuous-flow (CSTR-like) operation: Implementing the network in a\ncontinuous-flow reactor could maintain a constant supply of reagents and remove waste\nproducts, thereby reducing substrate-depletion-driven drift and extending steady-state\nlifetimes relative to the closed-batch operation.\nFeedback-based stabilization: Engineering explicit positive or negative feedback\nmotifs into the network could deepen attraction basins or introduce corrective\nregulatory dynamics, diminishing stochastic drift and potentially increasing the number\nof experimentally accessible states."}, {"round": 2, "reviewer_comments": {"reviewer_1": "THe authors have done convincing and in-depth revisions. I also appreciate clarifications on conceptual novelty! For me the paper is ready for publication.", "reviewer_3": "The authors have addressed most of my concerns, and thanks for the theoretical clarifications regarding scalability (Linear Stability Analysis) and experimental optimizations (dNTP titration, enzyme ratios). However, ​the central claim of “arbitrarily programmable” multistability lacks experimental substantiation, as detailed below:\n1.​Theoretical claims ≠ experimental innovation in DNA programming:\nBreakthroughs in DNA-programmed systems mustprioritize empirical validation. While LSA suggests \"thousands of steady states\" (Response §i), only ​26 resolvable states​ (Response §ii) were experimentally achieved—a marginal figure that contradicts \"arbitrary\" programmability. This gap between theory and practice undermines the innovation’s core promise.\n2.​​\"Arbitrary\" requires experimental scalability:\nProgrammable biological systems are judged by engineered outcomes, not theoretical potential. The reported ​5–92% compositional range​ (Response §ii) with ​high drift (0.07%/h)​​ confirms severe biochemical constraints, not arbitrary control. DNA computing should transcend limitations (e.g., via error correction), not rationalize them.\n3.​Terminological rigor concerns:\nRedefining \"arbitrary\" as \"context-dependent\" (Revised Manuscript p. 11) dilutes scientific precision. Equating \"arbitrary programmability\" with \"​26 states under narrow conditions​\" (Response §ii) misrepresents scalability and sets a misleading precedent for the field."}, "author_response": ""}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have satisfactorily addressed all my concerns, and I therefore believe this manuscript is now publishable."}, "author_response": ""}]

Decision letter