Synthesis and materialization of a reaction-diffusion French flag pattern

TL;DR

This paper demonstrates the creation of a synthetic French flag pattern using DNA reaction networks, enabling autonomous pattern formation and control over colloidal organization in non-living matter.

Contribution

It introduces a novel DNA-based reaction-diffusion system that synthesizes and materializes a French flag pattern, reprogrammable via network mutation.

Findings

Successful synthesis of a sharp, long-lasting French flag pattern.

Control of colloidal aggregation patterns through concentration gradients.

Potential for testing reaction-diffusion models and fabricating autonomous soft materials.

Abstract

During embryo development, patterns of protein concentration appear in response to morphogen gradients. These patterns provide spatial and chemical information that directs the fate of the underlying cells. Here, we emulate this process within non-living matter and demonstrate the autonomous structuration of a synthetic material. Firstly, we use DNA-based reaction networks to synthesize a French flag, an archetypal pattern composed of three chemically-distinct zones with sharp borders whose synthetic analogue has remained elusive. A bistable network within a shallow concentration gradient creates an immobile, sharp and long-lasting concentration front through a reaction-diffusion mechanism. The combination of two bistable circuits generates a French flag pattern whose 'phenotype' can be reprogrammed by network mutation. Secondly, these concentration patterns control the macroscopic organization of DNA-decorated particles, inducing a French flag pattern of colloidal aggregation. This experimental framework could be used to test reaction-diffusion models and fabricate soft materials following an autonomous developmental program.

Tables

Table 1. Supplementary Table 1: Experimental conditions used for the experiments presented in the Main Text Figures.

Fig.	Network topology	Initial conditions
1		T${}_{A_1}{}^{}=25$ nM, $A_1=1$ nM, $R_1=0-400$ nM
2		$T_H=0-200$ nM, $T_K=20$ nM, $T_{H2R_K}=20$ nM, $T_{K2R_H}=20$ nM, $I_K=0$ nM, $I_H=0$ nM, $H=0.5$ nM, $K=0.5$ nM
3A		$T_{A_2}=75$ nM, $T_{A_3}=25$ nM, $A_1=10$ nM$†$, $A_3=10$ nM, $R_2-R_3=0-400$ nM, $T_{A_2}=25$ nM
3B		$R_2=50$ nM, $T_{A_3}=25$ nM, $A_1=10$ nM$†$, A3=10 nM, $T_{A_2}-R_3=0-75$ nM
4C-E		$T_{A_1}=25$ nM, $T_L=50$ nM, $A_1=10$ nM, $R_1=0-400$ nM, $L_1=0$ nM
4F		$R_1=0-400$ nM, $R_3=0-150$ nM, $A_1=A_3=10$ nM, $T_{A_1}=T_{A_3}=25$ nM, $T_{L1}=T_{L2}=50$ nM

Table 2. Supplementary Table 2: Sequences of oligonucleotides used in this work. Asterisks stand for phosphorothioate backbone modifications (3 or 4 phosphorothioates on the 5’ end protect the strand from being degraded by the exonuclease). Nicking enzyme recognition sites are in bold, Nb.BsmI recognizes the site 5’- GAATG ′ C N-3’, and Nt.BstNBI the site 5’-N ′ NNNN GACTC -3’, where N is any nucleotide and ′ indicates the position in the complementary strand where a nick is created. The presence of a U in the nicking binding site of Nt.BstNBI, strongly reduces its binding affinity ( 26 ). bt indicates a biotin, P indicates a phosphate. int.Spacer18 corresponds to a hexaethyleneglycol spacer. CY3.5, DY530, FAM and TexasRed are the used fluorophores.

Network	Name	Sequence $\left(5^{\prime }-3^{\prime }\right)$
	A1$†$	CATTCAGGATC
	A2$†$	CATTCAGGATCG
	T${}_{A_1}$	CY3.5-C*G*A*T*CCTGAATG′CGATCCTGAA
	R1	T*T*T*T*TCGATCCTGAATG-P
	R2$‡$	bt-AAAAAACGATCCTGAATG-P
1-species	T${}_{A_2}$$‡$	bt-*A*A*G*ATCCTGAATG′CGATCCTGAAT
bistable	R2-R3$‡$	bt-AAAAAACGATCCTGAATG-int.Spacer18- T*T*T*T*CTCGTCAGAATG-P
	T${}_{A_2}$-R3$‡$	bt-*A*A*G*ATCCTGAATG′CGATCCTGAAT-int.Spacer18-T*T*T*T*CTCGTCAGAATG-P
	A3	CATTCTGACGAG
	T${}_{A_3}$	DY530-*C*T*C*GTCAGAATG′CTCGTCAGAA
	R3	T*T*T*T*CTCGTCAGAATG-P
	H	CTGAGTCTTGG
	K	TCGAGTCTGTT
	TH$\S$	TexasRed-C*C*A*AGACUCAG′CCAAGACTCAGTTTTT-bt
2-species	TK	A*A*C*AGACUCGA′AACAGACTCGA-P
bistable	IH	GTCTTGGCTGAGTAA
	IK	GTCTGTTTCGAGTAA
	T${}_{H2R_K}$	T*T*A*CTCGAAACAGAC′CCAAGACTCAG-FAM
	T${}_{K2R_H}$	T*T*A*CTCAGCCAAGAC′AACAGACTCGA-DY530
	B${}_{}{}^r{}_1{}^{}$	G*G*A*TGAAGATGAGCATTACTTTCCGTCCCGAGAGACCTAACTGACACGCTTCCCATCGCTA-bt
	B${}_{}{}^l{}_1{}^{}$	bt-AGCATTACTTTCCGTCCCGAGAGACCTAACTGACACGCTTCCCATCGCTAGGATGAAGATG-P
	TL1	T*T*G*GATGAAGATGGGATGAAGATGGAATG’CGATCCTGAATG-P
bead	L1	CATCTTCATCCCATCTTCATCCAA
	S	T*A*G*CGATGGGAAGCGTGTCAGTTAGGTCTCTCGGGACGGAAAGTAATGC-P
sequences	B${}_{}{}^r{}_2{}^{}$	A*A*G*TAGGAGTAAGCATTACTTTCCGTCCCGAGAGACCTAACTGACACGCTTCCCATCGCTA-bt
	B${}_{}{}^l{}_2{}^{}$	bt-AGCATTACTTTCCGTCCCGAGAGACCTAACTGACACGCTTCCCATCGCTAGATGTGGAGAG-P
	TL2	T*T*G*ATGTGGAGAGAAGTAGGAGTAGAATG’CTCGTCAGAATG-P
	L2	TACTCCTACTTCTCTCCACATCAA

Table 3. Table 3: Numerical values of the parameters used to numerically solve Equations ( 6 – 7 ). cau means concentration arbitrary units. All the parameters were extracted from a manual fit to a series of independent well-mixed experiments except for D 𝐷 D that was measured in ref. [ 4 ] .

Parameter	Value
$k_p^H$	0.012 nM-1min-1
$k_p^K$	0.15 min-1
$k_p^{I_H}$	0.005 min-1
$k_p^{I_K}$	0.02 min-1
$k_{deg}$	0.05 min-1
$k_{deg}^I$	0.005 min-1
$\lambda$	0.01 cau-1
${\lambda }_I$	10 cau-1
$D$	0.018 mm2/min
$ϵ$	${10}^{-7}$ cau/min
$H_0$	${10}^{-3}$ cau
$K_0$	${10}^{-3}$ cau
Integration time	900 min
Integration distance	30 mm

Equations

$${\partial }_{t}M(x,t)=D{\partial }_{xx}M(x,t)$$

$${\partial }_{t}M(x,t)\le D{\partial }_{xx}M(x,0).$$

$${\partial }_{t}M(x,t)\le \frac{D}{l^2}M(x,0).$$

$$\frac{\Delta M(x)}{M(x,0)}\le \frac{D\Delta t}{l^2},$$

$$\frac{\Delta M(x)}{M(x,0)}\le \mathrm{0.1.}$$

$$\frac{\partial H}{\partial t}$$

$$\frac{\partial K}{\partial t}$$

$$\frac{\partial I_H}{\partial t}$$

$$\frac{\partial I_K}{\partial t}$$

$$H(x,0)=H_0,K(x,0)=K_0,I_K(x,0)=I_K(x,0)=0,$$

$$\frac{\partial A(x,t)}{\partial t}=A(x,t)\left(\frac{k_p}{K_p+A(x,t)}-\frac{k_{d}R(x)}{K_{d}R(x)+A(x,t)}-\frac{k_r}{K_r+A(x,t)}\right)+D\frac{{\partial }^{2}C(x,t)}{\partial x^2}+ϵ$$

Full text

Synthesis and materialization of a reaction-diffusion French flag pattern

Anton S. Zadorin1,2, Yannick Rondelez3,4, Guillaume Gines3, Vadim Dilhas1,2,

Georg Urtel5, Adrian Zambrano1,2, Jean-Christophe Galas1,2,∗, André Estevez-Torres1,2,∗

1Université Pierre et Marie Curie, Laboratoire Jean Perrin, 4 place Jussieu, 75005 Paris, France.

2CNRS, UMR 8237, 75005, Paris, France. 3LIMMS/CNRS-IIS, University of Tokyo,

Komaba 4-6-2 Meguro-ku, Tokyo, Japan. 4Ecole supérieure de physique et chimie industrielle,

Laboratoire Gulliver, 10 rue Vauquelin, 75005, Paris, France. 5 Ludwig-Maximilians-Universität München,

Fakultät für Physik, Amalienstraße 54, 80799 Munich, Germany

∗To whom correspondence should be addressed;

E-mail: jean-christophe.galas@upmc.fr, andre.estevez-torres@upmc.fr

**During embryo development, patterns of protein concentration appear in response to morphogen gradients. These patterns provide spatial and chemical information that directs the fate of the underlying cells. Here, we emulate this process within non-living matter and demonstrate the autonomous structuration of a synthetic material. Firstly, we use DNA-based reaction networks to synthesize a French flag, an archetypal pattern composed of three chemically-distinct zones with sharp borders whose synthetic analogue has remained elusive. A bistable network within a shallow concentration gradient creates an immobile, sharp and long-lasting concentration front through a reaction-diffusion mechanism. The combination of two bistable circuits generates a French flag pattern whose ’phenotype’ can be reprogrammed by network mutation. Secondly, these concentration patterns control the macroscopic organization of DNA-decorated particles, inducing a French flag pattern of colloidal aggregation. This experimental framework could be used to test reaction-diffusion models and fabricate soft materials following an autonomous developmental program. **

From a chemist’s perspective, biological matter has the astonishing capability of self-constructing into shapes that are predetermined, robust to varying environmental conditions and remarkably precise in size and chemical composition at different scales. Living embryos, for instance, develop from a simple form into a complex one through a reproducible process called embryogenesis. The embryo is first structured chemically through pattern formation, a process during which out-of-equilibrium molecular programs generate highly ordered concentration patterns of $\mu$m to mm size[1]. Subsequent developmental steps involve morphogenesis, where the embryo changes its shape, cell differentiation, in which cells become structurally and functionally different, and finally growth, resulting in an increase of the mass of the embryo[1]. The emulation of such processes in non-living chemical systems adresses two important goals. First, in a reductionist perspective, it allows testing theoretical models describing the emergence of out-of-equilibrium material order in simplified experimental conditions[2]. Second, from a synthetic standpoint, it enables the conception of a new way of making soft materials[3, 4] that one may call ’artificial development’: materials that build themselves following a pre-encoded molecular program.

The first developmental step, pattern formation, is currently interpreted in the light of two archetypal scenarios: Turing’s instability [5] and Wolpert’s processing of positional information [6]. The Turing scenario implies an initially homogeneous reaction-diffusion (RD) system that spontaneously breaks the symmetry, generating repetitive structures of wavelength $\sqrt{D\tau}$, where $D$ is a diffusion coefficient and $\tau$ a characteristic reaction time. In Wolpert’s scenario, instead, the symmetry is already broken by a pre-existing morphogen gradient that is subsequently interpreted in a threshold-dependent manner, producing several chemically-distinct regions with sharp borders. Lewis Wolpert named this the French flag problem to illustrate the issue of creating three distinct regions of space —the head, the thorax and the abdomen— from an amorphous mass and a shallow concentration gradient (Fig. 1a). Historically, Turing’s and Wolpert’s scenarios have been considered mutually exclusive[7], the first needing diffusion in contrast with the second. However, recent hypothesis suggest that both scenarios could be combined during development[7] and that diffusion could make the processing of positional information more robust[8, 9]. The experimental proof of the Turing mechanism in vivo has met constant criticism[7], although recent work has provided new evidence[10, 11]. In contrast, Wolpert’s scenario is accepted in living embryos, possibly because it is more loosely defined. Well-known examples are provided by the gap gene system in Drosophila[12] and by the sonic hedgehog-induced patterning of the vertebrate neural tube[13]. Concerning non-living systems, Turing patterns were first demonstrated experimentally in 1990[14] whereas synthetic systems capable of interpreting a morphogen pre-pattern have remained elusive.

Here, we draw inspiration from pattern formation during early Drosophila development to address a synthetic challenge: can one build a French flag pattern outside of a living organism? In a learning-by-doing approach[15, 16, 17, 18, 19, 20] to this question our synthetic route combines reaction-diffusion with positional information, Turing and Wolpert ideas, suggesting that these two mechanisms are not antinomic. Furthermore, by considering pattern formation in the context of development as a blueprint for cell differentiation we ask a second question: may a self-organized chemical pattern influence the final structure of an initially homogeneous material? In doing this, we demonstrate that a soft material can be autonomously structured through an artificial developmental program.

Results

The conversion of a shallow morphogen gradient into a concentration boundary that is both sharp and immobile requires a reaction network that interprets the gradient in a non-linear fashion. Diverse evidence [21, 22, 23, 8] suggests that bistability is an essential property of such networks and that coupling bistability with diffusion provides immobile and robust fronts[23, 8, 9]. Our design for building a French flag pattern thus consists of two bistable networks that generate two RD fronts of concentration pinned in a gradient of a bifurcation parameter.

DNA oligonucleotides are particularly well-suited to construct pattern-forming reaction networks.[18, 17, 24, 4] The reactivity of the hybridization reaction obeys simple rules and a wide array of methods coming from biotechnology renders their synthesis, analysis and modification straightforward. We thus engineered a series of bistable networks using the PEN DNA toolbox, a molecular programming language designed to construct networks analogous to transcriptional ones but using only simple biochemical reactions.[26] This technology has recently been applied to construct out-of-equilibirum networks displaying oscillations[26, 27], bistability[7] and traveling concentration waves.[17, 4, 29] Fig. 1b depicts the simplest bistable network used here, with a first-order positive feedback loop and a non-linear repressor[8] (Supplementary Figs. 3-4). The nodes of the network, A1 and R1, are respectively 11 and 15-mer single-stranded DNAs (ssDNAs). Self-activation is set by T${}_{A_{1}}$, a 22-mer ssDNA template. Repression is encoded by promoting the degradation of A1 with a threshold given by $R_{1}$[8] —italized species names indicate concentration throughout the text. Three enzymes —a polymerase, an exonuclease and a nicking enzyme— catalyze the three basic reactions —DNA polymerization, ssDNA degradation and the nicking of double stranded DNA (dsDNA) in the presence of the correct recognition sequence— and dissipate free energy from a reservoir of deoxynucleoside triphosphates (dNTPs). In comparison with gene regulatory networks, template T${}_{A_{1}}$ plays the role of a gene, encoding information, A1 and R1 are analogous to transcription factors, as they promote or repress the activity of the template, and the three enzymes provide the metabolic functions homologous to the transcription-translation machinery. Moreover, A1 is continuously produced and degraded but the total template —gene— concentration is fixed. In contrast to networks in vivo, molecular interactions are well-known and the mechanism and kinetic rates can be precisely determined[26].

Our first goal was to create an immobile concentration front in the presence of a morphogen gradient, which we have called a Polish flag pattern. We performed patterning experiments within 5 cm-long sealed glass microchannels of $4\times 0.2$ mm2 cross-section (Fig. 1c). A gradient of morphogen R1 was generated along the longitudinal axis of the channel, noted $x$, by partially mixing two solutions with different $R_{1}$ by Taylor dispersion (Fig. 1e and Supplementary Fig. 5). The gradient was well approximated by an exponential decay $e^{(x-x_{0})/l}$ with characteristic length $l=2$ cm. Because diffusion is slow over long distances the gradient was stable over 50 h (within 10% at the center of the channel, as expected for the diffusion of a 17-mer ssDNA, see Supplementary Fig. 5). Initially, the channel contained homogeneous concentrations of the three enzymes, dNTPs, $T_{A_{1}}=25$ nM, $A_{1}=1$ nM, and a gradient of $R_{1}$ in the range $0-200$ nM. Throughout this work, the concentration of the network nodes, here A1, was related to fluorescence intensity using two methods (see Methods). Either by adding the DNA intercalator EvaGreen, for which the fluorescence signal is proportional to the concentration of dsDNA, or by labeling one template with a fluorophore that is quenched upon hybridization. In order to facilitate the interpretation of the data we represent in all figures the fluorescence shift, which is the absolute value of the difference between the fluorescence intensity at a given time and at initial time. As a consequence, at low concentration, the fluorescence shift is proportional to the concentration of the node species (see Methods). The fluorescence inside the channel was measured by recording time-lapse images with an inverted microscope. Fig. 1d displays the spatio-temporal dynamics of the patterning process. A short, purely reactional initial phase generated a sharp profile of $A_{1}$ at a location corresponding to low morphogen concentration (Supplementary Fig. 6). This profile later moved to the right through a RD mechanism, progressively decelerating until it stopped at the center of the channel at a position where $R_{1}^{RD}=30\pm 5$ nM. The front remained immobile up to 15 h and its characteristic width, defined as the decay length of a sigmoidal fit to the data, was $\lambda=2$ mm, 10-fold sharper than the morphogen gradient (Fig. 1d and Supplementary Fig. 6). When, instead of the repressor, the autocatalyst template was used as the morphogen, the complementary Polish flag pattern was obtained (Supplementary Figs. 8-9).

The gap gene network, which interprets the Bicoid morphogen gradient during the development of the Drosophila blastoderm, is not composed of a unique self-activating and repressed node but of a series of them[6] (Supplementary Fig. 10). To demonstrate that our approach is capable of emulating the most basic type of such networks, we used one with two self-activating nodes, H and K, that repress each other[7] (Fig. 2a and Supplementary Fig. 11) and recorded the patterning dynamics in a gradient of the Bicoid analogue, TH (Fig. 2b-d, Movie S1), the template corresponding to autocatalyst H. At short times, a purely reactional phase created two independent and sharp fronts of H and K. Subsequently, during an RD phase, the fronts traveled in opposite directions until they collided in the middle of the channel. At this time, the two profiles partially overlapped and a slow phase made the two fronts go backwards, repelling each other, until reaching a steady-state where the overlap disappeared. 1-dimensional simulations with a 4-variable model (Supplementary Section 3.1 and Supplementary Fig. 12) displayed a similar behavior and suggested that this last phase was due to a slow synthesis of the repressors (Fig. 2e-g). This patterning process was highly reproducible (Supplementary Figs. 13-14) and compatible with the immobilization of the morphogen gradient onto a surface (Supplementary Fig. 15).

To implement a French flag pattern with three chemically-distinct zones (Fig. 3) we combined two orthogonal bistables, A2 and A3, (Supplementary Fig. 16) into a single network using two different approaches. We used a bifunctional morphogen bearing either both repressors, $\mathrm{R}_{2}-\mathrm{R}_{3}$, or one autocatalyst template and a repressor, $\mathrm{T}_{A_{2}}-\mathrm{R}_{3}$. In a gradient of $\mathrm{R}_{2}-\mathrm{R}_{3}$, a channel containing a uniform concentration of $\mathrm{T}_{A_{2}}$ and $\mathrm{T}_{A_{3}}$ generated a French flag pattern that divided space into three regions with different composition, $\textrm{A}_{2}+\textrm{A}_{3}$, $\textrm{A}_{2}$ and $\emptyset$, for 100 min (Fig. 3a, Supplementary Figs. 17-18 and Supplementary Movie 2). By contrast, with a gradient of $\mathrm{T}_{A_{2}}-\mathrm{R}_{3}$ and a uniform concentration of $\mathrm{R}_{2}$ and $\mathrm{T}_{A_{3}}$ a different pattern separated the space into $\textrm{A}_{3}$, $\textrm{A}_{2}+\textrm{A}_{3}$ and $\textrm{A}_{2}$ (Fig. 3b).

In embryogenesis, pattern formation can induce tissue differentiation by providing localized chemical cues to pluripotent cells [1]. This strategy could be used for the autonomous fabrication of artificial materials. As a proof of concept we coupled the patterns obtained above to the conditional aggregation of 1 $\mu$m diameter beads[32] (Fig. 4). Streptavidin-labeled beads were decorated with two types of biotin-labeled DNAs that had two different 12-mer ssDNA dangling ends, B${}_{i}^{l}$ and B${}_{i}^{r}$; for the pair of beads $i$, one has a left and the other a right strand. In the working buffer, the beads aggregated only in the presence of a linker strand Li complementary to both left and right ssDNA portions (Supplementary Fig. 19). A capillary containing i) a homogeneous dispersion of beads B${}_{i}^{l}$ and B${}_{i}^{r}$, ii) a bistable network with node Aj coupled to the linear production of Li and iii) a gradient of $\textrm{R}_{j}$, produced a Polish flag of bead aggregation (Fig. 4c-e for $(i=1,j=1)$ and Supplementary Fig. 20 for $(i=2,j=3)$).

It is straightforward to generate two pairs of beads, (B${}_{1}^{l}$, B${}_{1}^{r}$) and (B${}_{2}^{l}$, B${}_{2}^{r}$), each pair aggregating independently of the other in the presence of its own linker (L1 or L2, Supplementary Fig. 21-23). The French flag pattern in Fig. 3a was thus materialized into three zones of space with distinct degree of aggregation (Fig. 4f and Supplementary Fig. 24): both pairs of beads aggregated / one pair of beads aggregated / no bead aggregated. In addition, bead aggregation brought two interesting properties. It allowed the visualization of RD patterns by eye (Supplementary Figs. 20 and 26) without the need of a fluorescence microscope and it froze the patterns into a state that was stable for at least one month, because aggregates precipitated to the bottom of the capillary (Supplementary Fig. 26). In other words: a steady-state dissipative structure of DNA was converted into a kinetically-trapped stable structure of beads.

Discussion

The object of synthesis in far from equilibrium chemistry is not anymore a molecular structure with desirable physico-chemical properties but a network of chemical reactions with particular dynamics. When these networks are coupled with some kind of transport —diffusion is just one possibility[33]— a length scale naturally emerges, which allows to structure space chemically. The first implication of our work is to provide an experimental framework for the synthesis of far from equilibrium chemical structures. Indeed, among the limited amount of frameworks that have been proposed to synthesize RD patterns[34, 35], few of them are both biocompatible and programmable. By programmable we mean that one may rationally choose which reactants will react with what mechanism to create a desired pattern. The PEN DNA toolbox used here is naturally biocompatible and DNA sequence complementarity makes it programmable —both for the topology of the network[26, 7, 17], as shown in Fig. 1-4 and Supplementary Fig. 2, and for the reaction and diffusion rates[4]. Biocompatibility opens the way for the structuration of biomolecules and living cells with RD patterns. Programmability makes it suitable to synthesize new patterns, such as two-dimensional RD solitons[36, 37], to probe experimentally the degree of complexity that may emerge from RD patterning. The versatility of this method allows conceiving a new way of processing materials inspired from embryo development: embed them with an autonomous developmental program and let them construct themselves.

We may indeed consider embryogenesis as an extremely precise and autonomous procedure to structure matter. In contrast with current fabrication methods, it is able to position chemicals with multiscale precision —from tens of nm for cellular organelles to 10 $\mu$m for cells—, outstanding reproducibility and controlled dynamics without mechanical parts. In this regard, RD mechanisms have already been used to process materials, notably for micro- and nanofabrication using precipitation reactions[38, 39]. Our approach significantly enlarges the complexity of the underlying reaction networks and, because information is encoded in DNA sequence, it could be coupled with directed evolution. If an autonomous fabrication method of this sort is once to find a real application it will need to build materials that are chemically diverse with a resolution at least comparable with the one in Drosophila patterning, which is 10 $\mu$m. Recent work demonstrates that it is possible to couple DNA species with interesting chemistry such as DNA nanostructures[40], aptamers[41], hydrogels[42], chemical synthesis[43] or gene delivery[44]. Although the current $2$ mm resolution of our patterning chemistry is low, it is far from its theoretical limit. Indeed, RD patterning has been used for fabricating 300 nm objects[38]. The spatial resolution is determined by $\lambda\sim\sqrt{D/k}$, where $k$ is the degradation rate (Supplementary Fig. 27). Here $\lambda=2$ mm and taking $D=1.8\times 10^{4}\leavevmode\nobreak\ \mu$m2/min for a 12-mer ssDNA[4] at $45^{\circ}$C one finds $k\approx 5\times 10^{-3}$ min*-1*. This value is in good agreement with the measured degradation kinetics of species W2 in the presence of R2, with rates in the $3\times 10^{-2}-3\times 10^{-3}$ min*-1* range (Supplementary Fig. 28). Improving the resolution down to 10 $\mu$m thus requires decreasing $D$ and increasing $k$ each by a factor 100, which is not unreasonable (Supplementary Fig. 28). In addition, coupling RD with other sharpening mechanisms could increase resolution further: the bead aggregation front was 4-fold sharper than the RD front, which may be due to the cooperativity of bead aggregation[45]. This multi-level sharpening recalls hierarchical patterning mechanisms in vivo. Incidentally, it appears plausible to make 100 $\mu$m-scale gradients by either using surface microprinting (Supplementary Fig. 15) or a self-generating diffusion-degradation mechanism[46, 3].

The second implication of our work is that it may help understanding the physics of chemical patterning in living systems[48, 49]. Although our experimental system is an extreme simplification —no cells, no molecular crowding, no forces—, it integrates non-trivial microscopic interactions that are present in vivo —intra and inter-molecular forces, DNA hybridization and protein-DNA recognition— and all the time and spatial scales of reaction and diffusion are naturally included with actual molecules. It could thus provide an interesting framework to perform experimental simulations instead of computer ones when these have limitations, for instance when noise becomes important.

Conclusion

Our results demonstrate that DNA-based molecular programming is well-suited to engineer concentration patterns that emulate those observed during early embryo development. They indicate that the combination of a bistable reaction network with diffusion provides a simple engineer’s solution to generate immobile concentration fronts that are both sharp and long lasting. Importantly, the simplicity of the method allowed us to record the patterning dynamics in real time, showing that a purely reactional initial phase is followed by a reaction-diffusion one. Our experimental model may help understanding the role of regulative and diffusive processes during development and suggests that relatively simple networks may have enabled patterning at an early stage of evolution. Finally, by coupling programmable patterns with matter we have demonstrated a primitive autonomous chemical structuration of a material in one dimension. This approach could be exploited to fabricate soft materials following an autonomous developmental program.

Methods

DNA strands were purchased from Biomers (Ulm, Germany). The Bst DNA polymerase large fragment and the two nicking enzymes (Nb.BsmI and Nt.BstNBI) were purchased from New England Biolabs. The recombinant Thermus thermophilus RecJ exonuclease was produced as described[1]. Experiments were performed at 45 °C for the Polish and French flag generating networks and at 42 °C for those in Fig. 2. Details on the the sample preparation, the DNA sequences, the experimental conditions, the data analysis and the simulations are provided in the Supplementary Information.

Bead suspension.

1 $\mu$m diameter, streptavidin-coated, paramagnetic beads (Dynabeads MyOne C1, Invitrogen) were functionalized with two types of biotinylated DNA constructs, as described[2], making a pair of beads, (B${}_{i}^{l}$, B${}_{i}^{r}$). Each construct consisted of a 49 bp-long dsDNA backbone terminated with a 12 bases-long single stranded sticky end. The construct corresponding to B${}_{i}^{l}$ (resp. B${}_{i}^{r}$) was biotin-labeled on the 5’ end (resp. 3’ end) and the corresponding sticky end was on the 3’ side (resp. 5’ side).

Measurement of DNA concentrations.

DNA concentrations were measured by fluorescence. To measure the concentration of autocatalytic nodes Ai, H and K we used two different strategies. i) All experiments involving Ai contained EvaGreen (Biotium), an intercalating dye which fluorescence is proportional to the concentration of dsDNA. ii) Some template strands were labeled with fluorescent dyes on their 3’ or 5’ ends (see Supplementary Table 2 for details) which fluorescence was quenched specifically by the corresponding complementary strands. For experiments with networks involving a single node (Fig. 1) we used EvaGreen. For networks involving two nodes we either used two orthogonal N-quenching dyes (Fig. 2) or EvaGreen and an orthogonal N-quenching dye (Fig. 3). The concentration of morphogen was measured either by using a fluorescently-labeled strand (Fig. 2) or by adding a cascade blue-dextran M${}_{w}=3000$ Da (Thermo Fisher Scientific) (Fig. 1 and 3-4), and recording its fluorescence in real time. This fluorescent dextran has a diffusion coefficient similar to the DNA templates (Supplementary Fig. 5). All kymographs and plots display the fluorescence shift, i.e the absolute value of the difference between the fluorescence intensity at a given time and at initial time. The fluorescence shift is proportional to the concentration of the species of interest if this remains below the concentration of its complementary strand.

Generation of the morphogen gradient.

Spatiotemporal experiments were performed within 50 mm $\times$ 4 mm $\times$ 0.2 mm glass capillaries (Vitrocom USA). One half of the capillary was loaded with the reaction solution without the morphogen strand and the other half with the same solution supplemented with the morphogen at 200 nM. To create the gradient the two solutions were mixed by applying a hydrodynamic flow along the capillary axis using a micropipette and a custom-made PDMS connector. 15 up-and-down pumps of 12.5 $\mu$L yielded a shallow morphogen gradient along the entire capillary length spanning between 0 and $\sim$100 nM. Subsequently, the capillaries were sealed with 5-minutes Araldite epoxy and glued over a $5\times 7.5$ cm glass slide.

Microscopy.

The fluorescence along the capillary was recorded on a Zeiss Axio Observer Z1 automated epifluorescence microscope equipped with a Tokai Hit thermo plate, an Andor iXon Ultra 897 EMCCD camera and a $2.5\times$ objective and controlled with MicroManager 1.4. For each capillary, 16 contiguous 3.17$\times$3.17 mm2 images were recorded automatically every 1 to 10 minutes. Multi-color fluorescence microscopy was used to record the concentration of different DNA species over time. Images of the beads were acquired in bright field with a 10 or a $40\times$ objective.

Data treatment.

The raw data were treated with ImageJ / Fiji (NIH) and Matlab (The Mathworks). The 16 images making one capillary were stitched together and corrected from inhomogeneous illumination. To obtain the kymographs, the images were averaged along the width of the capillary ($y$ axis) and the corresponding profiles stacked over time. These profiles were further smoothed by performing a moving average along $x$, normalized and fitted to a sigmoid function $f(x)=\frac{1}{1+e^{(x-x_{0})/\lambda}}$. The front position corresponds to $x_{0}$ and its width to $\lambda$. To determine the front velocity, the data of front position over time were fitted by a polynomial function to reduce noise. The velocity was calculated as the time derivative of this polynomial fit. To plot the bead aggregation profile along the capillary, the corresponding brightfield images were first binarized and the number of particles counted in areas 100 $\mu$m wide. A particle here indicates either an aggregate or a bead: the number of aggregates was thus inversely proportional to the number of detected particles.

Data availability.

All data generated or analysed during this study are included in this published article (and its supplementary information files). The raw datasets and the plasmid for expressing ttRecJ are available from the corresponding authors on reasonable request.

Acknowledgments

We thank E. Frey for insightful discussions, A. Vlandas for help with gradient generation and B. Caller and D. Woods for comments on the text. Supported by European commission FET-STREP (Ribonets), by ANR jeunes chercheurs program (Dynano), by C’nano Ile-de-France (DNA2PROT) and by Ville de Paris Emergences program (Morphoart). Correspondence and requests should be addressed to A.E.-T. (andre.estevez-torres@upmc.fr) or J.-C.G (jean-christophe.galas@upmc.fr).

Author contributions

A.S.Z., J.-C.G. and A.E.-T. performed most experiments and analyzed the data. Y.R. and G.G. designed the network in Fig. 1 and J.-C.G. and A.E.-T. designed the networks in Figs. 3-4. A.Z. and V.D. set up the bead experiments. G.U. performed critical control experiments. All the authors discussed the results. J.-C.G., A.S.Z., Y.R. and A.E.-T. designed research and J.-C.G. and A.E.-T. wrote the manuscript.

Competing financial interests

The authors declare no competing financial interests.

Table of contents summary

During embryogenesis patterns of protein concentration appear in response to morphogen gradients, providing spatial and chemical information that directs the fate of the underlying cells. Here, this process is emulated with DNA-based non-living matter and the autonomous structuration of a synthetic material is demonstrated.

Figure captions

Supplementary materials

1 Materials and methods

1.1 Preparation of solutions

Two types of bistable networks were used throughout the Main Text. A single species production with degradation network (Figures 1, 3 and 4) and a two species mutual inhibition network (Figure 2). They will be called, respectively, 1-species and 2-species bistable networks. The 1-species and 2-species bistable networks depend on two different nicking enzymes and thus the associated buffers slightly differ.

1.1.1 Components common to both buffers

1$\times$ thermopol buffer (NEB, containing 20 mM Tris-HCl, 10 mM (NH4)2SO4, 10 mM KCl, 2 mM MgSO4, 0.1% Triton X-100, pH 8.8@25°C), 1 g/L synperonic F108 (Sigma Aldrich), 3 mM Dithiothreitol (Sigma Aldrich), 50 mg/L Bovine Serum Albumin (NEB), 2$\mu$M Netropsin (Sigma Aldrich), 50 $\mu$M NaCl.

1.1.2 1-species bistable buffer

In addition to the components common to both buffers, the 1-species bistable buffer contains:

• •

MgSO4 6 mM,

• •

deoxynucleotidetriphosphates (0.8 mM of each dNTP) (NEB),

• •

Bst DNA polymerase large fragment (pol) (NEB) 0.1% of the 8,000 units/mL stock solution,

• •

Thermus thermophilus RecJ exonuclease (ttRecJ) expressed in house as described in [1], 12.5 nM (1% of the stock solution),

• •

Nb.BsmI nicking enzyme (nick) (NEB), 5% of the 10,000 units/mL stock solution.

• •

If templates bearing a biotin in the 5’ end were used, the buffer was supplemented with streptavidin at 200 nM in binding sites.

All experiments involving the 1-species bistable network were performed at 45°C.

1.1.3 2-species bistable buffer

In addition to the components common to both buffers, the 2-species bistable buffer contains:

• •

Tris-HCl 25 mM, pH 8,

• •

MgSO4 5 mM

• •

dNTPs (0.2 mM each) (NEB),

• •

Bst DNA polymerase large fragment (pol) (NEB) 0.2% of the 8,000 units/mL stock solution,

• •

ttRecJ 16 nM (1.33% of the stock solution),

• •

Nt.BstNBI nicking enzyme (nick) (NEB), 1% of the 10,000 units/mL stock solution.

All experiments involving the 2-species bistable network were performed at 42°C.

Note that the activity of both nicking enzymes occasionally changed from batch to batch, thus their concentrations were adjusted according to independent assays.

1.1.4 Particle suspension

The DNA-functionalized colloids used in this work are a slight variation of the design described in [2]. We used 1 $\mu$m diameter, streptavidin-coated (5 $\mu$M biotin binding sites and 10 mg/mL of beads) paramagnetic beads (Dynabeads MyOne C1, Invitrogen). The beads were decorated with two different biotinylated DNA constructs making two different beads, called B${}_{i}^{l}$ and B${}_{i}r$ (Supplementary Figure 1). Each construct consists of a 49 bp-long dsDNA backbone (S-S*) terminated with a 12 bases-long single stranded sticky end. The construct corresponding to B${}_{i}^{l}$ (resp. B${}_{i}^{r}$) was biotin-labeled on the 5’ end (resp. 3’ end) and the corresponding sticky end was on the 3’ side (resp. 5’ side). Such a design implies that the DNA-decorated beads will be stable in solution when mixed, unless the linker strand complementary to the sticky ends of each bead type is present. The subscript $i$ indicates that the pair of beads aggregate with linker strand Li. Li has 2 extra A bases on its 3’ end to impede polymerase extension and the subsequent strand displacement of S on B${}_{i}^{l}$ beads. The grafting protocol of DNA to get the two types of beads presented in Supplementary Figure 1 requires the following steps:

Two solutions with the two DNA constructs (B${}_{i}^{l}$ + S) and (B${}_{i}^{r}$ + S) are prepared at 8 $\mu$M final concentration of strands in the suspension buffer (10 mM phosphate, 50 mM NaCl and 0.1% w/w Pluronic surfactant F127, Sigma-Aldrich) and kept for at least 15 min at room temperature before use. 2. 2.

During this time, the beads are rinsed 3 times in the suspension buffer and split into two aliquots. 3. 3.

The bead supernatant is removed and the corresponding DNA construct solution is added to each bead aliquot and incubated at room temperature for 30 min under gentle mixing. The final concentration of the beads is 10 mg/mL. 4. 4.

The two beads solutions are mixed and rinsed 3 times in the suspension buffer, then maintained at 55°C for 30 min and rinsed again 3 times, to eliminate the non-grafted strands. The bead stock concentration is kept at 10 mg/mL. The solution is stored at 4°C.

Note that before each experiment, the bead stock solution is heated again at 55°C for 30 min. The bead solution is used at a concentration ranging from 0.1 to 0.5 mg/mL in the final reaction mix.

1.1.5 DNA concentrations for the experiments presented in the Main Text

$\dagger$ The node of the network corresponds well to species A2. However it was initiated by species A1 because the input side of templates T${}_{A_{1}}$ and T${}_{A_{2}}$ are closely related and both can be triggered by species A1.

1.2 DNA sequences and topology of the networks

The DNA oligonucleotides were purchased from Biomers with HPLC purification. All the reaction networks are presented in Supplementary Figure 2 and the sequences listed in Supplementary Table 2.

$\dagger$ Species A2 corresponds to the sequence of A1 where one base on the 3’ end has been removed. As a result, A1 can initiate the autocatalysts templated by T${}_{A_{1}}$ and T${}_{A_{2}}$.

$\ddagger$ bt in 5’ + streptavidin also protects the corresponding strand from exonuclease degradation. In this case the last two bases in 5’ (AA) are not replicated by the polymerase (see ref. [3] for details). The presence of both bt and phosphorothioate in 5’ is due to historical reasons.

$\mathsection$ The presence of a TTTTT-bt on the 3’ side allowed to attach this species to the surface through biotin-streptavidin (see Section 3.3). This modification was not needed for experiments in solution such as those presented in Figure 2 in the MT. However, for convenience, and because the modification does not play a role in the reaction network, the same sequence was used.

1.3 Measurement of DNA concentrations

DNA concentration over space and time was measured by fluorescence. We used three different strategies:

• •

Recording the green fluorescence from EvaGreen dye (1x EvaGreen DNA binder, 20x dilution of the manufacturer’s stock solution, Biotium). In this case, fluorescence is proportional to the concentration of double stranded DNA.

• •

Recording the fluorescence from dyes attached to the 3’ or 5’ end of template strands (see Supplementary Table 2 for details). When the corresponding inputs or outputs hybridize on these templates, fluorescence is quenched.

• •

The concentrations of morphogens were measured with two different methods. For the 2-species bistable (Figure 2 of the MT) template TH was labeled with a Texas-Red fluorophore. For the 1-species bistable we added 1 $\mu$M of cascade blue-dextran M${}_{w}=3000$ Da (Thermo Fisher Scientific, Molecular Probes) to the solution with high concentration of morphogen used to generate the gradient, and measured its fluorescence. The cascade blue-dextran has a molecular weight similar to the DNA templates and thus one expects a similar diffusion coefficient (Supplementary Figure 5e-f). In a control experiment we compared the fluorescence profile obtained with cascade blue-dextran with a fluorescent DNA template, they were similar (Supplementary Figure 5d). Note that in some experiments the gradient was imaged with ROX, a fluorophore that diffuses significantly faster than the morphogen. In those cases only the gradient at initial time is provided. In all other cases the gradient profile was recorded in real time.

The use of cascade-blue-dextran instead of a fluorescence-labeled oligonucleotide in 1-species experiments was due to constraints in the number of flurescence channels and the choice of fluorophores capable of N-quenching with the used sequences in the French flag experiments. The green channel was reserved for EvaGreen fluorescence. The yellow and red channels were used for DY530 and Cy3.5 dyes for monitoring species A1 and A3. Only the blue channel was left. Unfortunately, the signal of blue-fluorescent dyes attached to oligonucleotides was too low. We thus used cascade-blue-dextran at a much higher concentration than the morphogen to have enough signal.

In all kymographs and plots we represented the absolute value of the fluorescence shift, which is proportional to the concentration of fluorescent species in a concentration range between the detection limit and the concentration of the DNA strand complementary to the species of interest.

1.4 Temporal experiments

Temporal experiments were performed on a BioRad CFX or Qiagen Rotor-Gene Q PCR machine (both used as temperature-controlled fluorimeters) using 20 $\mu$L of solution in 150 $\mu$L PCR tubes. Fluorophores ROX or Cy5-A15 at 15 nM, where A15 is an oligonucleotide with 15 adenines, were used as internal standards. The raw fluorescence intensity was divided by the intensity of the reference internal standard at each time point.

1.5 Spatio-temporal experiments

Spatio-temporal experiments were performed within 50 mm$\times$4 mm$\times$0.2 mm glass capillaries (Vitrocom, USA). The capillary was loaded using a micropipette and a custom-made PDMS connector.

Generation of the morphogen gradient.

The capillary was filled with 45 $\mu$L of the reaction solution without morphogen strand and then a pipette was inserted being in the ’push’ position into the PDMS connector. The other end of the capillary was dipped into the reaction solution containing the morphogen DNA strand. 15 up-and-down pumps of 12.5 $\mu$L of solution were performed to create the gradient. Finally, the pipette was left in the ’pulled’ position and the pipette together with the PDMS connector were removed.

Microscopy.

Once the gradient was formed, the capillary was laid on a $5\times 7.5$ cm glass slide. 5-minutes Araldite epoxy was used, both to seal the capillary ends and to glue them to the glass slide. No evaporation at all was observed for 48 h at 45°C. The fluorescence along the capillary was recorded on a Zeiss Axio Observer Z1 fully automated epifluorescence microscope equipped with a CoolLED pE-2 fiber-coupled illuminator, DAPI, GFP, YFP and RFP filter sets, a Marzhauser XY motorized stage, a Tokai Hit thermo plate, and Andor iXon Ultra 897 EMCCD camera, a shutter and a $2.5\times$ objective. These instruments were controlled with MicroManager 1.4. For optimal thermal conduction, mineral oil was added between the glass slide and the thermoplate on the microscope, and also between the glass slide and the capillary. For each capillary, 16 contiguous 3.17$\times$3.17 mm2 ($128\times 128$ pixel2) images were recorded automatically every 1 to 10 minutes. Multi-color fluorescence microscopy was used to record the concentration of different DNA species over time. Images of the beads were acquired in bright field with a 10 or a $40\times$ objective.

Image treatment.

The raw data were treated with ImageJ / Fiji (NIH) and Matlab (The Mathworks). Prior to data analysis, the 16 images were stitched together without overlapping. Subsequently, the inhomogeneous illumination was corrected by two different protocols. When the initial concentration was flat (for all species except for the morphogen) a division by the first frame was performed. When the initial concentration was not homogeneous, for example for the morphogen, a polynome was fitted to the raw data of one of the 3.17$\times$3.17 mm2 of the first frame and each image was divided by it.

The counting of the particles (aggregates or single bead without distinction) has been done on images acquired in bright field with a 10 or a $40\times$ objective. Contiguous images were stitched together without overlapping. Images were first binarized. Then, the Fiji Particle Analysis plugins was used to determine the number of particles in contiguous areas (100$\mu$m to 250$\mu$m wide) along the capillary.

1.6 Data treatment

Kymographs.

Considering the 50 mm $\times 4$ mm capillary as a 1-dimensional reactor, we first averaged the corrected fluorescence images over the width of the capillary (along the $y$ axis). The kymographs were obtained by stacking these profiles over time.

Front position and width.

The profiles averaged along $y$ were further averaged along the $x$ axis by performing a moving average over 25 pixels and subsequently normalized between 0 and 1. A sigmoidal function $f(x)=\frac{1}{1+e^{(-x+x_{0})/\lambda}}$ was fitted to these data. The front position corresponds to $x_{0}$ while the width of the front is defined as $\lambda$. When the sigmoidal fit failed, the position of the fronts was determined by finding the $x$ coordinate for which the value of the normalized profile was closest to 0.5.

Front velocity.

The data of front position over time was fitted by a polynomial function to reduce noise. The velocity was calculated as the time derivative of this polynomial fit.

Morphogen concentration at the front position.

The morphogen concentration was either obtained directly by measuring the fluorescence of species TH (Figure 2 of the Main Text) or indirectly by measuring the fluorescence of the cascade blue-dextran gradient (Figures 1 and 3 of the Main Text). In all cases, two control channels containing a 0 and constant concentration of the fluorescent species were used in each experiment to calibrate the intensity to concentration response. The fluorescence signal of the morphogen was shifted to 0 where the concentration of the morphogen was expected to be 0 nM (i. e. at $x=0$) and then converted to concentration by using the intensity to concentration relation previously determined.

2 Data associated to Figure 1 in the Main Text

2.1 The 1-species bistable in well-mixed conditions

2.2 Characterization of the morphogen gradient and its variation over time

During the 10 to 40 hours-long experiments the morphogen gradient diffuses slightly. Here, we estimate an order of magnitude of the variation of the gradient over time. We have

[TABLE]

where $M$ is the concentration of morphogen, $D$ its diffusion coefficient, $x$ the spatial coordinate along the channel and $t$ the time. To get an upper limit to the change of $M(x,t)$, we take the second derivative at $t=0$,

[TABLE]

If $M(x,0)=M_{max}e^{(x-x_{0})/l}$, as we show in Supplementary Figure 5b, we have

[TABLE]

We can now integrate the differential equation over time $\Delta t$ and we get

[TABLE]

where $\Delta M(x)=M(x,t)-M(x,0)$. Taking $D=10^{-2}$ mm2/min [4], $l\approx 10$ mm and $\Delta t=10^{3}$ min (16 h) we have

[TABLE]

In conclusion, the maximum concentration change over 16 h is lower than 10%.

To test this estimation experimentally, Supplementary Figure 5f shows the profile of a gradient of fluorescein-labeled species R1 recorded for 50 h. Over 50 h the spatial drift in the middle of the channel was 2 mm and the concentration drift was 10 %. The maximum concentration drift (observed on the right side of the capillary) was 30%, in agreement with Equation 5.

2.3 Data on the emergence of the parasite

In Supplementary Figure 6 a very high increase of fluorescence appears on the kymograph after 2000 min of experiment. We call the species, or mixture of species, associated to this high fluorescence a ’parasite’. Such a ’parasite’ is intrinsic to the isothermal exponential amplification reaction scheme (EXPAR) —that couples a polymerase and a nicking enzyme and that is at the core of the PEN DNA toolbox— and it corresponds to a late-phase non-specific amplification [5]. Indeed, even in the absence of an autocatalytic template, the combination of a polymerase, a nicking enzyme and dNTPs, results in the emergence of one or several autocatalytic DNA sequences that are highly fluorescent. Sequencing of parasite species has revealed that these species are mainly poly(AT) sequences with occasional T/G and A/C substitutions, and occasional palyndromic stretches containing the nicking enzyme recognition sequence [5].

We used denaturing PAGE-analysis (12.5%-PAA gel, 19:1, 50% urea, TBE-buffer, SYBR Gold staining) to follow the time-evolution of the amplification reaction of A3 in the presence of T3 (Supplementary Figure 7).

2.4 The autocatalyst template also acts as a bifurcation parameter

3 Data associated to Figure 2 in the Main Text

3.1 Spatio-temporal simulations in Figure 2

Following [7] the 1-dimensional spatio-temporal dynamics of the 2-species network were simulated with a four variable model for species H, K, IH and $I_{K}$, the first two being the autocatalysts and the last two, respectively, the inhibitors of H and K.

[TABLE]

with initial conditions

[TABLE]

and null Neumann boundary conditions.

$M(x)=30e^{-(x-20)/6}+0.1$ is the morphogen gradient, $k_{p}$ is the rate constant for production, $k_{deg}$ and $k^{i}_{deg}$ the degradation rate constants for autocatalysts and inhibitors, respectively, $\lambda$ and $\lambda_{I}$ the factors accounting for the inhibition of the species by its activator and its inhibitor, respectively. Although mathematically similar, these two terms have different chemical origin. An excess of activator (e.g. K) for a given species (e.g. IH) increases the hybridization on the corresponding template strand, which is the substrate of the polymerase and saturates it. An excess of inhibitor (e.g. IK) for a given autocatalytic species (e.g. K) increases its hybridization on the corresponding template and prevents the hybridization of the activation species (see Supplementary Figure 11 for more details on the mechanism). $D$ is the diffusion coefficient for all species. $\epsilon$ is a constant accounting for the leak reaction that generates an autocatalytic species in the absence of its activator.

Equations (6–7) with the parameters detailed in Supplementary Table 3 were solved in Mathematica using the method of lines, space was discretized with a minimum of 400 points.

To determine the parameters in Supplementary Table 3 we fitted Equation 6 with $D=0$ to a series of independent experiments in well-mixed conditions (Supplementary Figure 11b) where the concentration of TH (i.e. $M$ in Equation 6) changed. The results of the fit are shown in Supplementary Figure 12. Due to the simplicity of the model only some features of the dynamics for $t<100$ min are captured. Notably the TH concentration value at which bifurcation occurs and the non-monotonous kinetics close to this bifurcation point at $T_{H}=7-9$ nM. The simple model is in very good agreement with the data for $t>100$ min.

3.2 Reproducibility of the pattern generated by the 2-species bistable network

The results of three experiments with the 2-species bistable (network III) are presented to illustrate: i) their high reproducibility and ii) that shifting the morphogen gradient along $x$ shifts the position of the front by the same distance. We generated morphogen gradients of similar curvature but shifted along the $x$ axis (Supplementary Figure 13) for identical compositions of the reaction network. The fluorescent profiles due to species H and K were identical except for a lateral shift. The dynamics of the position of all three fronts in the morphogen concentration domain are very similar in all three experiments (Supplementary Figure 14).

3.3 Immobile morphogen gradients attached to the surface

We tested the possibility of immobilizing the morphogen gradient by attaching it to the surface using biotinylated DNA strands. Functionalization of the glass substrate was done as follows:

Chemicals.

50 mL of concentrated H2SO4, 50 mL of 30% H2O2, 100 $\mu$L of 100% CH3COOH, 150 mg of biotin-PEG-silane (Laysan Bio Inc, USA), toluene, propan-2-ol.

Protocol.

Prepare the ’piranha’ solution by slowly adding 50 ml of 30% H2O2 to concentrated H2SO4 with constant mixing. Note ’piranha’, being both strongly acidic and oxidizing, violently reacts with organics, and is thus dangerous. 2. 2.

Wipe the surface of $2.5\times 7.5$ mm2 glass slides with propan-2-ol or acetone to get rid of majority of organics, put the glass slides into a slides box and pour 100 ml of freshly prepared ’piranha’ solution. Keep for 1h. 3. 3.

Rinse the slides by sonicating them in water several times. 4. 4.

Prepare 100 mL of 0.44 mM solution of biotin-PEG-silane in toluene, add 100 $\mu$L of CH3COOH and mix. Pour this mix on the glass slides to cover them and keep for 2 hours. 5. 5.

Rinse with toluene, then with propan-2-ol and finally with water and dry with nitrogen. 6. 6.

The device finally consists of the treated glass-PEG-biotin slide, a double layer of parafilm cutted using a Graphtec CE6000 plotter cutter to define the channel geometry and a polystyrene cover slide with embedded access holes (27). It was thermally assembled at 75°C for 5 minutes. Typical channels are $30\times 2\times 0.25$ mm3, which corresponds to a volume of 15 $\mu$L. 7. 7.

Before preparing the morphogen gradient, the channels were treated with streptavidin by filling them with a 1 $\mu$M solution of streptavidin and incubated at room temperature for 10 min before washing. Morphogen gradient profiles are then prepared by using 500 nM solution of TH. 10 min incubation at room temperature is then necessary for DNA binding. To perform the experiment, the channels are rinsed and filled with the network solution.

Immobile morphogen gradients generated Polish flag patterns with network III (Supplementary Figure 15). However, it was difficult to obtain reproducible gradients on the surface. For this reason and because gradients in solution were very robust they were preferred throughout this work.

4 Data associated to Figure 3 in the Main Text

4.1 Bifurcation of network IVb in a well-mixed reactor

To be able to observe the extent of reaction of each node independently in temporal kinetic experiments we modified slightly the network. In Figure 3A of the Main Text we used network IVa (with A2 growing on non-fluorescently labeled T${}_{A_{2}}$) while here we used network IVb (with A1 growing on T${}_{A_{1}}$ fluorescently labeled with Cy3.5). A1 and A2 only differ by one base.

4.2 The position of the borders of the French flag pattern can be independently controlled

5 Data associated to Figure 4 in the Main Text

5.1 Polish flag with B2 beads

5.2 Titration of the aggregation of the two pairs of beads B1 and B2 with their linkers L1 and L2

5.3 Orthogonality of the linkers

5.4 French flag network with beads in well-mixed conditions

5.5 Materialization of a French flag pattern with conditional B1 and B2 particle aggregation

5.6 Long-term stability of the aggregated bead pattern

6 Data associated to the Discussion in the Main Text

6.1 Dependence of the front width on $D$ and on the thresholded degradation kinetics

Dimensional analysis indicates that a reaction-diffusion process governed by a diffusion coefficient $D$ and a first-order kinetic rate $k$ should generate a pattern with a characteristic length $\lambda\sim\sqrt{D/k}$. We thus expect that the width of the Polish flag pattern is given by $\lambda$.

To illustrate this intuition we used a simple one-variable model of the 1-species bistable[8]

[TABLE]

where $A(x,t)$ is the concentration of the autocatalytic species that is produced with rate $k_{p}$ and saturation $K_{p}$, repressed by the gradient of repressor $R(x)$ at rate $k_{d}$ with saturation $K_{d}$ and degraded directly by the exonuclease with rate $k_{r}$ and saturation $K_{r}$. We integrated this adimensional equation with $k_{p}=1$, $k_{r}=30$, $K_{p}=1$, $K_{d}=0.1$ and $K_{r}=100$ for $x\in[0,10]$ and for different values of $k_{d}$ and $D$ in the presence of a gradient $R(x)=10e^{(x-10)}$. We obtained Polish flag-type patterns at steady state which profiles are shown in Supplementary Figure 27AB. We extracted the width of these profiles by fitting a decaying exponential $e^{-(x-x_{0})/\lambda}$. As expected from the scaling law, $\lambda^{2}$ depended linearly on both $D$ and $1/k_{d}$ (Supplementary Figure 27C-D).

6.2 Degradation kinetics of species W2 in the presence of R2 and exonuclease

7 Supplementary movies

7.1 Movie S1 - A DNA-based network with two self-activating nodes that repress each other generates two immobile fronts that repel each other

7.2 Movie S2 - Synthesis of a reaction-diffusion French Flag pattern