Engineering Gene Networks to Emulate Drosophila Embryonic Pattern Formation

Pattern formation is essential in the development of higher eukaryotes. For example, in the Drosophila embryo, maternal morphogen gradients establish gap gene expression domain patterning along the anterior-posterior axis, through linkage with an elaborate gene network. To understand the evolution and behaviour of such systems better, it is important to establish the minimal determinants required for patterning. We have therefore engineered artificial transcription-translation networks that generate simple patterns, crudely analogous to the Drosophila gap gene system. The Drosophila syncytium was modelled using DNA-coated paramagnetic beads fixed by magnets in an artificial chamber, forming a gene expression network. Transient expression domain patterns were generated using various levels of network connectivity. Generally, adding more transcription repression interactions increased the “sharpness” of the pattern while reducing overall expression levels. An accompanying computer model for our system allowed us to search for parameter sets compatible with patterning. While it is clear that the Drosophila embryo is far more complex than our simplified model, several features of interest emerge. For example, the model suggests that simple diffusion may be too rapid for Drosophila-scale patterning, implying that sublocalisation, or “trapping,” is required. Second, we find that for pattern formation to occur under the conditions of our in vitro reaction-diffusion system, the activator molecules must propagate faster than the inhibitors. Third, adding controlled protease degradation to the system stabilizes pattern formation over time. We have reconstituted transcriptional pattern formation from purified substances, including phage RNA polymerases, ribonucleotides, and an eukaryotic translation extract. We anticipate that the system described here will be generally applicable to the study of any biological network with a spatial component.

separate file (text file: simul.pl) from which it can be run using any machine with Perl installed. We also provide the parameter sets (param.txt) used for generating data for In Appendix A, quantitated data are provided for the Western blots in Figure 2 of the main text. In Appendix B, step-by-step instructions are given on how to build a transcription-translation chamber. Step-by-step construction of a transcription-translation chamber reaction . 15 Implementation 1

. General system description
The gene network consists of X and Y activator proteins (representing T7 and SP6 in the experimental system), and 3 target genes coding for repressor proteins A, B and C. The time-dependent distribution of all activator proteins and gene products is simulated using a system of difference equations, one for each relevant species (ie all mobile species with a heterogeneous spatial distribution as found in column 2 of Table S1).

Spatio-temporal resolution of the system
For the simulation of the experimental system, each of the 9 slabs composing the 'egg' chamber was divided into a fixed number of minislabs (typically 10) to get a higher spatial resolution of the peaks of protein forming as the system evolves, compared to that of the experimental data ( 1 per slab; see Fig 1 in main text).
The temporal resolution was always chosen to avoid numerical instabilities that are likely to arise when the following inequality does not hold : ∆t<4*∆z^2/D, where D is the effective diffusion constant for the species considered, ∆t is the timestep, and ∆z is the unit space (Euler's method; Institutionum Calculi Integralis, 1768)

System parameters
System parameters include initial spatial distributions of species, first-order rates for transcription, translation and degradation processes, matrix of binding constants for transcription factor-promoter complexes. For a complete description of parameters including default values for the repressed network, see Table S1.
We determined the binding affinities of the 3-finger DNA-binding subdomains of the repressors (proteins A, B and C), both by ELISA and gel retardation assays (Nat Biotechnol. 2001 Jul;19(7):656-60). For example, the 3-finger DNA-binding domain of Repressor C (previously called HIV-F in the Nature Biotech. paper) has an apparent Kd = 13 nM (Affinity = 77 uM-1). This value is used in the model. However, we merely estimated the binding affinities of the two six-finger fusions, Repressors A and B, based on the affinities of their component 3-finger domains and on our previous experience in

Simulation equations
For each species E, the dynamics of its concentration at time t and position s are governed by the following classical difference equation: where the terms diff, prod and deg ( accounting respectively for diffusion, production and degradation of each species) depend on the behaviour of the species. These are shown explicitly in Table S2, below. Table S2: Difference equations for each species in the repressed network ktl*Am(s,t) kdegp*E(s,t) Bp " ktl*Bm(s,t) " Cp " ktl*Cm(s,t) " # terms valid as long as s corresponds to an inner position. For boundary positions, see Section 5.

Diffusion at boundaries and T7 diffusion
At the two boundary positions of the system (first and last minislabs), the diffusion term diff, from equation 1, is modified to account for the mobility constraint imposed on the species by the presence of a wall at each end.
Thus, the diffusion factor E(s+1,t)+E(s-1,t)-2*E(s,t) used for all inner slabs (see column 1 of Table 2) is replaced by E(s+1,t)-E(s,t) for the left end position (as defined by s=0), and E(s-1,t)-E(s,t) for the right end position (as defined by the $slablength parameter). Note that this applies for all species in the system. In adjusting the model to the experimental system and accounting for the fact that T7 polymerase is the only component added dynamically to a system that has settled, we adopted an empirical approach to modelling "flow' of T7, based on our experimental observations for apparent diffusion (see Figs. S1, S2, below).
As shown in equation 2, we substituted the constant diffusion term of T7 for an exponentially-decreasing apparent diffusion factor DX, dropping from a very high initial value (DX0) to a physiological residual value (Dres).
(2) DX = ( DX0 -Dres ) * exp( -kX * t ) + Dres The basis for this approach arose from our initial observations that the T7 gradient formation appeared too rapid to be explained by diffusion, although the distribution was reproducible (indirectly observed by protein expression in our Western blots). Therefore, to investigate this phenomenon, we carried out a series of diffusion experiments in replica chamber conditions (containing agarose / in vitro translation mix / etc.) but replacing the indirectly-observed T7 pol. activity with a directly-observed coloured species; either of bromophenol blue (~0.7 kD) or the red protein, cytochrome C (12.5 kD); note that T7 polymerase itself is 99 kD. 1 ml of saturated aqueous solution of either coloured species was injected at the chamber edge and 'diffusion' was allowed to occur ( Fig S1A).
Snapshots of the cuvettes were taken at regular time intervals using a scanner, producing images in .tif format which were processed with Matlab. The images were converted into matrices of unsigned 8-bit integers (values from 0 to 255) and after collapsing on one dimension, subtracting background and normalizing to peak values, the images gave the profiles shown in Fig S1B. Freely diffusing molecules moving away from the water/dense coloured solution interface are expected to yield a spreading Gaussian distribution of concentration values, where the following ratio: x 2 /[log(1/cnorm)*t] provides an estimate of the diffusion constant value. This reads as the slope of the linear regression on the data points with x-coordinates [log(1/cnorm)*t] and y-coordinates x 2 . Fitting our experimental profiles to simple Gaussians yielded an 'apparent diffusion constant' for each time point of our data. These values were far too high for simple diffusion (completely unphysical), but could be conveniently described by 2 parameters: an initial 'apparent diffusion' constant DX0 around 2.5e4 mm 2 /s, and an exponential decay constant k of 1e-3 s -1 , describing the slowing-down of injected species over time ( Fig  S1 C,D,E).
In our modelling of the experimental system, we introduced a third parameter: the residual 'apparent diffusion constant.' This was given the approximate physiological value of 80 mm 2 /s to account for the mobility of the injected species after settling of the system (when simple diffusion becomes a dominant form of mobility).
A question remains, namely what possible mechanism could account for the rapid, decaying 'diffusion' patterns observed? We suspected that gravity/density-differences provided at least some of the driving force to the injected species: all the species tested (including T7) appear to be denser than the chamber medium. Therefore, after injection, it appears that the species sink and 'spread sideways' along the plastic chamber floor. To test this assertion, we carried out further diffusion experiments in a modified vertical chamber, where the coloured species were injected at the bottom and only allowed to diffuse against gravity (limited through the chamber boundaries; see Fig. S2).
In contrast to the horizontal situation, the vertical experiment obtained slower 'apparent diffusion' constants (in the order of 5e3 mm 2 /s), although these were still too rapid for expected values for diffusion in water (self-diffusion constant of water molecules in water: ~2e3 mm 2 /s 1 ). One possible cause is the mixing effects of convection currents set up by our desktop scanner (surface temperature ~35˚C; ambient room temperaturẽ 22˚C), although other effects might operate as well.
In summary, we applied the experimentally-derived apparent diffusion function from the horizontal experiments ( Supplementary Fig. S1) in our computer model and,   The simulation script is provided as an accompanying file (simul.pl). The script is annotated, and may be viewed and/or edited using any text editor.

Vertical 'diffusion' scans
To run, download the file from Supporting Information and either double-click the icon of the file or call it from a dos-type window (command prompt) by typing 'perl simul.pl'. To run it with the accompanying 'param.txt' file also provided, make sure both files are in the same folder.
Before any simulation is carried out, the user is prompted first for the network model to use, the parameter values to be loaded, the output content, finally the destination of output results.

Choice of network model
As described in the paper, several network connectivities of increasing complexity were tested experimentally. The strength of interaction between network component is specified in the connectivity matrix: a 3 by 3 matrix composed of binding constants where rows correspond to proteins A,B and C while columns correspond to promoters of genes A, B and C (such that, for instance, the matrix element in the first row, second column contains the binding constant between protein A and gene B).
While running the script, the user is first prompted to choose between running a simulation of one of the three pre-defined networks below: repressed network unrepressed network mutually repressed network and a fourth option, 'custom network' where the user supplies his own connectivity matrix by storing it in a simple tabulated text file (e.g. 'mat.txt').

Runmodes: 'd', 'm', 'p' or '<filename>'
One can run a simulation using the default set of parameter values for the repressed network (see Table S1) by selecting 'd' at the prompt, as a test run.
To run a single simulation with a defined set of parameter values, one can choose the manual mode of entry by typing 'm', and entering each parameter value individually, when prompted.
Alternatively, one can generate a text file where each line contains a series of tabdelimited parameter values (in the order defined in Table S1), and save it as 'param.txt' (or any other name) in the same folder as that containing the program script. To run this mode, type either 'p', for param.txt, or the full name of the file. Sample parameter sets Note: these parameter values were used for generating results for the repressed network in Fig. 3 and 4 of the main text. To run these values, first download the accompanying 'param.txt' file in Supplementary Information. Next, double-click the 'simul.pl' file containing the script to execute (or call it from a command prompt), and follow the prompts ('p' for uploading the values from the file, and twice <ENTER> to get default outputs and output location).

Output options: contents and location
After defining the input data, the user is asked to specify which species require outputs by typing their names, separated by commas. As a reminder to the user, a complete list of variables is provided. The default setting for this is selected by pressing <ENTER> and generates outputs for repressor proteins A, B and C (equivalent to typing 'Ap,Bp,Cp').
The amount of data to be generated for each output species is decided by choosing the period of output. By default, only the data for the last timestep is output, giving a snapshot of the system at the end of the simulation, but the user can specify any other period.
For each species and each timestep where output should be generated, two datasets are produced: one carries the concentration values of the given species at each position in the system, the other carries a normalized value of the previous dataset (the sum of the values over all minislabs is equal to 1). This second dataset is useful in comparing profiles for sets where absolute values differ by more than an order of magnitude.
Each dataset is stored in a separate text file carrying the species' name after the type of data the file contains (normalized 'norm' or absolute ''). Within each of these text files, each line of data consists of a tab-delimited series of concentration values of the given species corresponding to increasing positions in the system (from first to last minislab). The corresponding timelapse between two lines is equal to the period of data retrieval specified by the user.
In addition, and independently of the mode used to input parameter values (manual, default or read from a file), all parameter sets used to run the simulations are copied into a file called 'param.txt' in the same location as the previous files. This file can then be edited and used to re-run with modified parameter values.
All output results (species concentrations and parameter file) are placed in a single common folder for each time the script is run. The user is prompted to give a name to the folder, which is then placed in the same location as the perl script.

Note: For OS X, the output folders created may have insufficient privileges. Make sure folders have 'read and write' privileges from "Show Info" (cmd-I).
However, if the user chooses to use a default location (by typing <ENTER>), a series of folders are created under a single 'sim-results' folder, and the output folder for this run is placed inside according to the number of parameter sets specified (single-run or multiple-runs), and the period of output (single-time or timecourse). The naming scheme is as follows: an X followed by a random number (between 0 and 10000) (see example below).

Appendix B
Step-by-step construction of a transcription-translation chamber reaction The section that follows describes how to construct a transcription-translation chamber for gene network reactions from PCR-coated magnetic beads. Reactions are carried out on sterile petridishes (a fresh one for each new experiment). The Petri dish has a reaction chamber built on its base, using adhesive hybriwell strips. This chamber is built over a plasticine-filled dish, containing magnets. It should be noted, that we have also successfully built chambers by drilling holes for magnets in an acrylic base; this alternative can provide a more robust and regular support. Chamber construction is guided by a printed template, placed over the magnets. The chamber is filled with a transcription-translation reaction, containing ultra-low melting point agarose gel. DNAcoated beads are dispensed into this mixture, over discrete magnetic positions. Finally, T7 polymerase is pipetted at either end of the chamber to begin localised transcription.

Materials:
Sterile Petri dishes treated for cell culture (Nunclon, Nalge-Nunc, Cat. #150350) pair of magnets should be added at each end of the array to avoid edge effects (edge spots are otherwise more diffuse as the magnetic field spreads outwards, away from the line of magnets).
4. Replace the acetate template over the magnet array and fix with adhesive tape. 5. Take a sterile Petri dish with lid, place it over the plasticine-magnet base and fix in place with adhesive tape. Note the final layout of dishes and chamber in the schematic figure at the beginning of this section ("Overall chamber design").
6. Construct the edges of the chamber using strips cut from adhesive Hybriwell chambers (see picture below). Two layers of strips are sufficient to make a chamber 1 mm deep. Note: all materials and surroundings (e.g. gloves, scissors, bench) should be RNase-free and is it useful to wipe clean with RNaseZAP.
Cut the Hybriwell coverslip into approximately 3 mm wide strips with scissors.
Add Hybriwell strips to form a chamber: A then B then C then D. 7. Coat paramagnetic beads with PCR DNA (one primer biotinylated) using a Dynabeads Kilobase Binder Kit. Typically, gene A was used at 800 fmol per 10 ml beads, resuspended in 8 ml water; 200 fmol gene B and 140 fmol gene C were combined with 20 ml beads, and resuspended in 20 ml water.
8. Prepare a transcription-translation mixture using a TNT wheat germ kit: 2.5 ml water 20 ml TNT wheatgerm extract 1.2 ml TNT Reaction Buffer 0.6 ml amino acid mixture (1 mM) 1.2 ml RNasin (Promega; optional) 0.5 ml SP6 polymerase. 28 ml ultra low melting point agarose solution (melt 1.5 % (w/v) in boiling water, equilibrate in a 30˚C water bath for at least 10 min, mixing well before use).
Add the agarose last, mix and go to Step 9. immediately. 9. Dispense 54 ml of transcription-translation mixture per 18 x 3 x 1 mm chamber. Immediately add resuspended (vortexed) magnetic beads by pipetting 1 or 0.5 ml volumes gently over the base of the chamber, as close to each magnetic spot as possible.
The pipetting action is similar to loading a gel.
10. After the beads are arrayed, simultaneously inject 0.5 ml of T7 polymerase at either chamber end, using an appropriately loaded 0-10 ml multipipette. This is the start point of the experiment (see picture that follows). The chamber may be incubated at room temperature (25˚C). To minimise evaporation, tissues moistened with distilled water may be placed inside the Petri dish, avoiding contact with the transcription-translation chamber itself.
11. At the chosen end point (e.g. 15 to 60 min), transfer the entire chamber to 4˚C, to gel for 35 min. The ultra-low melting point agarose will gel at temperatures below 17˚C.
12. Gel slices may be cut with a razorblade (make a single cut, following a line on the acetate template as a guide). Aspirate the gel with a P10 Gilson pipette, set to a nominal 10 ml volume. Mix the gel samples with 10 ml SDS-loading buffer and analyse by SDS-PAGE, Western-blotting and ECL, with anti-M2 FLAG antibody.
Note:The wheat germ extract is autofluorescent and prevents direct fluorescence readout (eg GFP). However, we anticipate that with E. coli extract and using fluorescein digalactoside/lacZ as a reporter it should be possible to use a fluorescent scanner to obtain output signals in situ.