DNA dynamics and computation based on toehold-free strand displacement

We present a simple and effective scheme of a dynamic switch for DNA nanostructures. Under such a framework of toehold-free strand displacement, blocking strands at an excess amount are applied to displace the complementation of specific segments of paired duplexes. The functional mechanism of the scheme is illustrated by modelling the base pairing kinetics of competing strands on a target strand. Simulation reveals the unique properties of toehold-free strand displacement in equilibrium control, which can be leveraged for information processing. Based on the controllable dynamics in the binding of preformed DNA nanostructures, a multi-input-multi-output (MIMO) Boolean function is controlled by the presence of the blockers. In conclusion, we implement two MIMO Boolean functions (one with 4-bit input and 2-bit output, and the other with 16-bit input and 8-bit output) to showcase the controllable dynamics.

is that duplex D1 binds with corresponding blocker to form a blocked-D1(BD1). (b) Reaction (4) is that duplex D1 binds with duplex D2 to form a dimer D1D2. (c) Pre-reaction initial condition is reaction (3) plus reaction (4). (d) Post-reaction initial condition is reaction (3) minus reaction (4). The intermediates are omitted but considered in simulation. Figure 9. Native polyacrylamide gel electrophoresis results of switch from ON to OFF at pre-reaction initial condition in the dual-pair system at 25℃. Lane L: 50 bp DNA ladder. Lane 0×: sample with binding switched ON; lanes n× (n=1, 3, 5, …, 300): samples with binding switched OFF. Percentages above dimer bands indicate dimer yields (mean ± SD, N=3). Figure 16. Schematic of transition events. The blue arrow represents the target strand (D1); the orange arrow represents the substrate (D2 or B); kop, kcl are the kinetic constants for forward and backward reactions of end fraying; kon and koff are those for primary binding. Preformed origami cuboid IV with blockers were excessively added to J-tetromino and L-tetromino, respectively, and pentamer only formed from L-tetromino but not from J-tetromino. (b) Schematic diagram of S-/Ztetrominoes identification. Preformed origami cuboid IV with blockers were excessively added to S-tetromino and Z-tetromino, respectively, and pentamers of a specific shape formed from either tetromino. (c-f) Native agarose gel electrophoresis and TEM results. Left: native agarose gel electrophoresis results; right: TEM images of pentamers. Red arrows point to tetramer bands. Yellow arrows point to pentamer bands. Scale bars: 50 nm. nm. Supplementary Table 1

Supplementary Note 2. Probabilistic modelling of DNA base pairing kinetics
Master equations describe a continuous-time Markov process 1 , which are used here to model the DNA base pairing kinetics of the toehold-free displacement system. Each state of the model describes the probability of observing the system in a specific configuration as a function of time. The dynamic development of the probability of each discrete state can be described by a set of master equations, which requires the definition of the possible states and their transition rates. The master equations can be solved simultaneously as a set of ordinary differential equations (ODEs) when the number of states is of limited size. However, a Monte Carlo approach such as the stochastic simulation algorithm 1, 2 is needed when the number of states is large, which samples solution trajectories to obtain the statistical characteristics of the probability distribution described by the model. Various models based on master equations have been used for simulating DNA hybridization processes involving the kinetics of DNA breathing and the self-assembly of DNA origami. [3][4][5][6] The nature of DNA strand hybridization is discrete and stochastic on a molecular level. The configuration of a double-helix structure can be classified based on the status of individual nucleotides. Transitions between different configurations are dictated by the kinetics of base pairing, which leads to commonly observed processes such as end fraying and branch migration. 7,8 We use these concepts here to construct a probabilistic model of the toehold-free displacement system based on master equations, which allows for better understanding of the kinetic pathways that lead to displacement. Furthermore, the novel toehold-free displacement system can be compared to conventional toehold-mediated displacement systems with different toehold lengths to illustrate the importance of toehold-free displacement for the proposed switch. One of the challenges of modeling the kinetics of DNA self-assembly is the potentially large difference in time scales between the various kinetic events. This causes stochastic methods to sample predominantly fast processes such as end fraying (typically every 10 -6 s) and only rarely the events of interest that change the structural nature of the system such as a new primary binding from solution.
where P  is the probability of finding D1 in configuration  ,    is the propensity factor describing the rate of transition from configuration  to β. Writing Equation (1) for all configurations leads to where P is the state vector consisting of the probabilities of every configuration and ε is a () NN   propensity matrix with elements as follows: if , , where Nα is the number of total configurations. The number of states in the model (Nα) will be large when considering either hybridization with D2 or B at each position (i.e., Nα = 3 16 ), which would complicate the proposed direct solution procedure. Therefore, model reduction is applied. We assume that primary binding of strands B or D2 from solution only occur at the vacant and correct positions of D1 (Supplementary Figure 15.a). Furthermore, at most one strand of D2 and B can be hybridized to the target strand, i.e., the impact of any competition between the same type of strand or displacement from a third one is neglected ( Supplementary Figure 15

END IF END IF END WHILE END IF END WHILE
To construct the propensity matrix, two types of events are considered, i.e., primary binding of D2 or B from solution, edge fraying of a double-strand domain (illustrated in Supplementary Figure 16.). The possible transitions for a given configuration are stored by searching the matrix I with specific conditions. For example, to identify the possible transitions of a D1 configuration ([ld, d1, 0, 0]) by primary binding from a B strand, the searching condition will be: (1) length and position (ld, d1) of D2 of the given configuration; (2) one base pair of D1 binding to blocker (lb=1). To identify the transitions by edge fraying of D2 domain of a given configuration ([ld, d1, lb, b1]), searching conditions will be: (1) where ΔGini is the free energy change for primary binding, ΔGNN is the free energy change for neighboring pairs. In our approach, ΔGbp for base pairing is defined as the averaged value of ΔGNN at a given temperature. Therefore, ΔG for the pairing of strands will be ΔGini+(l-1) ΔGbp. The reaction kinetic constants for pairing steps (kon and kcl) are obtained from the literature 7,8 , and kinetic constants for unpairing steps are calculated based on the ΔG of given events. The propensities for the different events are shown in Supplementary Table 2. The master equations (Eqn. 2) for systems with 4 types of toehold domain were numerically solved by the variable-step, variable-order ODE solver ode15s in Matlab (The MathWorks, Inc.). The probability of observing D1 in a dimer configuration, i.e., with at least one base pair connected to D2, is shown in Supplementary Figure 18. These model simulations predict that the reversibility of the displacement system largely vanishes with an increasing toehold domain length. For instance, a 50% probability of observing the system in a dimer configuration corresponds to a concentration ratio of B to D2 of 1× for the toehold-free strand displacement system, which increases to 100× when there is a toehold of 2-nt. The nearest-neighbor thermodynamic model predicts that the difference of ΔG for the formation of BD1 and D1D2 is a function of the toehold length (Lt) and temperature, which is reflected by the gap of probabilities for the two tested temperatures for different toehold lengths (Lt > 0). While for Lt = 0, an identical ΔG is calculated from our model since D2 and B share the same base sequences, therefore, the same equilibrium probability distributions are obtained regardless of the temperature difference. The dynamic trajectories of the probabilities for the displacement process (Supplementary Figure 19.) demonstrate that the system reaches equilibrium within hundreds of seconds, which is consistent with the experimental observations ( Supplementary Figure 20.). This suggests that the duration of the experiments of one hour was enough to reach the equilibrium state of the displacement system.

Supplementary Note 3. Second-order reaction model of the proof-of-concept single-pair system
The reversible elementary reactions in the proof-of concept displacement system with a single pair of binding partners and their Gibbs free energy change are described by: , (R3) where D1 and D2 are the duplexes. B1 is the blocker fully complementary to D1, B2 is a blocker that is 1-nt shorter than D1. D1D2 is the dimer structure (defined as 'ON' state), B1D1 and B2D1 are the blocked D1 structures (defined as 'OFF' state). The combination of the dimerization (R1) and blocking (R2) reactions on D1 creates a toehold-free strand displacement (TFSD) system, as D2 and B1 are of identical length and have the same base pair sequence. In contrast, when combining the dimerization reaction (R1) with the blocking reaction (R3), a toehold-mediated strand displacement (TMSD) system is created due to the toehold domain on B2D1 resulting from the difference in length between B2 and D1. ΔGD, ΔGB1 and ΔGB2 are the Gibbs free energy change for the dimerization and blocking reactions on a single sticky end. At equilibrium, the concentrations of the species are described by the reaction equilibrium constants as follows, where C~ is the equilibrium concentration of a given species in solution. KD, KB1 and KB2 are the equilibrium constants for dimerization and blocking reactions, respectively, which can be obtained from the corresponding ΔG as follows: where T is the temperature and R is the ideal gas constant. Material balances for the species in the displacement systems complete the model, where C~,0 is the initial concentration of a given species.

Supplementary Note 4. Parameter estimation of ΔG in single-pair system
For the analytical calculation of the concentrations of the target species at a given temperature, we assume that the reaction systems only have one type of blocker at the same time, which is consistent with the experiments, i.e., the concentration of either B1 or B2 is equal to zero. Therefore, CB represents the concentration of either B1 or B2, whichever is nonzero. The equilibrium constants KD, KB1 and KB2 are obtained from equations 7 to 9 for a given temperature T. By substitution of equation 4 into 11, the concentration of D2 as a function of the concentration of D1 is obtained, Substitution of 5 and 6 into 12 and 13 leads to an expression for the concentration of B as follows: From substitution of equation 11, 12 and 13 into 10, expression of 1 D C is derived as Given that for our experimental conditions, the initial concentration of D1 is equal to D2 (i.e., 1 2 D ,0 D ,0 CC  ). Therefore, by substituting equation 14 and 15 into 16, an implicit expression in which the concentration of D1 is the only unknown can be obtained, The coefficients of equation 16 compared to the standard cubic equation The equilibrium state of the displacement system at a given temperature is only governed by the initial concentrations of D1, D2, B1 and B2, and the ΔG for the elementary reactions. We aim to fit the experimental data involving the fractions of different species (i.e., D1 D2, B1D1, B2D1 and D2, fractions of given species in S1) at different initial blocker concentrations to the model presented in this section by estimating the values of ΔGD, ΔGB1 and ΔGB2 from the following optimization problem (P1) based on maximum likelihood estimation: where Nr is the total number of tested ratios of the blocker concentration to the D2 concentration, and M is the total number of measurements for each experiment. is the measured fraction for a given species from experiments, and ~,i F is the fraction predicted from the model at the same conditions. 2 ,i  is the variance of experimental data. The value of the variance is replaced with an averaged variance of the measurements of other species from the same experiment in case the variance is equal to zero. The lower bound for ΔG of -30kcal/mol is chosen sufficiently low compared to the value predicted from the nearest-neighbor thermodynamics model 9 (around -22kcal/mol) to allow for realistic solutions while avoiding potential numerical complications resulting from an unrestricted search space. Optimization problem P1 was implemented in Matlab and was solved with the fmincon solver based on the experimental results from two temperatures, i.e., 25℃ and 40℃. The experimental and simulation results are illustrated in Supplementary Figure 21. The model simulations consistently describe the decreasing trends of dimer fractions with the increasing blocker concentrations. The model simulations well fit the fractions of D1D2, B1D1 and D2 in the TFSD system, while lower fractions of D1D2 and higher fractions of D2 from TMSD system when N < 25 are obtained from the reaction model. The TFSD and TMSD share same strands with different length of blocker strand, so the dimer fractions from experiments without blocker should be identical. But an increment of fractions of D1D2 and a decrement of fractions of D2 at N = 0 from TMSD compared to the TFSD were obtained from the experiment results, which is unexpected. Such differences of fractions could attribute to unreplicated experiment of TMSD. Our reaction model only involves ΔG and concentrations of species as variables, therefore, it could not capture the differences from unexpected effect. The estimated values for ΔGD, ΔGB1, ΔGB2 at different temperatures are listed in Supplementary Table 3. The estimated ΔGD, ΔGB1 and ΔGB2 are higher than those predicted from the nearest-neighbor thermodynamics model, which could point to a reduced entropy change during the pairing reaction of binding partner of D1 compared to a standard duplex formation. Focusing on the ΔGD, ΔGB1 from TFSD, the pairing of D1 by D2 appears more stable than pairing with a blocker at 25℃ and 40℃, as the estimated ΔGD is lower than ΔGB1. The reason why dimer is more stable could be some internal interactions (e.g., Van der Waals forces) between two duplexes. The effect of temperature on the dimerization process is not substantial based on our model simulations, as the experimental fractions without blockers at 25℃ and 40℃ do not differ substantially. ΔGB1 decreases with increasing temperature suggesting that B1D1 is more stable at a higher tested temperature, which could be resulted from some unexpected effects. For instance, a higher temperature will promote the entanglement of the poly-T domain of binding partner of B1D1 with the duplex, which could stabilize the blocked duplex. Such entanglement is difficult to occur in a dimer structure as the binding partner connected to two duplexes is less flexible. Therefore, a low temperature favors the 'ON' state.

Supplementary Note 5. Second-order reaction model of the proof-of-concept dual-pair system
In the dual-pair system, two types of blockers are introduced to complement the dual binding partners. The binding partners and blockers all share the same length; therefore, the displacement process is toehold-free. We assume that the paring of any pair of binding partners involves identical free energies (ΔGD), and blocking of any binding partner involves identical ΔGB. The dual-pair system can be described by four reaction equations with corresponding free-energy changes as follows: , (R4) , (R7) where DI and DII are the duplexes, BI and BII are the blockers fully complementary to the corresponding sticky ends on DI. DIDII is the dimer structure (defined as 'ON' state), and DIBIBII is the fully blocked structure ('OFF' state). DIBIDII and DIBIIDII are the intermediate structures with one of the binding partners blocked and the other connected to DII. The equilibrium concentrations of the reactions are described by: where C~ is the equilibrium concentration of a given species. KD and KB are the equilibrium constants described by:  I  I  I II  I I II  I II II  I I II  II  I II  I I II  I II II where C~,0 is the initial concentration of a given species.

Supplementary Note 6. Parameter estimation of ΔG in dual-pair system
The model has an explicit solution for the concentration of DI, which can be derived as follows. Since BI and BII share identical total concentration ( B,0 C ), it can be shown that (27) We assume that the concentrations of the intermediate species (i.e., DIBIDII and DIBIIDII) are negligibly small compared to the concentrations of the dimer structure and the fully 51 blocked species. Under this assumption, the model can be reduced to two reaction equations, R4 and R7, and only the equilibrium (17 and 20) will be involved. The mass balances are reduced as below (30) By combining equations 28 to 30 with the experiment condition that DI and DII share identical total concentration, the material balance of DI can be written as: (33) Note that CB should be nonnegative value. Therefore, the only feasible solution of quadratic equation 33 is given by: An approximation to the square root by continued fraction is applied: We use the second iteration of 35 to obtain: The approximation will result in a 5% difference between the LHS and RHS of 36 by substituting the obtained I D C from equation 38. The concentration of the blockers as a function of concentration of DI is given by   I I II  I I II  I II  I I II  II where α is a weighing factor for improving the quality of the fit at higher blocker concentrations, as the experimental results of dimer fractions are around 100% at low concentrations of blockers (N < 25). We used two sets of α for two temperatures, where BDI represents fully blocked species. Noting that equation 40 and 41 are under the condition that all the blockers share identical total concentration. In case of the total concentration CB,0 < 1M, increasing binding partners will reduce the magnitude of B p n C . Therefore, when the total concentrations of blocker and DII are close (e.g., N < 5 in the dual-pair experiment), is substantially larger than B p n C , which is the reason why DIDII is the dominant structure when N < 5. The threshold of blocker concentration is the concentration which shifts the equilibrium composition from dimer structure to blocked species. Such shifting of equilibrium composition needs to be identified by PAGE separation. The estimated values for the free energy are listed in Supplementary Table 4. Opposite to the findings from the single-pair system, ΔGB is smaller than ΔGD, which is consistent with a blocked state that is more stable than the dimer. This difference could be resulted from an entropy penalty from the internal loop formation of the dimer with dual binding partners.
To conclude, increasing binding partners on the duplexes will raise the requirement of concentration of blockers to obtain a switch between ON and OFF states. OFF state will be more stable with the increased pairs of binding partners. In the case of sufficient concentration of blocker, switch from ON to OFF state will be benefit from multiple binding partners.

Supplementary Note 8. Reactions for switch from ON state to OFF state
The preformed origami rectangles I and II without purification were directly mixed and subjected to an isothermal annealing program (40℃ for 40 hours) before native agarose gel electrophoresis purification. The purified dimers (ON state) were turned to be switched to monomers (OFF state) with blocker set N' at gradient concentrations (0×, 1×, 3×, 5×, 25×, 100× and 300×) in an isothermal incubation (48.6℃ for 40 hours). Samples after displacement reaction were subjected to native agarose gel electrophoresis.

Supplementary Note 10. Reactions for the 4-bit input/2-bit output system
The preformed origami cuboid I (with sticky faces N and M) and origami cuboid II (with sticky faces N* and M*) after purification were mixed with corresponding blocker sets N' and M' at an excess concentration in 0.5× TAE buffer supplemented with 35 mM MgCl2. Samples were subjected to an isothermal incubation (42℃ for 14 hours) before native agarose gel electrophoresis.

Supplementary Note 11. Definition and the MIMO Boolean functions of a dualunit system of origami cuboids
The Boolean function contains inputs (system state and controller) and output (new system state). Binary variables (0/1) of the system state are sticky faces, and that of the controller are blockers. We use 0 to represent sticky face not being blocked or blockers not being added. And we use 1 to represent sticky face being blocked or blockers being added with an excess amount. Boolean function formulation may describe iterative computations. In the manuscript's Fig. 3(a) the Boolean function of 2-bit input and 1-bit output is described. Let sk = sk(1) be the system state at the k-th iteration, where sk(1) represents whether N is blocked by N'. Let ck = ck(1) be the controller at the k-th iteration, where ck(1) represents whether blocker N' is added at an excess amount. Then sk+1 = sk+1 (1) represents the new system state after the k-th iteration. We have sk+1 (1) = f(sk (1), ck (1)), where function f is defined by the truth table in Supplementary Table 5. In the manuscript's Fig. 3(b), the Boolean function of 4-bit input and 2-bit output are described. Let sk = (sk(1), sk(2)) be the system state at the k-th iteration, where sk(1) and sk (2) represent whether N and M are blocked by N' and M', respectively. Let ck=(ck (1), ck (2)) be the controller at the k-th iteration, where ck(1) and ck(2) represent whether the blockers N' and M' are added at an excess amount, respectively. Then sk+1 = (sk+1(1), sk+1(2)) represent the system state after the k-th iteration. Since the dynamics at the first 1-bit system state sk(1) and the second 1-bit system state sk (2) are independently controlled by ck (1) and ck (2), respectively. We have sk+1 = (f1(sk (1), ck (1)), f2(sk (2), ck (2))), where function f1 and f2 are defined by the truth table in Supplementary Table 6. We can see that f1 and f2 are the same as the aforementioned f. So we have sk+1 = (f(sk (1), ck (1)), f(sk (2), ck (2))). This example shows that we may use the Boolean function formulation presented in this work to describe internal structures of the function, such as iterations and compositions.

Supplementary Note 12. Reactions for tetromino assembly
The preformed origami cuboids I, II, III and IV after purification were mixed with certain combinations of blocker sets at an excess amount in 1× TAE buffer supplemented with 35 mM MgCl2. Samples were subjected to an isothermal incubation (42℃ for 14 hours) before native agarose gel electrophoresis.

Supplementary Note 13. The MIMO Boolean functions of a quadruple-unit system of origami cuboids
The tetromino assembly is a Boolean function of 16-bit input and 8-bit output. Among the 16-bit input in this quadruple-unit system, the first 8 bits represent whether sticky faces are being blocked, and the rest 8 bits represent whether the blockers are added at an excess amount. Then the 8-bit output define the new system state. This Boolean function as f3, relative to tetromino shapes, is shown in the truth table at Supplementary  Table 8. Let sk = sk(i), (i = 1, 2, 3, 4, 5, 6, 7 and 8), where sk(i) represents whether sticky faces (A, B, C, D, E, F, G and H) are blocked by the blockers (A', B', C', D', E', F', G' and H'), respectively. Then sk+1 = (sk+1(i), i=1, …, 8) represent the new system state after the k-th iteration. So we have sk+1 = f3(sk, ck) where for each i = 1, …, 8 we have sk+1(i) = (f (sk(i), ck(i))) in which function f is defined by the truth table in Supplementary Table 8. The graycolored rows in Supplementary Table 8 represent ones to be chosen as experimental groups shown in manuscript's Fig. 5.