Hierarchical control of enzymatic actuators using DNA-based switchable memories

Inspired by signaling networks in living cells, DNA-based programming aims for the engineering of biochemical networks capable of advanced regulatory and computational functions under controlled cell-free conditions. While regulatory circuits in cells control downstream processes through hierarchical layers of signal processing, coupling of enzymatically driven DNA-based networks to downstream processes has rarely been reported. Here, we expand the scope of molecular programming by engineering hierarchical control of enzymatic actuators using feedback-controlled DNA-circuits capable of advanced regulatory dynamics. We developed a translator module that converts signaling molecules from the upstream network to unique DNA strands driving downstream actuators with minimal retroactivity and support these findings with a detailed computational analysis. We show our modular approach by coupling of a previously engineered switchable memories circuit to downstream actuators based on β-lactamase and luciferase. To the best of our knowledge, our work demonstrates one of the most advanced DNA-based circuits regarding complexity and versatility.

and FAM fluorophores which is 0 in the absence of template's input primer and 1 at the steady-state value of primer β and α respectively. The results confirm that the bistable switch with asymmetric amounts of autocatalytic templates and inhibitory templates can keep its steady-states for an extended time period and, therefore, is indeed bistable until accumulation of waste, limiting amount of dNTP's and enzyme activity increasingly dominate the dynamics of the network.

Supplementary Figure 5. Coupling of the translator module to a two-input bistable switch. a,
Schematic illustration of the system, in which the translator module is coupled to β or α of the PENbased bistable switch. The core of the bistable switch consists of four templates including the autocatalytic templates αtoα and βtoβ and the inhibitory templates αtoiβ and βtoiα. The network switches between states upon injection of γ and δ which are received by templates γtoα and δtoβ.
The dynamics of the bistable switch are followed via N-quenching using templates βtoiα and αtoiβ which are 3'-end labeled with a DY530 and FAM fluorophore respectively. b, Results of the experiments in which a concentration range of translator template was coupled to β or α. The sequence of the translator template was adapted for coupling to β enabling the translator module to receive β as input primer while producing σ as output strand. Experiments were carried out as described in the Methods (PEN-based experiments) using 20 nM βtoiα, 15 nM αtoiβ, 24 nM βtoβ, 10 nM αtoα, γtoα and δtoβ, 10 U/mL Bst 2.0 warmstart DNA polymerase, 10 U/mL Nt. bstNBI and 200 nM ttRecJ. The switch was initiated with 1 nM of α when the network was switched from the α-to the β-state, while the network was initiated with 1 nM of β when switched from the β-to the αstate. The dotted lines indicate the time point at which 30 nM δ (from αto β-state) and γ (from β-to α-state) was injected. The charge level is the normalized fluorescence of the signal of DY530 and FAM fluorophores which is 0 in the absence of template's input primer and 1 at the steady-state value of primer β and α respectively. The results of the translator coupled to α or β show the production rate of σ follows the dynamics of the switch instantaneously. Furthermore, the effect of retroactivity from the translator coupled to α or β was analyzed by plotting the data in the phase plane of αtoiβ and βtoiα (Fig. 4). c, Results of simulations performed using the heuristic model

Determination of the thermodynamic dissociation constant of DNA hybridization
The thermodynamic dissociation constants of DNA hybridization were determined from the melting curves ( Supplementary Fig. 18 and Supplementary Table 2), obtained using JASCO V-650 spectrophotometer and a 1 cm path length cuvette with a volume of 200 μL. UV absorbance of all possible partial duplexes in the DNA-based network was measured, at a wavelength of 260 nm, as a function of temperature. A temperature gradient of 1°C min -1 was used, since melting and cooling profiles were significantly similar at this gradient. The melting curves were converted from absorbance to fraction associated (DNA hybrid) using two baselines ( Supplementary Fig. 18a). This curve was used for non-linear least squares analysis using the following equation to obtain the enthalpy and entropy: 4 where θ is the fraction of partially duplex, T is the temperature in Kelvin, R is the gas constant (kcal mol -1 K -1 ), and C 0 is the initial concentration of the duplex DNA (M) divided by the molarity of water (M). Nonlinear least square optimization was performed using the Matlab routine lsqnonlin which uses a subspace trust-region method based on the interior-reflective Newton method. The enthalpy (kcal mol -1 K -1 ) and entropy (kcal mol -1 ) obtained from this fitting were used to determine the standard Gibb's free energy (kcal mol -1 ) of DNA hybridization (from which the thermodynamic dissociation constant can be calculated) at a temperature of 42 °C: The thermodynamic binding constants of the measured hybrids are shown in Supplementary Table   2.  Table 2 and vide supra). Finally, the mean of the multiple experiments was determined for further analysis. To obtain the second-order rate constant (k a ) non-linear least squares multiple-curve fitting was performed (black) using the Matlab routine lsqnonlin with a subspace trust-region method based on the interior-reflective Newton method. The hybridization association equilibrium constant (K β ) was determined experimentally (vide supra) and was fixed during the non-linear least-square analysis (Supplementary Fig. 19).

Determination of the kinetics of exonuclease
To When this was observed the initial data was removed, and an estimate of these datapoints was made by extrapolation, i.e. by performing linear fitting to the first part of the raw data curve an estimate could be made of the RFU of the first datapoints. The RFU at the timepoint at which its values became constant, corresponding to the timepoint at which all non-protected primer was degraded by exonuclease, was subtracted from the datacurve, so that at t = t end an RFU of zero was obtained. Then, the RFU was converted to the concentration of ssDNA by dividing the RFU with a factor obtained by dividing the initial RFU by the initial concentration of degradable ssDNA. Using the datacurve obtained from the experiment without competitor ( Supplementary Fig. 20a)

Determination of the kinetics of the molecular beacon
To characterize the kinetics of the molecular beacon, experiments in triplo were performed (Supplementary Fig. 21 and Supplementary Table 2) in which the molecular beacon in master mix without enzymes was added to a cuvette which was preheated at 42 ᵒC in a fluorescence spectrophotometer (Carry Eclipse equipped with a Varian Peltier Multicell Holder and a Cary Temperature Controller). Subsequently, the fluorescence of the molecular beacon was measured to obtain the baseline fluorescence. Thereafter, DNA strand σ was added giving a volume of 120 µL after which the mixture was suspended and measurement started at 42 ᵒC. Data handling was done for each single curve in the same way. First, the baseline fluorescence was determined by taking the mean of a one minute measurement of the fluorescence from the molecular beacon. After this, the fluorescence from the measurement was subtracted by this baseline value. Subsequently, the missing points due to suspending (10-15 seconds) were estimated by extrapolation. Thereafter, the fluorescence was converted to concentration of opened beacon by assuming the reaction was equilibrated in 18 minutes. An ODE model based on the bimolecular reaction model of DNA strand displacement was developed to describe the kinetics of this step. 5 To obtain the second-order rate constant (k rep ) of the toehold mediated strand displacement reaction non-linear least-square optimization of the experimental kinetic traces (depicted in Supplementary Fig. 21b) to the mathematical model was performed. The Matlab routine lsqnonlin with a subspace trust-region method based on the interior-reflective Newton method was applied, yielding a k rep of 5.3 x 10 7 M -1 min -1 .

Mathematical model of Bistable Switch: A heuristic model
Following [1], we assume that the activation of a primer y on a template T and an activator x x + y + T can be described by the Michaelis-Menten derivation: where V is proportional to the total concentration of template. Furthermore, the second step with reaction rate k 2 comprises two enzymatic reactions including polymerization and nicking. Therefore, we assume k 2 << k −1 and K becomes the dissociation constant of x on T which can be experimentally determined (vide inf ra). Here the square bracket notation [ ] is used to denote the concentration of the corresponding species. (We remark that t is an independent variable that only resembles time.) If iy is a competitive inhibitor acting on y, then where constant λ y roughly depends on the dissociation constant of iy and the dissociation constant of x.
The degradation of a primer x, x + E xE E will be described by a competitive Michaelis-Menten mechanism: with the Michaelis-Menten constant K M,exo as experimentally determined (videinf ra). The parameter V exo was estimated to be lower as experimentally determined in isolation, since sequestration of exonuclease by the templates is present. We manually adapted this value.

Switch with translator
We describe the model of the translator for the case that it is coupled to the α-side of the switch.
Coupling to the β-side of the switch is done analogously. The chemical scheme of the translator is as follows: αT L σ .
∅ and binds to a reporter We model the activation of σ via α again via a Michaelis-Menten mechanism: The activation of σ changes the concentration of α by: where the first term at the right-hand-side represents the amount of α that is used by the translator for the production of σ, and the second term accounts for the reproduction of α due to the dissociation of α from the nicked translator module. The rate of reproduction is given by θ representing the concentration of nicked translator module, linearly scaled with the dissociation rate constant of α. Furthermore, the ordinary differential equation of θ is given by: where k α models the dissociaton rate constant of α and ρ models the fraction of translator module being in the nicked state. The dissociation rate constant of α is determined by the equilibrium dissociation constant K Lα as the association rate constant is invariable for primers with lengths exceeding five bases. 4 The constant depends on the timescale of nicking the duplex relative to the timescale of the polymerase strand displacement reaction. In the extreme case of ρ = 1 the equilibrium of the two states of the translator module is shifted to the nicked state (minimal inherent retroactivity) and, hence, the amount of α reproduced depends on K Lα . In the other extreme case ρ = 0 no α is reproduced (maximal inherent retroactivity) independent on K Lα . Note that for being at steady-state it is required that the rate of polymerase is identical to the rate of nickase, which implies that ρ depends, to a large extend, on the reaction rates of polymerase and nickase.
The binding of σ to the reporter M B is assumed to be governed by mass-action-kinetics and can be described by the system of ODEs Letting we ensure that [iβ] > 0 for [α] > 0. Likewise, for SS2 we find such that for we have [iα] > 0 if [β] > 0. Inequalities: V exo > V j , j = α, β, iα, iβ ensure that solutions are bounded in forward time. Furthermore, it is interesting to note that imply that the trivial steady-state solution is unstable. (This can easily be verified by determining the Jacobian matrix around this steady-state).
We remark that there also exists an other steady-state solution with [α] > 0 and [β] > 0. This steady-sate can be shown to be unstable when λ α and λ β are sufficiently large ( Supplementary Fig. 27).

Numerical simulations
Given the constraints on the parameters that are obtained with the steady-state analysis we choose the remaining parameters such that there is a good qualitative correspondence between numerical simulations and experiments. Initial conditions are chosen as