Visual synchronization of two 3-variable Lotka-Volterra oscillators based on DNA strand displacement.

DNA strand displacement as a theoretic foundation is helpful in constructing reaction networks and DNA circuits. Research on chemical kinetics is significant to exploit the inherent potential property of biomolecular systems. In this study, we investigated two typical reactions and designed DNA strands with a fluorophore and dark quencher for reaction networks using a 3-variable Lotka-Volterra oscillator system as an example to show the convenience of and superiority for observation of dynamic trajectory using our design, and took advantage of the synchronization reaction module to synchronize two 3-variable Lotka-Volterra oscillators. The classical theory of chemical reaction networks can be used to represent biological processes for mathematical modeling. We used this method to simulate the nonlinear kinetics of a 3-variable Lotka-Volterra oscillator system. This journal is © The Royal Society of Chemistry.

Introduction

Nonlinear kinetics exists widely in mass chemical reaction equations. Chemical reactions with DNA strand displacement (DSD) also have the property of being nonlinear and cascading. We can use the property of cascading to construct complex chemical reaction networks (CRNs) as well as to compile chemical equations into circuit systems, for example, in pulse counter and incrementer state machines.[1-3]

With the aid of DNA strand displacement, massive DNA circuits, including digital and analog circuits, have been designed and validated by experiments.[4-13] The outputs of the DNA circuit are represented by the concentrations of the DNA strands, but the detection of concentration is difficult or even impossible in some cases. For example, the structure of a DNA polymer is massive, complex, and unstable, such that detection of the concentration of the DNA polymer is impossible.[14] If the DNA strand is designed with a fluorophore and dark quencher, observation of the output is convenient and fast. Luminous intensity is enhanced with an increase in the concentration of DNA strands.[15] Thus, it is helpful to measure the concentration of DNA strands by fluorescence intensity. This property is particularly useful for digital circuits because these circuits are encoded with binary bits, where 2 possible values of each bit are “0” and “1”. If the luminous intensity is strong enough, the value of the output is “1”, otherwise the value of the output is “0”. For an analog circuit, the output is represented qualitatively by luminous intensity; for polymerization of a hybridization chain reaction,[14] characterization of the hyperbranched chain reaction can be recognized by fluorescence intensity.

Synchronization is an important dynamical behavior and a general and ubiquitous phenomenon that has been investigated in a variety of fields such as complex networks, chemical reaction networks,[16-18] secure communication, and genetic regulatory networks. Recently, it was found that synchronization is beneficial to DNA nanotechnology; for example, synchronization is applied in self-assembly,[19] gene expression[20] and DNA replication.[21] Synchronized systems can enjoy common dynamical behavior and isochronous signal transition or transmission. Thus, synchronization has extensive application potential in digital and analog DNA circuits.

As compared to previous studies, we have proposed catalysis, degradation and synchronization reaction modules with a fluorophore and dark quencher, and identified different DNA strands with different colors of fluorophores, where the synchronization reaction module is a new reaction module. Then, we compiled these reaction modules into the 3-variable Lotka–Volterra oscillator and realized synchronization of two 3-variable Lotka–Volterra oscillator systems. From the simulations of these reaction modules, it was demonstrated that the DNA implementation of CRN can imitate the dynamical behavior of an ideal formal reaction, and there is a relationship between the reaction rates of DNA CRNs, intermediate CRNs and ideal formal CRNs.

We used visual DSD software to design the DNA strand and test the feasibility of the chemical reaction with DNA strand displacement, and utilized chemical reaction network (CRN) software to numerically simulate the dynamic behavior of the 3-variable Lotka–Volterra oscillator.

Recent work

DNA strand displacement with fluorophore and dark quencher

We presented a reversible, bimolecular chemical reaction, as shown in Fig. 1, where A1 and A2 are reactants and A3 and A4 are products; k1 denotes the forward reaction rate, and k2 denotes the backward reaction rate. The reaction begins with an invader strand A1 binding to the target strand A2 based on the Watson-Crick model.[22] The position, number, and length determine the feasibility, reversibility, and reaction rate of DNA strand displacement, respectively. The reaction rate first increases and then decreases with an increase in the length of the toehold,[23] such that we can adjust the reaction rate by changing the length of the toehold. Reversible and irreversible reactions are 2 types of DNA strand displacement reactions.

Fig. 1

DNA strand displacement with the fluorophore and quencher. The arrow in a DNA strand indicates the 3′ end. The domain a2 is 12 nucleotides, and the toehold domain is usually less than 10 nucleotides, so the toehold binding reaction is reversible.

The fluorophore and quencher exist at the end of the strand. When the fluorophore and quencher appear in pairs, the fluorophore strand does not fluoresce, but when the fluorophore strand releases from the dark quencher strand, it can be detected by fluorescence intensity.

DNA strand displacement with hairpin structure

Hairpin[24-26] is a special structure for single strand DNA (ssDNA), which can protect some domains. As shown in Fig. 2, domains a2, b1, and b2 are protected before the hairpin is open. When an invader DNA strand X reacts with H, the hairpin is open and then, domains a2, b1, and b2 can be used for further reactions.

Fig. 2

DNA strand displacement with a hairpin structure. The toehold b1 is as long as the other toeholds. The hairpin cannot open spontaneously as a result of the dissociation of hybridization at the end of the branch migration, the loop portion a2 – b1 – b2 of the hairpin is protected if it does not open.

Results

CRN modules

In DNA strand displacement as a type of chemical reaction, there is a corresponding relationship between DNA CRNs and the ideal formal CRNs, where there are intermediate CRNs is the transition between DNAs and ideal formal CRNs. In order to endue CRNs modularity, we propose three types of reaction modules by DNA CRNs and ideal CRNs, where [X]0, [XM]0 and [X1]0 represent the initial concentration of species X, XM and X1 in DNA CRNs, intermediate CRNs and ideal CRNs, respectively, and satisfied [X]0 = [XM]0 = [X1]0. [Y]0, [YM]0 and [Y1]0 represent the initial concentration of species Y, YM and Y1 in DNA CRNs, intermediate CRNs and ideal CRNs, respectively, and satisfied [Y]0 = [YM]0 = [Y1]0.

In the catalysis reaction module (1), as shown in Fig. 3(b), is produced by XM in (5), and the increase in will catalyze (6), while the producer XM can catalyze the production of . According to (1), d[X]/dt = −qC[X] and d[XM]/dt = −θ1[XM] by (5); if DSD (1) and (2) can approximate to the intermediate (5), q = θ1/C is obtained. According to (3), d[X′]/dt = −qC[X′], while based on (6). If DSD (3) and (4) can approximate to the intermediate (6), q = θ/C. d[X]/dt = −θ2[X] and d[XM]/dt = −qC [XM], according to (4) and (6); if (4) is equivalent to (6), q = θ2/C. Similarly, we can obtain k1 = θ1 based on (5), (6) and (7). Initial concentrations of the auxiliary species and are set to C, where [X]0 ≪ C, and [X]0, [Y]0 are the initial concentrations of X and Y. As shown in Fig. 3(c), the evolutions of X in DNA CRNs, intermediate CRNs and ideal CRNs are approximately coincident. This shows that DNA CRNs, intermediate CRNs and ideal formal CRNs in the catalysis reaction module (1) are equivalent, but there are nuances between DNA CRNs, intermediate CRNs and ideal formal CRNs.

Fig. 3

Catalysis reaction module (1) X → 2X. DNA reactions are listed in (a); (b) exhibits the DNA CRNs and appropriate ideal CRNs, where intermediate CRNs is the transition between DNAs and ideal formal CRNs; (c) shows the result of the variable X in DNA CRNs, media CRNs and ideal formal CRNs respectively, where C = 104 nM and [X]0 = [XM]0 = [XI]0 = 5 nM; q = 10−8 nM s−1, q = 10−7 nM s−1 and q = 10−5 nM s−1. DNA implementation of formal species X′ is gained by species X. Species X reacts with to displace in (1), releases two X′ in (2). In (3), species X′ reacts with hairpin to produce , and and releases X in (4) on reaction with . This reaction replaces the ssDNA strand without the fluorophore to the ssDNA strand with the fluorophore, as adapted from.[4]q represents a maximum strand displacement rate constant under the assumption of equal binding strength for all full toeholds; when q and k are satisfied q,q ≪ q, k ≪ qC; [X]0 and C are satisfied [X]0 ≪ C, as adapted from.[1]

In catalysis reaction module (2), as shown in Fig. 4(b), is produced by XM and YM in (13), and the increase in will catalyze (14), while the producer YM can catalyze the production of and consumption of XM. According to ref. 1, we have. q = θ3. Similar to catalysis reaction module (1), q = θ4/C, q = θ/C and k2 = θ3. Initial concentrations of the auxiliary species A, B, H, P, Sp, and Sq are set to C, where [X]0, [Y]0 ≪ C. As shown in Fig. 4(c), the evolutions of Y in DNA CRNs, intermediate CRNs and ideal CRNs are approximately coincident. This indicates that DNA CRNs, intermediate CRNs and ideal CRNs in catalysis reaction module (2) are equivalent, but there are nuances between DNA CRNs, intermediate CRNs and ideal CRNs.

Fig. 4

Catalysis reaction module (2) X + Y → 2Y. DNA reactions are listed in (a) and (b) exhibits the DNA CRNs and corresponding ideal formal CRNs, where intermediate CRNs is the transition between DNAs and ideal CRNs; (c) shows the result of the variable Y in DNA CRNs, media CRNs and ideal formal CRNs, respectively, where C = 104 nM, [X]0 = [XM]0 = [XI]0 = 5 nM, and [Y]0 = [YM]0 = [YI]0 = 1 nM; q = 10−5 nM s−1, q = 10−6 nM s−1 and q = 10−3 nM s−1. DNA implementation of formal species is gained by species XM and YM in intermediate CRN (13). Species X reacts with A to displace T. In (9), C is produced when Y irreversibly displaces T, and H releases 2 Y′ in reaction (10), as adapted from.[1] In (11), species Y′ reacts with hairpin P to produce S, and S releases X in (12) on reacting with Sp, as adapted from.[4]

In the synchronization reaction module, it is assumed that X and Ma attain instant equilibrium through (18), with [X]/[Ma] = q/q; when XM and Ma attain instant equilibrium through (24), we obtain [XM]/[Ma] = θ5/θ. According to (19) and (25), we have d[Ma]/dt = −Cqs[Ma] = −θ6[Ma]; as a result, qs = θ6/C. Taking (25) and (28) into account, we obtain d[YM]/dt = θ6[Ma] and d[YI]/dt = k4[XI], respectively; if intermediate reaction (25) can be approximated to the ideal formal reaction (28), θ6[Ma] = k4[XI]. Because [X]/[Ma] = q/q when X and Ma attain instant equilibrium and [XI] = [X]+[Ma], we can obtain k4 = qθ6/(q + q) = qsqC/(q + q). The initial concentrations of the auxiliary species Na, , Nb, Nc, Nd, , Ne, and Nf are set to C, where [X]0,[Y]0 ≪ C. Fig. 6(c) shows that X and Y can share the same dynamic behavior in DNA CRNs, intermediate CRNs and ideal formal CRNs, but their synchronous values are different. The reasons for this phenomenon are as follows:

Fig. 6

Synchronization reaction module. (a) A list of DNA reactions; (b) exhibits the DNA CRNs and appropriate ideal formal CRNs, where intermediate CRNs is the transition between DNAs and ideal formal CRNs; (c) represents the process of synchronization between X and Y in DNA CRNs, intermediate CRNs and ideal formal CRNs, respectively, where C = 104 nM, [X]0 = [XM]0 = [XI]0 = 4 nM, and [Y]0 = [YM]0 = [YI]0 = 8 nM; qs = 10−6 nM s−1, q = 3 × 10−4 nM s−1 and q = 10−3 nM s−1. In (18), X react with Na producing Ma reversely; the producer Ma reacts with Nb to produce Mb in (19), and Y is produced by Mb and Nc, which is equivalent to YI being produced by XI in ideal CRNs (28). DSD implementation of ideal reaction (29) is similar to (28).

I. In DNA CRNs

When DNA CRNs approach reaction equilibrium, the following can be obtained:

Considering (30) and (31), we can obtain the synchronous value:

II. In intermediate CRNs

When intermediate CRNs approach reaction equilibrium, we can obtain

III. In ideal CRNs

Eqn (36) is satisfied when X1 and Y1 approach instant equilibrium through (27) and (28)

Therefore, the synchronous value of X1 and Y1 are obtained as follows:

3-Variable Lotka–Volterra oscillator

The Lotka–Volterra model is one of the earliest sound mathematical principles[27-29] to describe population dynamics, and the Lotka–Volterra oscillator is known to be an ecological food chain oscillator. A 2-variable Lotka–Volterra oscillator can be described by ideal formal CRNs as follows:[1]where X1 and X2 are reactants; r1, r2, and r3 represent the reaction rate of reaction (39a–c), respectively. Reaction (39a) leads to an increase in X1, and X1 can catalyze reaction (39a).

The product of X1 will increase the product X2, and X2 can catalyze reaction (39b). When the concentration of X1 increases to a certain degree, an increase in X2 will lead to consumption of X1 and the decrease in X1 will slow down reaction (39b), which can lead to a decrease in X2. Therefore, in the cycle, the concentrations of X1 and X2 will fluctuate over time.

According to the basic Lotka–Volterra oscillator, a 3-variable Lotka–Volterra oscillator model can be described by ideal CRNs as follows:

CRNs (40) can be represented by the following set of ODEs:

Fig. 8(a) shows the evolution of Lyapunov exponents λ (i = 1,2,3) of ideal formal CRNs (40) with the chemical reaction rate parameter, r1. It is known that the kinetic characteristics of the system can be characterized by the Lyapunov exponent, if exist in the range in which one Lyapunov exponent is positive and another Lyapunov exponent is negative, indicating that the system is at a chaotic state. Furthermore λ1 < 0, λ2 > 0 and λ3 > 0 when the parameter is taken as r1 = 5.5/300 s−1, and the phase map of the state variables X1, X2 and X3 is shown in Fig. 8(b).

Fig. 8

(a) The evolution of Lyapunov exponents λ (i = 1, 2, 3) with the chemical reaction rate k1. (b) The chaotic attractor in the (X1, X2, X3) phase space and its two-dimension projection.

Synchronization of two 3-variable Lotka–Volterra oscillator systems

Considering another 3-variable Lotka–Volterra oscillator systemwhere the initial concentrations of the above two systems are different. In this study, we used four types of fluorophores as reporters: ssDNA X1 and Y1 are identified by the TYE563 fluorophore (blue ball) with λEx = 549 nm, λEm = 563 nm; ssDNA X2 and Y2 are identified by the ROX fluorophore (red ball) with λEx = 588 nm, λEm = 608 nm; ssDNA X3 and Y3 are identified by FAM fluorophore (green ball) with λEx = 495 nm and λEm = 520 nm; and the dark quencher is Lowa Black RQ.[11]

In order to realize internal synchronization between two chemical reaction networks, we added coupling terms as follows:

ODEs of CRNs (45) and (46) can be described as follows:

The error state was defined as e = y − x (i = 1, 2, 3), and the dynamics of error can be described as follows:

In simulations, the two 3-variable Lotka–Volterra oscillator systems with different initial values, and the evolution of the maximum Lyapunov exponent of eqn (53) with the chemical reaction rate r6 is given in Fig. 9. When the maximal Lyapunov exponent λ is negative, the synchronization of the two 3-variable Lotka–Volterra oscillator systems can be obtained.[30]

Fig. 9

The evolution of maximum Lyapunov exponent with chemical reaction rate r6.

CRNs (45) and (46) can be described by the abovementioned DNA reaction modules as follows:

Visual DSD software can be used to realize the approximations of the CRNs (45) and (46).[22] The DSD implementation of the catalysis reactions (54) and (59) are shown in Fig. 3. The DSD implementation of the catalysis reactions (55), (56), (60) and (61) are shown in Fig. 4. The reactions (57), (58), (62) and (63) are degradation reactions. The DSD implementation of reactions (57), (58), (62) and (63) are shown in Fig. 5. Reactions (64)–(66) are used as the synchronization module. The DSD implementation of reactions (64)–(66) are shown in Fig. 6.

Fig. 5

Degradation reaction module Y → ∅: DNA implementation of formal species Y and degradation with auxiliary species F; in (16), species Y reacts with species F producing inert waste, where the initial concentration of auxiliary species F is set to C.

Simulation results

The parameters of the 3-variable Lotka–Volterra oscillator are given in Table 1, and [X1]0, [X2]0, [X3]0, [Y1]0, [Y2]0 and [Y3]0 represent the initial values of species X1, X2, X3, Y1, Y2 and Y3, respectively. [X1]0 = 20 nM, [X2]0 = 12 nM, [X3]0 = 8 nM, [Y1]0 = 1 nM, [Y2]0 = 50 nM, and [Y3]0 = 15 nM, where the values of chemical reaction rates r, (i = 1,⋯,5) are referenced[1].

Parameter Values of the Lotka–Volterra Oscillator

	Values		Values		Values
r 1	5.5/300 s−1	q 1	r 1 C −1 m	C m	2 × 10−5 M
r 2	106 M s−1	q 2	105 M s−1	q m	106 M s−1
r 3	1.008 × 106 M s−1	q 3	r 2
r 4	4/300 s−1	q 4	105 M s−1
r 5	4/300 s−1	q 5	r 3
r 6	0.8/300 s−1	q 6	105 M s−1
		q 7	r 4 C −1 m
		q 8	r 5 C −1 m
		q 9	2 × 105 M s−1
		qs	r 6(qm + q9)/(q9Cm)

The evolution of the 3-variable Lotka–Volterra oscillator based on ideal CRNs and the 3-variable Lotka–Volterra oscillator based on DNA CRNs are shown in Fig. 7 and 8. When the three coupling terms are added to the system, the evolution of the 3-variable Lotka–Volterra oscillator based on DNA CRNs and the 3-variable Lotka–Volterra oscillator based on ideal CRNs are given in Fig. 9. According to Fig. 10, it can be found that the error between Y and X (i = 1, 2, 3) approaches zero, which indicates that synchronization among X and Y with different initial concentrations is realized. It is clear that errors between X and Y in the ideal CRNs approach zero, which is similar to errors in DNA CRNs, as shown in Fig. 10 and 11.

Fig. 7

Evolution of concentrations of x1, x2, x3, y1, y2 and y3 of two 3-variable Lotka–Volterra oscillator systems based on ideal CRNs and DNA intermediate CRNs, respectively.

Fig. 10

Process of synchronization of two 3-variable Lotka–Volterra oscillator systems based on DNA and ideal CRNs.

Fig. 11

Evolution of errors between two 3-variable Lotka–Volterra oscillator systems based on DNA and ideal CRNs.

In CRNs, luminous intensity relates to the concentration of DNA strands. We can observe the process of synchronization of the two 3-variable Lotka–Volterra oscillator systems through luminous intensity. The luminous intensity of X and Y are different without adding coupling terms to the systems. The luminous intensity of X and Y reach unanimity when there is synchronization among X and Y.

Conclusion

We proposed a 3-variable Lotka–Volterra oscillator system as an example and realized synchronization of two 3-variable Lotka–Volterra oscillator systems with different initial concentrations. For more convenient observation and more accurate detection of synchronization, DNA strands were designed with fluorophores and quenchers. We demonstrated the contributions of fluorophores and quenchers using the CRNs modules, and analyzed the differences between DNA CRNs and ideal formal CRNs. We used Visual DSD software as a design tool for the DNA strands and CRN language to simulate the synchronization of two 3-variable Lotka–Volterra oscillator systems, one each in DNA CRNs and ideal CRNs. Numerical simulation demonstrates that the method of synchronization in this article is effective for ideal CRNs as well as DSD implementation CRNs.

Conflicts of interest

There are no conflicts to declare.