Numerical simulation of stochastic gene circuits

Abstract

Armed with increasingly fast supercomputers and greater knowledge of the molecular mechanisms of gene expression, it is now practical to numerically simulate complex networks of regulated biological reactions, or gene circuits. Using an exact stochastic simulation algorithm, we obtain an accurate time-evolution of the behavior of complex gene circuits, including the effects of fluctuations caused by highly dilute, but significant, regulatory proteins. Specifically, we examine an important gene circuit, the bistable switch, and use the stochastic simulation algorithm to develop design principles that will enable us to produce a fast and robust switch for use in future applications. We pay particular attention to different transcriptional control mechanisms and their effects on the amount of fluctuations in the system, emphasizing methods that increase certainty and improve the switching rate.

Introduction

Gene expression is the primary method that a biological organism converts the stored information within its DNA into functional molecules. The dynamic regulation of gene expression is critical to an organism's proper functioning because many of the encoded functional molecules, whether they are structural proteins, enzymes, regulatory proteins, or functional RNAs, have positive effects only when certain environmental stimuli are present or at specific points in the organism's life cycle. The production of unneeded molecules comes at significant cost to the energy reservoir of the organism, in addition to any possible negative effects the molecules may cause when untimely expressed. Regulation may occur at many of the different steps of gene expression, including transcriptional initiation, elongation, and termination (Borukhov & Nudler, 2003; Burgess & Anthony, 2001; Coulombe & Burton, 1999; Dai & Rothman-Denes, 1999; Orphanides & Reinberg, 2002) as well as translational initiation and termination (Ramakrishnan, 2002; Romby & Springer, 2003). Regulatory proteins or RNA may indirectly or directly alter the binding interactions of specific sequences of DNA or mRNA, either positively or negatively affecting transcriptional and translational efficiency (Franch & Gerdes, 2000; Masse, Majdalani, & Gottesman, 2003; McLeod & Johnson, 2001; Rojo, 2001; Xu & Hoover, 2001). The stabilities of regulatory proteins and mRNAs are also regulated by other proteins, adding an additional layer of complexity (Jenal & Hengge-Aronis, 2003).

Using quantitative modeling, it is now practical to simulate large-scale systems of biological processes, including gene expression, in order to design organisms that possess specific, desired functions. By inserting a system of novel genes into an organism with a desired connectivity, one may create a variety of interesting mini-functions, or motifs.

Some proposed motifs include switches, oscillators, amplitude filters, noise filters or amplifiers, and memory storage (Wolf & Arkin, 2003). Combining these individual motifs, one may create potent applications that have significant potential in both the industrial and medical fields. One significant application is a biosensor, which may be used for the detection of both environmental toxins and disease-indicative chemical species within the human body. The ability to “program” a biological organism to perform specific functions, either with a desired frequency or in response to environmental stimuli, combines the capability of an analog signal processor with the organism's innate ability to chemically and physically alter itself and its environment. The construction of a designed system of genes that performs a specific function has been termed genetic circuit engineering (Hasty, McMillen, & Collins, 2002) with each motif considered a gene circuit. Currently, many of the proposed designs for gene circuits remain inarticulate, lacking the detail to be directly applied to an experimental system. However, with the use of highly descriptive models and relatively new mathematical formulations, it will be possible to create and design novel gene circuits in silico, guiding their successful construction in an organism.

Two functional gene circuits, the bistable switch and the oscillator, have already been constructed inside biological organisms, demonstrating the ability of multiple genes, through their interconnections, to perform specific tasks (Elowitz & Leibler, 2000; Gardner, Cantor, & Collins, 2000). The quality of the designs, however, may be improved. These gene circuits utilize genetic components from both the lac operon of E. coli and the lambda phage virus (Ptashne, 1996). Genetic components may be considered “building blocks” because of the relative ease in which they may be extracted, replicated, altered, and spliced into new biological organisms. As more systems become as well characterized as the lac operon and lambda phage virus, more building blocks will become available for use, making the design of gene circuits easier and more flexible. In addition, through common genetic engineering techniques, already available genetic components may be altered to fit the necessities of a particular design. Novel molecules may also be designed to fit a specific requirement (Ansari, Mapp, Nguyen, Dervan, & Ptashne, 2001).

The lac operon of E. coli has been well studied (Beckwith & Zipser, 1970; Kercher, Lu, & Lewis, 1997; McKnight & Yamamoto, 1992), including quantitative measurements. In this work, the genetic components, or “building blocks”, are extracted from the lac operon with the corresponding kinetic parameters, and are used in order to construct a highly descriptive representational model of gene expression. The main goal of our work is to use the model of gene expression to construct a design for a high-certainty bistable switch genetic circuit, consisting of two genes, each producing a repressor that represses the expression of the other gene (see Fig. 1). The addition of an inducer, a chemical species that disrupts the repression of a gene, toggles the switch from one state to another.

The major components of the lac operon consist of a promoter sequence, three different operator sites, and three coding genes (see Fig. 2). The Lac repressor, a symmetric tetramer protein, binds specifically to the three operator sites with different affinities and, because of its symmetry, may bind to two operators at once, looping the DNA around itself (Friedman, Fischmann, & Steitz, 1995). Repressors bound to an operator site become an obstacle to successful transcriptional initiation and elongation. Lactose, the inducer for the lac operon, or its more potent imitator isopropyl-thio-galactoside (IPTG) binds to the Lac repressor and disrupts its ability to bind to the operators in a process called induction. Induction results in a rapid increase in gene expression.

Biological systems are unique because they feature extremely small length scales, highly dilute chemical concentrations, and event driven processes. A single cell may contain only 10 molecules of a chemical species, making reactions that affect it probabilistic rather than deterministic. When such small numbers of chemical species interact, through reactions, we may no longer assume that the rates of the reactions are continuous, deterministic functions. We must represent the system as a Jump Markov stochastic process, where each reaction represents the transition of the system from one state to another, one where the reactant molecules are consumed and the product molecules are synthesized (for an introduction, see Gillespie, 1992). Reactions will occur discontinuously in time, with their probability of occurrence a function of the number of chemical molecules in the system. The stochastic simulation algorithm, an exact method of simulating a Jump Markov process, was developed by Gillespie (1976) and has been successfully used to represent the dynamics of gene expression (Arkin, Ross, & McAdams, 1998; Wolf & Arkin, 2002).

In the following sections, we define a genetic circuit bistable switch, list the reactions we use to model its gene expression and regulatory interactions, describe the kinetic parameters used, and explore the effects of changes in design on the performance output of the switch. These results lead us to develop a set of design rules for genetically engineering any high-certainty, fast and robust bistable genetic switch.