text:
# Lead Paragraph

## Lead Paragraph

Many fundamental cellular processes are governed by genetic programs which employ protein-DNA interactions in regulating function. Owing to recent technological advances, it is now possible to design synthetic gene regulatory networks. While the idea of utilizing synthetic networks in a therapeutic setting is still in its infancy, the stage is set for the notion of engineered cellular control at the DNA level. Theoretically, the biochemistry of the feedback loops associated with protein-DNA interactions often leads to nonlinear equations, and the tools of nonlinear analysis become invaluable. Here we utilize a naturally occurring genetic network to elucidate the construction and design possibilities for synthetic gene regulation. Specifically, we show how the genetic circuitry of the bacteriophage $`\lambda `$ can be used to design switching and oscillating networks, and how these networks can be coupled to cellular processes. This work suggests that a genetic toolbox can be developed using modular design concepts. Such advancements could be utilized in engineered approaches to the modification or evaluation of cellular processes.

## 1 Introduction

Remarkable progress in genomic research is leading to a complete map of the building blocks of biology. Knowledge of this map is, in turn, fueling the study of gene regulation, where proteins often regulate their own production or that of other proteins in a complex web of interactions. Post-genomic research will likely center on the dissection and analysis of these complex dynamical interactions. While the notions of protein-DNA feedback loops and network complexity are not new glass1 , savageau1 , kauffman , glass2 , glass3 , savageau2 , goodwin , tyson1 , ackers , palsson , moran , reinitz , experimental advances are inducing a resurgence of interest in the quantitative description of gene regulation novak1 , novak2 , hammond , keller , mcadams2 , arkin2 , smolen , wolf , elowitz , gardner , becskei , endy , sveiczer . These advances are beginning to set the stage for a modular description of the regulatory processes underlying basic cellular function novak1 , shapiro , mcadams1 , hartwell , lauffenburger , knight . In light of nearly three decades of parallel progress in the study of complex nonlinear and stochastic processes, the project of quantitatively describing gene regulatory networks is timely.

The concept of engineering genetic networks has roots that date back nearly half a century monod , novick . It is relatively recent, however, that experimental progress has made the design and implementation of genetic networks amenable to quantitative analysis. There are two dominant reasons for constructing synthetic networks. First, simple networks represent a first step towards logical cellular control, whereby biological processes can be manipulated or monitored at the DNA level chen . Such control could have a significant impact on post-genomic biotechnology. From the construction of simple switches or oscillators, one can imagine the design of genetic code, or software, capable of performing increasingly elaborate functions bray , knight . A second complementary motivation for network construction is the scientific notion of reduced complexity; the inherently reductionist approach of decoupling a simple network from its native and often complex biological setting can lead to valuable information regarding evolutionary design principles barkai .

Ultimately, we envision the implementation of synthetic networks in therapeutic applications. However, such a utilization depends on concurrent progress in efforts to uncover basic genomic and interspecies information. For example, broad applicability will only arise with detailed information regarding tissue-specific promoters, proteins, and genes. Likewise, quantitative network design is contingent on a firm understanding of cellular differentiation and fundamental processes such as transcription, translation, and protein metabolism. More crucially, delivery is a major hurdle; without identifiable cell-specific recognition molecules, there is no method for introducing a network to a specific type of cell. Since, in many regards, therapeutic applications are somewhat premature, we focus on the implementation of synthetic networks in less complicated organisms. The design of synthetic circuits and optimization of their function in bacteria, yeast, or other plant organisms should reveal nonlinear properties that can be employed as possible mechanisms of cellular control.

In this paper, we develop several models describing the dynamics of the protein concentration in small self-contained synthetic networks, and demonstrate techniques for externally controlling the dynamics. Although our results are general, as they originate from networks designed with common gene regulatory elements, we ground the discussion by considering the genetic circuitry of bacteriophage $`\lambda `$. Since the range of potentially interesting behavior is wide, we focus primarily on the concentration of the $`\lambda `$ repressor protein. We first show how bistability in the steady-state value of the repressor protein can arise from a single-gene network. We then show how an external noise source affecting protein degradation can be introduced to our model, and how the subsequent Langevin equation is analyzed by way of transforming to an equation describing the evolution of a probability function. We then obtain the steady-state mean repressor concentration by solving this equation in the long-time limit, and discuss its relationship to the magnitude of the external perturbation. This leads to a potentially useful application, whereby one utilizes the noise to construct a genetic switch. We next show how the addition of a second network component can lead to a genetic relaxation oscillator. We study the oscillator model in detail, highlighting the essential design criteria. We introduce a mechanism for coupling the oscillator to a time-varying genetic process. In the model equations, such coupling leads to a driven oscillator, and we study the resulting system in the framework of synchronization. We illustrate the utility of such driving through the construction of an amplifier for small periodic signals. Finally, we turn to the construction of a genetic toggle switch, and compare switching times for our network with those of a network constructed using negative regulation.

## 2 Background

Many processes involving cellular regulation take place at the level of gene transcription GenesVI . The very nature of cellular differentiation and role-specific interaction across cell types implicates a not yet understood order to cellular processes. Various modeling approaches have successfully described certain aspects of gene regulation in specific biological systems reinitz , novak1 , novak2 , arkin2 , endy , sveiczer , ackers , shea , crp\_cites . It is only recent, however, that designed network experiments have arisen in direct support of regulatory models elowitz , gardner , becskei . In this section, we highlight the results of these experimental studies, and set the stage for the discussion of the network designs described in this work.

For completeness, we first discuss the basic concepts of promoters and regulatory feedback loops mcclure , von\_hippel . A promoter region (or, simply, a promoter) denotes a segment of DNA where an RNA polymerase molecule will bind and subsequently transcribe a gene into an mRNA molecule. Thus, one speaks of a promoter as driving the transcription of a specific gene. Transcription begins downstream from the promoter at a particular sequence of DNA that is recognized by the polymerase as the start site of transcription. A chemical sequence of DNA known as the start codon codes for the region of the gene that is converted into amino acids, the protein building blocks. Feedback arises when the translated protein is capable of interacting with the promoter that drives its own production or promoters of other genes. Such transcriptional regulation is the typical method utilized by cells in controlling expression jacob , dickson , and it can occur in a positive or negative sense. Positive regulation, or activation, occurs when a protein increases transcription through biochemical reactions that enhance polymerase binding at the promoter region. Negative regulation, or repression, involves the blocking of polymerase binding at the promoter region. Proteins commonly exist as multi-subunits or multimers which perform regulatory functions throughout the cell or serve as DNA-binding proteins. Typically, protein homodimers (or heterodimers) regulate transcription, and this fact is responsible for much of the nonlinearity that arises in genetic networks smolen .

Recently, there have been three important experimental studies involving the design of synthetic genetic networks. All three employ the use of repressive promoters. In order of increasing complexity, they consist of (i) a single autorepressive promoter utilized to demonstrate the interplay between negative feedback and internal noise becskei , (ii) two repressive promoters used to construct a genetic toggle switch gardner , and (iii) three repressive promoters employed to exhibit sustained oscillations elowitz . We now briefly review the key findings in these three studies.

In the single gene study, both a negatively controlled and an unregulated promoter were utilized to study the effect of regulation on variations in cellular protein concentration becskei . The central result is that negative feedback decreases the cell-to-cell fluctuations in protein concentration measurements. Although the theoretical notion of network-induced decreased variability is not new savageau4 , this study empirically demonstrates the phenomenon through the measurement of protein fluorescence distributions over a population of cells. The findings show that, for a repressive network, the fluorescence distribution is significantly tightened, and that such tightening is proportional to the degree to which the promoter is negatively controlled. These results suggest that negative feedback is utilized in cellular design as a means for mitigating variations in cellular protein concentrations. Since the number of proteins per cell is typically small, internal noise is thought to be an important issue, and this study speaks to issues regarding the reliability of cellular processes in the presence of internal noise.

The toggle switch involves a network where each of two proteins negatively regulates the synthesis of the other; protein “A” turns off the promoter for gene “B”, and protein B turns off the promoter for gene A gardner . In this work, it is shown how certain biochemical parameters lead to two stable steady states, with either a high concentration of A (low B), or a high concentration of B (low A). Reliable switching between states is induced through the transient introduction of either a chemical or thermal stimulus, and shown to be significantly sharper than for that of a network designed without co-repression. Additionally, the change in fluorescence distributions during the switching process suggests interesting statistical properties regarding internal noise. These results demonstrate that synthetic toggle switches can be designed and utilized in a cellular environment. Co-repressive switches have long been proposed as a common regulatory theme monod2 , and the synthetic toggle serves as a model system in which to study such networks.

In the oscillator study, three repressible promoters were used to construct a network capable of producing temporal oscillations in the concentrations of cellular proteins elowitz . The regulatory network was designed with cyclic repressibility; protein A turns off the promoter for gene B, protein B turns off the promoter for gene C, and protein C turns off the promoter for gene A. For certain biochemical parameters, the “repressilator” was shown to exhibit self-sustained oscillations over the entire growth phase of the host E. coli cells. Interestingly, the period of the oscillations was shown to be longer than the bacterial septation period, suggesting that cellular conditions important to the oscillator network were reliably transmitted to the progeny cells. However, significant variations in oscillatory phases and amplitudes were observed between daughter cells, and internal noise was proposed as a plausible decorrelation mechanism. These variations suggest that, in order to circumvent the effects of noise, naturally-occurring oscillators might need some additional form of control. Indeed, an important aspect of this study was its focus on the utilization of synthetic networks as tools for biological inference. In this regard, the repressilator work provides potentially valuable information pertaining to the design principles of other oscillatory systems, such as circadian clocks.

These studies represent important advances in the engineering-based methodology of synthetic network design. In all three, the experimental behavior is consistent with predictions which arise from continuum dynamical modeling. Further, theoretical models were utilized to determine design criteria, lending support to the notion of an engineering-based approach to genetic network design. These criteria included the use of strong constitutive promoters, effective transcriptional repression, cooperative protein interactions, and similar protein degradation rates. In the immediate future, the construction and analysis of a circuit containing an activating control element (i.e., a positive feedback system) appears to be a next logical step.

In this work, we present several models describing the design of synthetic networks in prokaryotic organisms. Specifically, we will utilize genetic components from the virus bacteriophage $`\lambda `$. While other quantitative studies have concentrated on the switching properties of the $`\lambda `$ phage circuitry reinitz , arkin2 , ackers , shea , we focus on its value as a parts list for designing synthetic networks. Importantly, the biochemical reactions that constitute the control of $`\lambda `$ phage are very well characterized; the fundamental biochemical reactions are understood, and the equilibrium association constants are known ackers , ptashne , meyer , johnson1 , johnson2 , ohlendorf . In its naturally occurring state, $`\lambda `$ phage infects the bacteria Escherichia coli (E. coli). Upon infection, the evolution of $`\lambda `$ phage proceeds down one of two pathways. The lysis pathway entails the viral destruction of the host, creating hundreds of phage progeny in the process. These progeny can then infect other bacteria. The lysogenous pathway involves the incorporation of the phage DNA into the host genome. In this state, the virus is able to dormantly pass on its DNA through the bacterial progeny. The extensive interest in $`\lambda `$ phage lies in its ability to perform a remarkable trick; if an E. coli cell infected with a lysogen is endangered (i.e. exposure to UV radiation), the lysogen will quickly switch to the lysis pathway and abandon the challenged host cell.

The biochemistry of the viral “abandon-ship” response is a textbook example GenesVI  of cellular regulation via a naturally-occurring genetic switch. The lytic and lysogenic states are controlled by the cro and cI genes, respectively. These genes are regulated by what are known as the $`P_{RM}`$ (cI gene) and $`P_R`$ (cro gene) promoters. They overlap in an operator region consisting of the three binding sites OR1, OR2, and OR3, and the Cro and $`\lambda `$ repressor (“repressor”, the cI product) protein actively compete for these binding sites. When the Cro protein (product of cro gene) binds to these sites, it induces lysis. When repressor binds, lysogeny is maintained and lysis suppressed. When potentially fatal DNA damage is sensed by an E. coli host, part of the cellular response is to attempt DNA repair through the activation of a protein called RecA. $`\lambda `$ phage has evolved to utilize RecA as a signal; RecA degrades the viral repressor protein and Cro subsequently assumes control of the promoter region. Once Cro is in control, lysis ensues and the switch is thrown.

## 3 Bistability in a Single-Gene Network

In this section, we develop a quantitative model describing the regulation of the $`P_{RM}`$ operator region of $`\lambda `$ phage. We envision that our system is a DNA plasmid consisting of the promoter region and cI gene.

As noted above, the promoter region contains the three operator sites known as OR1, OR2, and OR3. The basic dynamical properties of this network, along with a categorization of the biochemical reactions, are as follows. The gene cI expresses repressor (CI), which in turn dimerizes and binds to the DNA as a transcription factor. This binding can take place at one of the three binding sites OR1, OR2, or OR3. The binding affinities are such that, typically, binding proceeds sequentially; the dimer first binds to the OR1 site, then OR2, and lastly OR3 footnote\_lambda . Positive feedback arises due to the fact that downstream transcription is enhanced by binding at OR2, while binding at OR3 represses transcription, effectively turning off production and thereby constituting a negative feedback loop.

The chemical reactions describing the network are naturally divided into two categories – fast and slow. The fast reactions have rate constants of order seconds, and are therefore assumed to be in equilibrium with respect to the slow reactions, which are described by rates of order minutes. If we let $`X`$, $`X_2`$, and $`D`$ denote the repressor, repressor dimer, and DNA promoter site, respectively, then we may write the equilibrium reactions

$`X+X`$ $`\stackrel{K_1}{}`$ $`X_2`$ (1)
$`D+X_2`$ $`\stackrel{K_2}{}`$ $`D_1`$
$`D_1+X_2`$ $`\stackrel{K_3}{}`$ $`D_2D_1`$
$`D_2D_1+X_2`$ $`\stackrel{K_4}{}`$ $`D_3D_2D_1`$
where $`D_i`$ denotes dimer binding to the ORi site, and the $`K_i=k_i/k_i`$ are equilibrium constants. We let $`K_3=\sigma _1K_2`$ and $`K_4=\sigma _2K_2`$, so that $`\sigma _1`$ and $`\sigma _2`$ represent binding strengths relative to the dimer-OR1 strength.

The slow irreversible reactions are transcription and degradation. If no repressor is bound to the operator region, or if a single repressor dimer is bound to OR1, transcription proceeds at a normal unenhanced rate. If, however, a repressor dimer is bound to OR2, the binding affinity of RNA polymerase to the promoter region is enhanced, leading to an amplification of transcription. Degradation is essentially due to cell growth. We write the reactions governing these processes as

$`D+P`$ $`\stackrel{k_t}{}`$ $`D+P+nX`$ (2)
$`D_1+P`$ $`\stackrel{k_t}{}`$ $`D_1+P+nX`$
$`D_2D_1+P`$ $`\stackrel{\alpha k_t}{}`$ $`D_2D_1+P+nX`$
$`X`$ $`\stackrel{k_x}{}`$
where $`P`$ denotes the concentration of RNA polymerase, $`n`$ is the number of repressor proteins per mRNA transcript, and $`\alpha >1`$ is the degree to which transcription is enhanced by dimer occupation of OR2.

Defining concentrations as our dynamical variables, $`x=[X]`$, $`x_2=[X_2]`$, $`d_0=[D]`$, $`d_1=[D_1]`$, $`d_2=[D_2D_1]`$, and $`d_3=[D_3D_2D_1]`$, we can write a rate equation describing the evolution of the concentration of repressor,

$$\dot{x}=2k_1x^2+2k_1x_2+nk_tp_0(d_0+d_1+\alpha d_2)k_xx$$
(3)
where we assume that the concentration of RNA polymerase $`p_0`$ remains constant during time.

We next eliminate $`x_2`$ and the $`d_i`$ from Eq. (3) as follows. We utilize the fact that the reactions in Eq. (1) are fast compared to expression and degradation, and write algebraic expressions

$`x_2`$ $`=`$ $`K_1x^2`$ (4)
$`d_1`$ $`=`$ $`K_2d_0x_2=(K_1K_2)d_0x^2`$
$`d_2`$ $`=`$ $`K_3d_1x_2=\sigma _1(K_1K_2)^2d_0x^4`$
$`d_3`$ $`=`$ $`K_4d_2x_2=\sigma _1\sigma _2(K_1K_2)^3d_0x^6`$
Further, the total concentration of DNA promoter sites $`d_T`$ is constant, so that

$$md_T=d_0(1+K_1K_2x^2+\sigma _1(K_1K_2)^2x^4+\sigma _1\sigma _2(K_1K_2)^3x^6)$$
(5)
where $`m`$ is the copy number for the plasmid, i.e., the number of plasmids per cell.

We next eliminate two of the parameters by rescaling the repressor concentration $`x`$ and time. To this end, we define the dimensionless variables $`\stackrel{~}{x}=x\sqrt{K_1K_2}`$ and $`\stackrel{~}{t}=t(k_tp_0d_Tn\sqrt{K_1K_2})`$. Upon substitution into Eq. (3), we obtain

$$\dot{x}=\frac{m(1+x^2+\alpha \sigma _1x^4)}{1+x^2+\sigma _1x^4+\sigma _1\sigma _2x^6}\gamma _xx$$
(6)
where $`\gamma _x=k_x/(d_Tnk_tp_0\sqrt{K_1K_2})`$, the time derivative is with respect to $`\stackrel{~}{t}`$, and we have suppressed the overbar on $`x`$. The equilibrium constants are $`K_1=5.0\times 10^7`$M<sup>-1</sup> and $`K_2=3.3\times 10^8`$M<sup>-1</sup> ackers , ptashne , johnson1 , johnson2 , so that the transformation from the dimensionless variable x to the total concentration of repressor (monomeric and dimeric forms) is given by \[CI\]$`=(7.7x+3.0x^2)`$ nM. The scaling of time involves the parameter $`k_t`$, and since transcription and translation are actually a complex sequence of reactions, it is difficult to give this lump parameter a numerical value. However, in Ref. hasty\_mcmillen , it is shown that, by utilizing a model for the lysogenous state of the $`\lambda `$ phage, a consistency argument yields a value for the product of parameters $`(d_tnk_tp_0)=87.6`$ nM min<sup>-1</sup>. This leads to a transformation from the dimensionless time $`\stackrel{~}{t}`$ to time measured in minutes of $`t(minutes)=0.089\stackrel{~}{t}`$.

Since equations similar to Eq. (3) often arise in the modeling of genetic circuits (see Refs. focus\_refs ) of this Focus Issue), it is worth noting the specifics of its functional form. The first term on the right hand side of Eq. (6) represents production of repressor due to transcription. The even polynomials in $`x`$ occur due to dimerization and subsequent binding to the promoter region. As noted above, the $`\sigma _i`$ prefactors denote the relative affinities for dimer binding to OR1 versus that of binding to OR2 ($`\sigma _1`$) and OR3 ($`\sigma _2`$). The prefactor $`\alpha >1`$ on the $`x^4`$ term is present because transcription is enhanced when the two operator sites OR1 and OR2 are occupied ($`x^2x^2`$). The $`x^6`$ term represents the occupation of all three operator sites, and arises in the denominator because dimer occupation of OR3 inhibits polymerase binding and shuts off transcription.

For the operator region of $`\lambda `$ phage, we have $`\sigma _12`$, $`\sigma _20.08`$, and $`\alpha 11`$ ackers , ptashne , johnson1 , johnson2 , so that the parameters $`\gamma _x`$ and $`m`$ in Eq. (6) determine the steady-state concentration of repressor. The parameter $`\gamma _x`$ is directly proportional to the protein degradation rate, and in the construction of artificial networks, it can be utilized as a tunable parameter. The integer parameter $`m`$ represents the number of plasmids per cell. While this parameter is not accessible during an experiment, it is possible to design a plasmid with a given copy number, with typical values in the range of 1-100.

The nonlinearity of Eq. (6) leads to a bistable regime in the steady state concentration of repressor, and in Figure 1A we plot the steady-state concentration of repressor as a function of the parameter $`\gamma _x`$. The bistability arises as a consequence of the competition between the production of $`x`$ along with dimerization and its degradation. For certain parameter values, the initial concentration is irrelevant, but for those that more closely balance production and loss, the final concentration is determined by the initial value.

Before turning to the next section, we make one additional observation regarding the synonymous issues of the general applicability of a synthetic network and experimental measurement. In experimental situations, a Green Fluorescent Protein (GFP) is often employed as a measurement tag known as a reporter gene. This is done by inserting the gene encoding GFP adjacent to the gene of interest, so that the reporter protein is produced in tandem with the protein of interest. In the context of the formulation given above, we can generalize Eq. (6) to include the dynamics of the reporter protein,

$`\dot{x}`$ $`=`$ $`f(x)\gamma _xx`$ (7)
$`\dot{g}`$ $`=`$ $`f(x)\gamma _gg`$
where $`f(x)`$ is the nonlinear term in Eq. (6), $`\gamma _g=k_g/(d_Tnk_tp_0\sqrt{K_1K_2})`$, and the GFP concentration is scaled by the same factor as repressor ($`\stackrel{~}{g}=g\sqrt{K_1K_2}`$). In analogy with the equation for $`x`$, $`k_g`$ is the degradation rate for GFP, and we have assumed that the number of proteins per transcript $`n`$ is the same for both processes. This ability to co-transcribe two genes from the same promoter and transcribe in tandem has two important consequences. First, since proteins are typically very stable, it is often desirable to substantially increase their degradation rate in order to access some nonlinear regime elowitz , gardner . Such a high degradation rate typically will lead to a low protein concentration, and this, in turn, can induce detection problems. The utilization of a GFP-type reporter protein can help to mitigate this problem, since its degradation rate can be left at a relatively low value. Second, and perhaps more importantly, are the significant implications for the generality of designer networks; in prokaryotic organisms, any protein can be substituted for GFP and co-transcribed, so that one network design can be utilized in a myriad of situations.

## 4 A Noise-Based Protein Switch

We now focus on parameter values leading to bistability, and consider how an external noise source can be utilized to alter the production of protein. Physically, we take the dynamical variables $`x`$ and $`g`$ described above to represent the protein concentrations within a colony of cells, and consider the noise to act on many copies of this colony. In the absence of noise, each colony will evolve identically to one of the two fixed points, as discussed above. The presence of a noise source will at times modify this simple behavior, whereby colony-to-colony fluctuations can induce novel behavior.

Noise in the form of random fluctuations arises in biochemical networks in one of two ways. As discussed elsewhere in this Focus Issue arkin\_chaos , internal noise is inherent in biochemical reactions, often arising due to the relatively small numbers of reactant molecules. On the other hand, external noise originates in the random variation of one or more of the externally-set control parameters, such as the rate constants associated with a given set of reactions. If the noise source is small, its effect can often be incorporated post hoc into the rate equations. In the case of internal noise, this is done in an attempt to recapture the lost information embodied in the rate-equation approximation. But in the case of external noise, one often wishes to introduce some new phenomenon where the details of the effect are not precisely known. In either case, the governing rate equations are augmented with additive or multiplicative stochastic terms. These terms, viewed as a random perturbation to the deterministic picture, can induce various effects, most notably, switching between potential attractors (i.e., fixed points, limit cycles, chaotic attractors) horsthemke .

In previous work, the effects of coupling between an external noise source and both the basal production rate and the transcriptional enhancement process were examined hasty . Here, we analyze the effect of a noise source which alters protein degradation. Since the mathematical formulation is similar to that of Ref. hasty , our goal here is to reproduce the phenomenology of that work under different assumptions. As in Ref.hasty , we posit that the external noise effect will be small and can be treated as a random perturbation to our existing treatment; we envision that events induced will be interactions between the external noise source and the protein degradation rate, and that this will translate to a rapidly varying protein degradation embodied in the external parameters $`\gamma _x`$ and $`\gamma _g`$. In order to introduce this effect, we generalize the model of the previous section such that random fluctuations enter Eq. (7) multiplicatively,

$`\dot{x}`$ $`=`$ $`f(x)(\gamma _x\xi _x(t))x`$ (8)
$`\dot{g}`$ $`=`$ $`f(x)(\gamma _g\xi _g(t))g`$ (9)
where the $`\xi _i(t)`$ are rapidly fluctuating random terms with zero mean ($`<\xi _i(t)>=0`$). In order to encapsulate the independent random fluctuations, we make the standard requirement that the autocorrelation be “$`\delta `$-correlated”, i.e., the statistics of the $`\xi _i(t)`$ are such that $`<\xi _i(t)\xi _j(t^{})>=D\delta _{i,j}(tt^{})`$, with $`D`$ proportional to the strength of the perturbation, and we have assumed that the size of the induced fluctuations is the same for both proteins.

Since, in Eqs. (8) and (9), the reporter protein concentration $`g`$ does not couple to the equation for the repressor concentration, the qualitative behavior of the set of equations may be obtained by analyzing $`x`$. We first define a change of variables which transforms the multiplicative Langevin equation to an additive one. Letting $`x=e^z`$, Eq. (8) becomes,

$`\dot{z}`$ $`=`$ $`{\displaystyle \frac{1+e^{2z}+22e^{4z}}{e^z+e^{3z}+2e^{5z}+.16e^{7z}}}\gamma _x+\xi _x(t)`$ (10)
$``$ $`g(z)+\xi _x(t)`$
Eq. (10) can be rewritten as:

$$\dot{z}=\frac{\varphi (z)}{z}+\xi _x(t)$$
(11)
where the potential $`\varphi (z)`$ is introduced:

$$\varphi (z)=g(z)𝑑z$$
(12)
$`\varphi (z)`$ can be viewed as an “energy landscape”, whereby $`z`$ is considered the position of a particle moving in the landscape. One such landscape is plotted in Fig. 1B. Note that the stable fixed points correspond to the minima of the potential $`\varphi `$ in Fig. 1b, and the effect of the additive noise term is to cause random kicks to the particle (system state point) lying in one of these minima. On occasion, a sequence of kicks may enable the particle to escape a local minimum and reside in a new valley.

In order to analyze Eq. (11), one typically introduces the probability distribution $`P(z,t)`$, which is effectively the probability of finding the system in a state $`z`$ at time $`t`$. Then, given Eq. (11), a Fokker-Planck (FP) equation for $`P(z,t)`$ can be constructed vankampen . The steady-state solution for this equation is given by

$$P_s(z)=Ae^{\frac{2}{D}\varphi (z)}$$
(13)
where $`A`$ is a normalization constant determined by requiring the integral of $`P_s(z)`$ over all $`z`$ be unity.

Using the steady-state distribution, the steady-state mean (ssm), $`<z>_{ss}`$, is given by

$$<z>_{ss}=_0^{\mathrm{}}zAe^{\frac{2}{D}\varphi (z)}𝑑z$$
(14)
In Fig. 1C, we plot the ssm value of $`z`$ as a function of $`D`$, obtained by numerically integrating Eq. (14). It can be seen that the ssm of z increases with $`D`$, corresponding to the increasing likelihood of populating the upper state in Fig. 1B.

Figure 1C indicates that the external noise can be used to control the ssm concentration. As a candidate application, consider the following protein switch. Given parameter values leading to the landscape of Fig. 1B, we begin the switch in the “off” position by tuning the noise strength to a very low value. This will cause a high population in the lower state, and a correspondingly low value of the concentration. Then at some time later, consider pulsing the system by increasing the noise to some large value for a short period of time, followed by a decrease back to the original low value. The pulse will cause the upper state to become populated, corresponding to a concentration increase and a flipping of the switch to the “on” position. As the pulse quickly subsides, the upper state remains populated as the noise is not of sufficient strength to drive the system across either barrier (on relevant time scales). To return the switch to the off position, the upper-state population needs to be decreased to a low value. This can be achieved by applying a second noise pulse of intermediate strength. This intermediate value is chosen large enough so as to enhance transitions to the lower state, but small enough as to remain prohibitive to upper-state transitions.

Figure 1D depicts the time evolution of the switching process for noise pulses of strengths $`D=1.0`$ and $`D=0.1`$. Initially, the concentration begins at a level of $`0.4\mu `$M, corresponding to a low noise value of $`D=0.01`$. At $`40`$ minutes, a noise burst of strength $`D=1.0`$ is used to drive the concentration to a value of $`2.2\mu `$M. Following this burst, the noise is returned to its original value. At $`80`$ minutes, a second noise burst of strength $`D=0.1`$ is used to return the concentration to its original value.

## 5 A Genetic Relaxation Oscillator

The repressillator represents an impressive step towards the generation of controllable in vivo genetic oscillations. However, there were significant cell-to-cell variations, apparently arising from small molecule number fluctuations elowitz , barkai . In order to circumvent such variability, the utilization of hysteresis-based oscillations has recently been proposed barkai . In this work, it was shown how a model circadian network can oscillate reliably in the presence of internal noise. In this section, we describe an implementation of such an oscillator, based on the repressor network of Section 3.

The hysteretic effect in Fig. 1A can be employed to induce oscillations, provided we can couple the network to a slow subsystem that effectively drives the parameter $`\gamma _x`$. This can be done by inserting a repressor protease under the control of a separate $`P_{RM}`$ promoter region. The network is depicted in Fig. 2A. On one plasmid, we have the network of Section 3; the repressor protein CI, which is under the control of the promoter $`P_{RM}`$, stimulates its own production at low concentrations and shuts off the promoter at high concentrations. On a second plasmid, we again utilize the $`P_{RM}`$ promoter region, but here we insert the gene encoding the protein RcsA. The crucial interaction is between RcsA and CI; RcsA is a protease for repressor, effectively inactivating its ability to control the $`P_{RM}`$ promoter region rozanov .

The equations governing this network can be deduced from Eq.(6) by noting the following. First, both RcsA and repressor are under the control of the same promoter, so that the functional form of the production term $`f(x)`$ in Eq.(6) will be the same for both proteins. Second, we envision our network as being constructed from two plasmids – one for repressor and one for RcsA, and that we have control over the number of plasmids per cell (copy number) of each type. Lastly, the interaction of the RcsA and repressor proteins leads to the degradation of repressor. Putting these facts together, and letting $`y`$ denote the concentration of RcsA, we have

$`\dot{x}`$ $`=`$ $`m_xf(x)\gamma _xx\gamma _{xy}xy`$ (15)
$`=`$ $`m_xf(x)\gamma (y)x`$
$`\dot{y}`$ $`=`$ $`m_yf(x)\gamma _yy`$
where $`\gamma (y)\gamma _x+\gamma _{xy}y`$, and $`m_x`$ and $`m_y`$ denote the plasmid copy numbers for the two species.

In Fig. 2B, we present simulation results for the concentration of repressor as a function time. The nature of the oscillations can be understood using Fig. 1A. Suppose we begin with a parameter value of $`\gamma (y)=4`$ and on the upper branch of the figure. The large value of repressor will then serve to activate the promoter for the RcsA, and thus lead to its increased production. An increase in the RcsA acts as an additive degradation term for repressor (see Eq. 15), and thus effectively induces slow motion to the right on the upper branch of Fig. 1A. This motion will continue until the repressor concentration falls off the upper branch at $`\gamma (y)5.8`$. At this point, with the repressor concentration at a very low value, the promoters are essentially turned off. Then, as RcsA begins to degrade, the repressor concentration slowly moves to left along the lower branch of Fig. 1A, until it encounters the bifurcation point at $`\gamma (y)3.6`$. It then jumps to its original high value, with the entire process repeating and producing the oscillations in Fig. 2B.

The oscillations in Fig. 2B are for specific parameter values; of course, not all choices of parameters will lead to oscillations. The clarification of the specific parameter values leading to oscillations is therefore important in the design of synthetic networks elowitz . For proteins in their native state, the degradation rates $`\gamma _x`$ and $`\gamma _y`$ are very small, corresponding to the high degree of stability for most proteins. For example, a consistency argument applied to a similar model for $`\lambda `$ phage switching hasty\_mcmillen  leads to $`\gamma _x0.004`$. However, using a temperature-sensitive variety of the repressor protein, $`\gamma _x`$ can be made tunable over many orders of magnitude. Other techniques, such as SSRA tagging or titration, can be employed to increase the degradation rate for RcsA. The copy numbers $`m_x`$ and $`m_y`$ can be chosen for a particular design, and the parameter $`\gamma _{xy}`$, which measures the rate of repressor degradation by RcsA, is unknown.

In Fig. 3A, we present oscillatory regimes for Eq. (15) as a function of $`\gamma _x`$ and $`\gamma _y`$, and for two fixed values of the parameter $`\gamma _{xy}`$. We see that the oscillatory regime is larger for smaller values of the parameter $`\gamma _{xy}`$. However, the larger regime corresponds to larger values of the degradation rate for RcsA. Interestingly, if we take the native (i.e., without tuning) degradation rates to be $`\gamma _x\gamma _y0.005`$, we note that the system is naturally poised very near the oscillatory regime. In Fig. 3B, we present the oscillatory regime as a function of $`\gamma _x`$ and $`\gamma _{xy}`$, and for two fixed values of $`\gamma _y`$. The regime is increased for smaller values of $`\gamma _y`$, and, in both cases, small values of $`\gamma _x`$ are preferable. Moreover, both Figs. 3A and 3B indicate that the system will oscillate for arbitrarily small values of the repressor degradation parameter $`\gamma _x`$. In Fig. 3C we depict the oscillatory regime as a function of the copy numbers $`m_x`$ and $`m_y`$, and for fixed degradation rates. Importantly, one can adjust the periodic regime to account for the unknown parameter $`\gamma _{xy}`$. Figure 3C indicates that, for oscillations, one should choose as large a copy number as possible for the plasmid containing the repressor protein ($`m_x`$). Correspondingly, one should design the RcsA plasmid with a significantly smaller copy number $`m_y`$.

We now turn briefly to the period of the oscillations. If designed genetic oscillations are to be utilized, an important issue is the dependence of the oscillation period on the parameter values. In Fig. 4A we plot the oscillation period for our CI-RcsA network as a function of the degradation parameter $`\gamma _y`$, and for other parameter values corresponding to the lower wedge of Fig. 3A. We observe that an increase in $`\gamma _y`$ will decrease the period of oscillations. Further, since the cell-division period for E. coli is $`3540`$ minutes, we note that the lower limit roughly corresponds to this period, and that, at the upper limit, we can expect four oscillations per cell division. The utilization of tuning the period of the oscillations to the cell-division time will be discussed in the next section. In Fig. 4B, we plot the period as a function of the copy number $`m_x`$. We observe that the period depends very weakly on the copy number.

## 6 Driving the Oscillator

We next turn to the utilization of an intrinsic cellular process as a means of controlling the oscillations described in the previous section. We will first consider a network design which exhibits self-sustained oscillations (i.e., with parameters that are in one of the oscillatory regions of Fig. 3), and discuss the driving of the oscillator in the context of synchronization. As a second design, we will consider a synthetic network with parameter values near, but outside, the oscillatory boundary. In that case, we will show how resonance can lead to the induction of oscillations and amplification of a cellular signal.

We suppose that an intrinsic cellular process involves oscillations in the production of protein $`U`$, and that the concentration of $`U`$ is given by $`u=u_0\mathrm{sin}\omega t`$. In order to couple the oscillations of $`U`$ to our network, we imagine inserting the gene encoding repressor adjacent to the gene encoding $`U`$. Then, since $`U`$ is being transcribed periodically, the co-transcription of repressor will lead to an oscillating source term in Eq. (15),

$`\dot{x}`$ $`=`$ $`m_xf(x)\gamma _xx\gamma _{xy}xy+\mathrm{\Gamma }\mathrm{sin}(\omega t)`$ (16)
$`\dot{y}`$ $`=`$ $`m_yf(x)\gamma _yy`$
We first consider parameter values as in Fig. 2, so that the concentrations $`x`$ and $`y`$ oscillate in the absence of driving. Here, we are interested in how the drive affects the “internal” oscillations. Although there are many interesting properties associated with driven nonlinear equations such as Eq. (16), we focus on the conditions whereby the periodic drive can cause the dynamics to shift the internal frequency and entrain to the external drive frequency $`\omega `$. We utilize the numerical bifurcation and continuation package CONT cont  to determine the boundaries of the major resonance regions. These boundaries are depicted in the parameter-space plot of Fig. 5A, where the period of the drive is plotted versus the drive amplitude. The resonance regions form the so-called Arnold tongues, which display an increasing range of the locking period as the amplitude of drive is increased. Without the periodic drive, the period of the autonomous oscillations is equal to 14.6 minutes. As one might expect, the dominant Arnold tongue is found around this autonomous period. Within this resonance region, the period of the oscillations is entrained, and is equal to the external periodic force. The second largest region of frequency locking occurs for periods of forcing which are close to half of the period of the autonomous oscillations. As a result of the periodic driving, we observe 1:2 locking, whereby the system responds with one oscillatory cycle, while the drive has undergone two cycles. Other depicted resonant regions (3:2, 2:1, 5:2, 3:1) display significantly narrower ranges for locking periods. This suggests that higher order frequency locking will be less common and probably unstable in the presence of noise. Outside the resonance regions shown in Fig. 5A one can find a rich structure of very narrow M:N locking regions with M and N quite large, together with quasiperiodic oscillations. The order of resonances along the drive period axis is given by the Farey sequence farey , i.e. in between two resonance regions characterized by rational numbers, M<sub>1</sub>:N<sub>1</sub> and M<sub>2</sub>:N<sub>2</sub>, there is a region with ratio (M<sub>1</sub>+M<sub>2</sub>):(N<sub>1</sub>+N<sub>2</sub>).

The preceding notions correspond to the driving of genetic networks which are intrinsically oscillating. We now turn to a network designed with parameter values just outside the oscillatory region, and consider the use of resonance in the following application. Suppose there is a cellular process that depends critically on oscillations of a given amplitude. We seek a strategy for modifying the amplitude of this process if, for some reason, it is too small. For concreteness, consider a cellular process linked to the cell-division period of the host for our synthetic network. For E. coli cells at a temperature of $`37`$ degrees C, this period is of order $`3540`$ minutes. Using Figs. 3A and 4A, we can deduce parameter values that will cause a CI-RscA network to oscillate, when driven, with this period. The lower wedge of Fig. 3A implies that, for $`\gamma _{xy}=0.1`$, we should design the network with values of $`\gamma _x`$ and $`\gamma _y`$ just below the lower boundary of the wedge. Fig. 4A implies that, for $`\gamma _x=0.1`$, a choice of $`\gamma _y=0.004`$ will yield oscillations with a period close to the cell-division period. In order to stay outside the oscillatory region, we therefore choose $`\gamma _y`$ just below this value. Taken together, these choices will yield a network whereby oscillations can be induced by cellular processes related to cell division. In Fig. 6A, we plot the drive versus response amplitudes ($`\mathrm{\Gamma }`$ vs $`\mathrm{\Gamma }_x`$) obtained from numerical integration of Eq. (16). We see that, depending on the proximity to the oscillatory region, oscillations are triggered when the drive reaches some critical amplitude. In Fig. 6B, we plot the gain $`g(\mathrm{\Gamma }+\mathrm{\Gamma }_x)/\mathrm{\Gamma }`$ as a function of the drive amplitude $`\mathrm{\Gamma }`$, and observe that, for certain values of the amplitude of the drive, the network can induce a significant gain.

## 7 Harnessing the Lambda Switch

The ability to switch between multiple stable states is a critical first step towards sophisticated cellular control schemes. Nonlinearities giving rise to two stable states suggest the possibility of using these states as digital signals to be processed in cellular-level computations (see, for example, bray , knight ). One may eventually be able to produce systems in which sequences of such switching events are combined to control gene expression in complex ways. In any such application, the speed with which systems make transitions between their stable states will act as a limiting factor on the time scales at which cellular events may be controlled. In this section, we describe a bistable switch based on the mechanism used by $`\lambda `$ phage, and show that such a system offers rapid switching times.

The genetic network of $`\lambda `$ phage switches its host bacterium from the dormant lysogenous state to the lytic growth state in roughly twenty minutes ptashnebook ). As discussed in Section 2, the regulatory network implementing this exceptionally fast switch has two main features: two proteins (CI and Cro) compete directly for access to promoter sites; and one of the proteins (CI) positively regulates its own level of transcription. Here, we compare a synthetic switch based on the $`\lambda `$ phage’s switching mechanism to another two-protein switch (the toggle switch described in Ref. gardner ), and numerically show that the $`\lambda `$-like system offers a faster switching time under comparable conditions.

To implement the synthetic $`\lambda `$ switch, we use the plasmid described in Section 3, on which the $`P_{RM}`$ promoter controls the expression of the $`\lambda `$ repressor protein, CI. To this, we add a second plasmid on which the $`P_R`$ promoter is used to control the expression of Cro. The operator regions OR1, OR2, and OR3 exist on each plasmid, and both proteins are capable of binding to these regions on either of the plasmids. On the $`P_{RM}`$-promoter plasmid, transcription of CI takes place whenever there is no protein (of either type) bound to OR3; when CI is bound to OR2, the rate of CI transcription is enhanced. On the $`P_R`$-promoter plasmid, Cro is transcribed only when operator site OR3 is either clear, or has a Cro dimer bound to it; either protein being bound to either OR1 or OR2 has the effect of halting the transcription of Cro.

Letting $`y`$ represent the concentration of Cro, the competition for operator sites leads to equations of the form $`\dot{x}=f(x,y)\gamma _xx`$, $`\dot{y}=g(x,y)\gamma _yy`$. We derive the form of these equations by following the process described in Section 3. As with the CI plasmid of that section, we have Eq. 1 describing the equilibrium reactions for the binding of CI to the various operator sites. To these, we add the reactions entailing the binding of Cro, and the reactions in which both proteins are bound simultaneously to different operator sites:

$`Y+Y`$ $`\stackrel{K_3}{}`$ $`Y_2`$ (17)
$`D+Y_2`$ $`\stackrel{K_4}{}`$ $`D_3^Y`$
$`D_3^Y+Y_2`$ $`\stackrel{\beta _1K_4}{}`$ $`D_3^YD_2^Y`$
$`D_3^Y+Y_2`$ $`\stackrel{\beta _2K_4}{}`$ $`D_3^YD_1^Y`$
$`D_3^YD_2^Y+Y_2`$ $`\stackrel{\beta _3K_4}{}`$ $`D_3^YD_2^YD_1^Y`$
$`D_2^XD_1^X+Y_2`$ $`\stackrel{\beta _4K_4}{}`$ $`D_3^YD_2^XD_1^X`$
$`D_1^X+Y_2`$ $`\stackrel{\beta _5K_4}{}`$ $`D_3^YD_1^X,`$
where $`Y`$ represents the Cro monomer, and $`D_i^p`$ represents binding of protein $`p`$ to the ORi site. For the operator region of $`\lambda `$ phage, we have $`\beta _1\beta _2\beta _30.08`$, and $`\beta _4\beta _51`$ ptashne , johnson1 , johnson2

The transcriptional processes are as follows. Transcription of repressor takes place when there is no protein (of either type) bound to OR3. When repressor is bound to OR2, the rate of repressor transcription is enhanced, and Cro is transcribed only when OR3 is either vacant, or has a Cro dimer bound to it. If either repressor or Cro is bound to either OR1 or OR2, the production of Cro is halted. These processes, along with degradation, yield the following irreversible reactions,

$`D+P`$ $`\stackrel{k_{tx}}{}`$ $`D+P+n_xX`$ (18)
$`D_1^X+P`$ $`\stackrel{k_{tx}}{}`$ $`D_1^X+P+n_xX`$
$`D_2^XD_1^X+P`$ $`\stackrel{\alpha k_{tx}}{}`$ $`D_2^XD_1^X+P+n_xX`$
$`D+P`$ $`\stackrel{k_{ty}}{}`$ $`D+P+n_yY`$
$`D_3^Y+P`$ $`\stackrel{k_{ty}}{}`$ $`D_3^Y+P+n_yY`$
$`X`$ $`\stackrel{k_{dx}}{}`$
$`Y`$ $`\stackrel{k_{dy}}{}`$
Following the rate equation formulation of section 3, we obtain

$`\dot{x}`$ $`=`$ $`{\displaystyle \frac{m_x(1+x^2+\alpha \sigma _1x^4)}{Q(x,y)}}\gamma _xx`$ (19)
$`\dot{y}`$ $`=`$ $`{\displaystyle \frac{m_y\rho _y(1+y^2)}{Q(x,y)}}\gamma _yy`$
where

$$Q(x,y)=1+x^2+\sigma _1x^4+\sigma _1\sigma _2x^6+y^2+(\beta _1+\beta _2)y^4+\beta _1\beta _3y^6+\sigma _1\beta _4x^4y^2+\beta _5x^2y^2.$$

The derivatives are with respect to dimensionless time, with scaling as in Section 3; $`\stackrel{~}{t}=t(k_{tx}p_0d_Tn_x\sqrt{K_1K_2})`$, where $`k_{tx}`$ is the transcription rate constant for CI, and $`n_x`$ is the number of CI monomers per mRNA transcript. The integers $`m_x`$ and $`m_y`$ represent the plasmid copy numbers for the two species; $`\rho _y`$ is a constant related to the scaling of $`y`$ relative to $`x`$. The parameters $`\gamma _x`$ and $`\gamma _y`$ are directly proportional to the decay rates of CI and Cro, respectively; we will tune these values to cause transitions between stable states. The system exhibits bistability over a wide range of parameter values, and we plot the null-clines in Fig. 7A.

For comparison, we now consider the co-repressive toggle switch briefly reviewed in section 2 gardner . This switch uses the CI and Lac proteins, where each protein shuts off transcription from the other protein’s promoter region. The experimental design was guided by the model equations,

$`\dot{u}`$ $`=`$ $`{\displaystyle \frac{\alpha _1}{1+v^\delta }}u`$ (20)
$`\dot{v}`$ $`=`$ $`{\displaystyle \frac{\alpha _2}{1+\left[\frac{u}{\left(1+\frac{[IPTG]}{K}\right)^\eta }\right]^\mu }}kv`$
where $`u`$ and $`v`$ are dimensionless concentrations of the Lac and CI proteins, respectively, and the time derivatives are with respect to a dimensionless time: $`\tau =k_dt`$, with $`k_d=2.52h^1`$ reinitz , arkin2  being the protein decay rate. The dimensionless parameters $`\alpha _1`$, $`\alpha _2`$, $`\delta `$, and $`\mu `$ define the basic model. The CI protein used in the experiments is temperature-sensitive, changing its rate of degradation with temperaturevillaverde , lowman ; we modify the original model slightly to include the factor $`k`$, which represents a varying decay rate for the CI protein. Switching is induced either by changing $`k`$, or by adjusting the concentration of isopropyl-$`\beta `$-D-thiogalactopyranoside (IPTG); the parameters $`K=2.9618\times 10^5`$ M and $`\eta =2.0015`$ define the effect of the inducer molecule IPTG on the Lac protein. Over a wide range of parameter values, the system has two stable fixed points; the null-clines are shown in Figure 7B.

The time courses of switching between stable states in the two models are shown in Figure 8; transitions are induced by eliminating the bistability, then restoring it. The precise time course of switching from one stable state to another is determined by the way in which the model parameters are adjusted to eliminate the bistability. In each case, some parameter is increased until the system passes through a saddle node bifurcation: two stable fixed points and one unstable fixed point collapse into a single stable point. In an effort to examine the behaviour of the two systems under analogous conditions, we eliminate the bistability in every case by setting the system to a parameter value 10% past the bifurcation point.

The transitions shown in Figure 8 are generated as follows. The system begins (0-1 hour) in its default bistable state, sitting at one of the two stable fixed points. Then (1–4 hours) the bistability is eliminated (as described above), with the only remaining fixed point being such that the other protein has a high concentration. Once the concentrations have switched, the default parameters are restored and the system moves to the nearby stable fixed point (4–7 hours). Finally, the system is rendered monostable again (7–10 hours), causing another transition, followed by a period (10–11 hours) during which the bistable parameters are restored.

Under the conditions shown, the $`\lambda `$ switch model displays significantly more rapid transitions between its stable states than those seen in the model of the toggle switch. The numerical results indicate that the properties of the $`\lambda `$ switch do offer an advantage in terms of the speed of transitions, indicating that it may be fruitful to study synthetic models based on this natural system. Future analytical work on models such as the one presented in this section may allow us to make more precise statements regarding the source of this advantage.

## 8 Conclusion

From an engineering perspective, the control of cellular function through the design and manipulation of gene regulatory networks is an intriguing possibility. Current examples of potential applicability range from the use of genetically engineered microorganisms for environmental cleanup purposes szafranski , to the flipping of genetic switches in mammalian neuronal cells harding . While the experimental techniques employed in studies of this nature are certainly impressive, it is clear that reliable theoretical tools would be of enormous value. On a strictly practical level, such techniques could potentially reduce the degree of “trial-and-error” experimentation. More importantly, computational and theoretical approaches will lead to testable predictions regarding the current understanding of complex biological networks.

While other studies have centered on certain aspects of naturally-occurring genetic regulatory networks reinitz , novak1 , novak2 , arkin2 , endy , sveiczer , ackers , shea , crp\_cites , an alternative approach is to focus on the design of synthetic networks. Such an engineering-based approach has significant technological implications, and will lead, in a complementary fashion, to an enhanced understanding of biological design principles. In this work, we have shown how several synthetic networks can be designed from the genetic machinery of the virus $`\lambda `$ phage. We have highlighted some of the possible behavior of these networks through the discussion of the design of two types of switches and a relaxation oscillator. Additionally, in the case of the oscillator, we have coupled the network to an existing cellular process. Such coupling could lead to possible strategies for entraining or inducing network oscillations in cellular protein levels, and prove useful in the design of networks that interact with cellular processes that require precise timing.

With regard to model formulation, there are several intriguing areas for further work. For one, the number of molecules governing the biochemistry of genetic networks is often relatively small, leading to interesting issues involving internal noise. Recent pivotal work mcadams2 , arkin2 , mcadams1  has led to a systematic modeling approach which utilizes a Monte Carlo-type simulation of the biochemical reactions gillespie . While this approach is impressively complete, its complexity makes analysis nearly impossible. An alternative approach could entail the use of Langevin equations, whereby the effects of internal noise are incorporated into stochastic terms whose magnitudes are concentration-dependent. Indeed, in the context of genetic switches, this approach has recently been suggested bialek . The advantage of this formulation is that stochastic effects can be viewed as a perturbation to the deterministic picture, so that analytic tools can be utilized.

A potentially important technical issue involves the implicit assumption that the reactions take place in three-dimensional space. While this assumption is perhaps the most natural, proteins have been observed sliding along a DNA molecule in search of a promoter region harada , so that protein-DNA reactions might effectively take place on a surface. While this would not alter the qualitative form of Eq. (6), the exponents on the variable $`x`$ could take on other values savageau3 , and this, in turn, could lead to significant quantitative differences.

It has been nearly 30 years since the pioneering theoretical work on interacting genetic networks glass1 -tyson1 . Due, in part, to the inherent complexity of regulatory networks, the true significance of these studies had to await technological advances. Current progress in the study of both naturally occurring and synthetic genetic networks suggests that, as the pioneers envisioned, tools from nonlinear dynamics and statistical physics will play important roles in the description and manipulation of the dynamics underlying cellular control.

ACKNOWLEDGEMENTS. We warmly acknowledge insightful discussions with William Blake, Michael Elowitz, Doug Mar, John Reinitz, and John Tyson. This work is supported by The Fetzer Institute (J. H.) and the Office of Naval Research.

## Figure Captions

FIG. 1. Results for additive noise with parameter value $`m=1`$. (A) Bifurcation plot for the steady-state concentration of repressor versus the model parameter $`\gamma _x`$. (B) The energy landscape. Stable equilibrium values of Eq. (10) (with $`D=0`$) correspond to the valleys at $`z=1.6`$ and $`0.5`$, with an unstable value at $`z=0.52`$. (B) Steady-state probability distributions for noise strengths of $`D=0.04`$ (solid line) and $`D=0.4`$ (dotted line). (C) The steady-state equilibrium value of $`z`$ plotted versus noise strength. The corresponding concentration will increas as the noise causes the upper state of (B) to become increasingly populated. (D) Simulation of Eqs. (8) and (9) demonstrating the utilization of external noise for protein switching. Initially, the concentration begins at a level of $`\text{[GFP]}0.4\mu `$M corresponding to a low noise value of $`D=0.01`$. After $`40`$ minutes, a large $`2`$-minute noise pulse of strength $`D=1.0`$ is used to drive the concentration to $`2.2\mu `$M. Following this pulse, the noise is returned to its original value. At $`80`$ minutes, a smaller $`10`$-minute noise pulse of strength $`D=0.1`$ is used to return the concentration to near its original value. The simulation technique is that of Ref. sancho .
FIG. 2. The relaxation oscillator. (A) Schematic of the circuit. The $`P_{RM}`$ promoter is used on two plasmids to control the production of repressor (X) and RcsA (Y). After dimerization, repressor acts to turn on both plasmids through its interaction at $`P_{RM}`$. As its promoter is activated, RcsA concentrations rise, leading to an induced reduction of repressor. (B) Simulation of Eqs. (15). Oscillations arise as the RcsA-induced degradation of repressor causes a transversal of the hysteresis diagram in Fig. 1A. The parameter values are $`m_x=10`$, $`m_y=1`$, $`\gamma _x=0.1`$, $`\gamma _y=0.01`$, and $`\gamma _{xy}=0.1`$.
FIG. 3. Oscillatory regimes for the relaxation oscillator. (A) The bifurcation wedge is larger for smaller values of the parameter $`\gamma _{xy}`$. This larger regime corresponds to larger values of the RcsA degradation parameter $`\gamma _y`$. Note that the native (i.e., without tuning) degradation rates of $`\gamma _x\gamma _y0.005`$ are very near the oscillatory regime. (B) Bifurcation diagrams as a function of $`\gamma _x`$ and $`\gamma _{xy}`$, and for two fixed values of $`\gamma _y`$. The oscillatory regime is increased for smaller values of $`\gamma _y`$, and, in both cases, small values of $`\gamma _x`$ are preferable for oscillations. (C) The bifurcation diagram as a function of the copy numbers $`m_x`$ and $`m_y`$, and for fixed degradation rates. Importantly, one can adjust the periodic regime to account for the unknown parameter $`\gamma _{xy}`$. The figure also indicates that, for oscillations, one should choose as large a copy number as possible for the plasmid containing the repressor protein ($`m_x`$). In (A) and (B), constant parameter values are $`m_x=10`$ and $`m_y=1`$, and in (C) $`\gamma _x=1.0`$ and $`\gamma _y=0.01`$.
FIG. 4. Parameter dependence of the oscillatory period. (A) An increase in $`\gamma _y`$ decreases the period of oscillations. (B) The period depends very weakly on the copy number. In (A) $`m_x=10`$ and in (B) $`\gamma _y=0.01`$, and for both plots, other parameter values are $`\gamma _x=0.1`$, $`\gamma _xy=0.1`$, and $`m_y=1`$.
FIG. 5. Dynamics of a periodically driven relaxation oscillator - Eq.(15); $`\gamma _x=0.1`$, $`\gamma _y=0.01`$, $`\gamma _{xy}=0.1`$. (A) Resonant regions in period-amplitude parametric plane. Solid lines (limit lines of periodic solutions) together with dashed (period doubling) lines define boundaries of stable periodic solutions for a given phase locking region M:N, where M is the number of relaxation oscillation and N is the number of driving sinusoidal oscillations. (B,C,D) Oscillations in periodically driven repressor (top curve) concentration together with the oscillations of the sinusoidal driving (bottom curve). (B) 1:1 synchronization; the 14.6 minute period of cI oscillations is equal to the driving period. (C) 1:2 phase locking; the 29.2 minute period of the cI oscillations is twice as long as the driving period. (D) 2:1 phase locking; the 7.6 minute period of the cI oscillations is equal to one half of the driving period.
FIG. 6. (A) As a function of the driving amplitude $`\mathrm{\Gamma }`$, the amplitude $`\mathrm{\Gamma }_x`$ of the induced network oscillations shows a sharp increase for a critical value of the drive. The critical value corresponds to a drive large enough to induce the hysteretic oscillations, and it increases as one decreases $`\gamma _y`$ and moves away from the oscillatory region in parameter space. The three curves denoted 1,2, and 3 are for $`\gamma _y`$ values of 0.0038, 0.0036, and 0.0034. (B) The gain as a function of the drive amplitude for $`\gamma _y=0.0038`$. Close to the oscillatory region, a significant gain in the drive amplitude can be induced. Parameter values for both plots are $`\gamma _x=0.1`$, $`\gamma _{xy}=0.1`$, $`m_x=10`$, and $`m_y=1`$. Note that, corresponding to these values, the network does not oscillate (without driving) for $`\gamma _y<0.004`$ (see the bottom wedge of Fig. 3A).
FIG. 7 Null-clines for the two-protein bistable switch systems. Stable fixed points are marked with circles, and unstable fixed points are marked with squares. (A) Null-clines for the synthetic $`\lambda `$ switch, Eqs. (19). Solid line: $`\dot{x}=0`$ cline. Dashed line: $`\dot{y}=0`$ cline. Parameter values: $`\gamma _x=0.004`$; $`\gamma _y=0.008`$; $`\rho _y=62.92`$; $`\alpha =11`$; $`m_x=m_y=1`$; $`\sigma _1=2`$; $`\sigma _2=0.08`$; $`\beta _1=\beta _2=\beta _3=0.08`$; and $`\beta _4=\beta _5=1`$. (B) Null-clines for the toggle switch, Eqs. (20). Solid line: $`\dot{u}=0`$ cline. Dashed line: $`\dot{v}=0`$ cline. Parameter values (from Ref. gardner ): $`\alpha _1=156.25`$; $`\alpha _2=15.6`$; $`\delta =2.5`$; $`\mu =1`$; $`\eta =2.0015`$; $`[IPTG]=0`$; $`k=1`$.
FIG. 8. Transitions between stable states for the two-protein bistable switch systems. The protein concentrations have been normalized (the trace for each protein is normalized relative to its own maximum value). The system parameters are varied over time, altering the stability of the system and causing transitions, as described in the text. Upper Plots: The switching of Lac (solid) and CI (dashed) in the synthetic toggle model gardner . The parameter values are as given in the caption to Fig. 7, except as follows. (1–4 hours): $`[IPTG]=2`$ mM, $`k=1.0`$. (7–10 hours): $`[IPTG]=0.0`$, $`k=50.81`$. Lower Plots: The switching of CRO (dashed line) and CI (solid) in the synthetic $`\lambda `$ model. The parameter values are as given in the caption to Fig. 7, except as follows. (1–4 hours): $`\gamma _x=0.004`$, $`\gamma _y=21.6`$. (7–10 hours): $`\gamma _x=18.0`$, $`\gamma _y=0.008`$.

Figure 1 - Hasty et al.

Figure 2 - Hasty et al.

Figure 3 - Hasty et al.

Figure 4 - Hasty et al.

Figure 5 - Hasty et al.

Figure 6 - Hasty et al.

Figure 7 - Hasty et al.

Figure 8 - Hasty et al.


format: text

text:
# 1 Introduction: GRBs and the blow up of the energy

## 1 Introduction: GRBs and the blow up of the energy

Gamma Ray Burst, GRB, has been associated in last decade with huge explosions many order of magnitude ($`10^{10}`$) more power-full than common Supernova. These bursts are so intense ( because of GRB cosmic distance origin) and energetic ( at the electron mass edges) as well as so sharp (within millisecond scale in time) that their own opacity (at corresponding small hundred kilometers sizes) make difficult their own shining outside (an over Eddington luminosity). This Eddington-opacity is reached also for smaller distance (and less power-full) sources, named Soft Gamma Repeaters SGRs, because softer and sometime repeating. For all last decade ruled Fireball: a huge isotropic explosive model where boosted relativistic electron pairs sea and their annihilation may have remitted soft KeV photons into harder MeV energy band in a wide shells. The GRBs structure arises in Fireball by internal-external shell scattering. We considered since 1993 the $`\gamma `$ Jet model for GRBs and SGRs both of them at largest galactic halo; after the rare GRB 980425 and its Supernova-GRB connection, we committed to explain both of them as a precessing Jet at different energy regime, both for GRBs at Supernova beamed power and SGRs at much lower X-pulsar beamed luminosity. Most of the scientific community remained firmly on Fireball model before and after the first cosmic identification of GRB at cosmic edges on February 1997. For me in that epoch was difficult to accept huge supernova jet and I preferred a galactic halo model; since April 25 GRB, however, I immediately believed in GRB and SGR unified model at huge different beamed powers. However the fireball believers, a formidable school all over the world, increased their excitement as much as the GRB distance (or redshift) and its apparent output power was increased by event to event. Their explosive model was feed by the apparent X-decay after-glow tail. Nevertheless the huge $`GRB990123`$ event blow up too much even for the most excited mind: its unbelievable energy puzzle (imagine, two or more solar masses totally converted in pure gamma within tens of seconds within a millisecond-hundred kilometer substructures) forced most authors to bend their collimate Fireball model toward en-longed beamed Jet, in order to economize the energy. Therefore with a polite compromise the Fireball models shifted toward a wide (ten degrees or so ) open cone or a Fountain Jet Models. Still one-shoot Fireball model. How was possible to hold slimy neutrino in such a Jet is a mystery. No general agreement on the Supernova-GRB connection was in general acceptable to Fireball believers. The new Jet-Fireball shared still the explosive nature and the shock wave time modulation ruled by external relic shells, but it relaxed the energy demand by the inverse of the solid angle beam. The word beam or jet has been often hidden under new fashionable names. As the time went on the beam of Fireball’s model shrinks and the GRBs energy reduces from ($``$ $`10^{54}`$ $`erg`$) first to ($``$ $`10^{51}`$ $`erg`$) and even toward ($``$ $`10^{48}`$ $`erg`$). However the smaller the (explosive) beam angle the higher is the needed GRB event rate: from nominal $`10^6`$ event a year for galaxy one must expect a more frequent $`10^3`$ or even $`1`$ GRB a year in galaxies. Therefore we must accept the co-existence in our Universe of two totally separate wild animals: Supernova whose explosive isotropic output reaches $`10^{44}`$ $`ergs^1`$ powers and a total energy $`10^{51}`$ $`erg`$ and a as common as animal whose jet output is $`10^{48}`$ $`ergs^1`$ and total energy $`10^{51}`$ $`erg`$. A possible GRB-Supernova connection call for an additional paradox: Why GRB Fireball-Jet power is 4-7 order of magnitude higher than the corresponding optical one in Supernova. Why not to respect at least energy equi-partition? In this Fireball scenario the rare GRB of April 1998 stand outside or better is often hidden or taken away because of the more puzzling apparent $`low`$ energy output and softer nature while being nearer and un-red-shifted. Of course off-axis Fireball event make the miracle but cannot solve the statistics.
We argued (1998-1999) on the contrary that GRBs and SGRs find a comprehensive theory within a very thin (tens of seconds) spinning and multi precessing $`\gamma `$ Jet, sprayed by a Neutron Star, NS, or Black Hole, BH, (Fargion 1994-1999, Fargion, Salis 1995-1998). Higher black holes mass ( $`10^1`$ $`10^4`$ solar masses) cannot help Fireball isotropic survival because contradictions with the sharp GRB time-scales. We have shown that the same rarity of GRB-SN detection and the established GRB980425-SN link favors very clearly the thin (tens of arc second or millisecond radiant) Jet Nature of GRBs. Such thin jet may accommodate the Supernova power into an apparent GRB luminosity in agreement with energy equi-partition. Moreover the Jet in precession explains Soft Gamma Repeaters. Even originally ($`197080`$) most authors unified GRB/SGR models , since last fifteen years are commonly separated by their repeater and spectra differences and finally their much nearer galactic distances; however we and other have shown that they, rarely, openly shared the same spectra, time and flux structures. Last SGR1900+14 (May-August-October 1998) events and SGR1627-41 (June-October 1998) events did exhibit at peak intensities hard spectra comparable with classical GRBs. Indeed the SGR1900+14 event BATSE trigger 7171 left an almost identical event comparable to a just following GRB (trigger 7172) on the same day, same detector, with same spectra and comparable flux (Fargion 1998-1999) . This Hard-Soft connection has been re-discovered and confirmed more recently by BATSE group ( Woods et all 1999) with an additional hard event of SGR 1900+14 recorded in GRB990110 event.
Additional GRB-SGR connection occur between GRB980706 event with an almost identical (in time, channel spectra, morphology and intensities) observed in GRB980618 originated by SGR 1627-41. Nature would be very perverse in mimic two signals, (even if at different distances and different powers), by two extreme different source engines. This points once again to their common Nature. At least SGRs offer the unique occasion for a study of a near-by laboratory GRB to a better understanding of both models. However they are both up present times usually described by catastrophic spherical explosions, Fireball and Magnetar Models. Their different distances, cosmic versus galactic ones, imply very different power source Jet, but their morphological similarity strongly suggest an unique process: the blazing of a spinning and multi-precessing gamma Jet, from either Neutron Star or Black Hole. The $`\gamma `$ Jet is born by high GeVs electron pairs Jet which are regenerating, via Inverse Compton Scattering, an inner collimated beamed $`\gamma `$ (MeVs) precessing Jet. The thin jet (an opening angle inverse of the electron Lorentz factor, a milli-radiant or below), while spinning, is driven by a companion and/or an asymmetric accreting disk in a Quasi Periodic Oscillation (QPO) and in a Keplerian multi-precessing blazing mode: its $`\gamma X`$ ray lighthouse trembling and flashing is the source of the complex and wide structure of observed Gamma Bursts. These $`\gamma `$ Jets share a peak power of a Supernova ($`10^{44}ergs^1`$) at their birth (during SN and Neutron Star formations), decaying by power law $`t^1`$ $``$ $`t^{(1.5)}`$ to less power-full Jets that converge to present persistent SGRs stages. Indeed these ones are blazing events from late relic X pulsar observable only at nearer distances. The $`\gamma `$ Jet emit in general at $``$ $`10^{35}`$ erg$`s^1`$ powers; both of GRB and SGR show an apparent luminosity amplified by the inverse square of the thin from $`10^3`$ to $`10^4`$ radiant angle Jet beaming: the corresponding solid angle $`\mathrm{\Omega }`$ spreads between $`10^7`$ and $`10^9`$ leading to the necessary amplification for both SN-GRB and X pulsar-SGR . Optical-Radio After-Glows are not the fading fireball explosion tails often observed in puzzling variable non monotonic decay, but the averaged external Jet tails moving and precessing and geometrically fading away. The rare optical re-brightening (the so called SN bump) observed in few afterglow has been, very probably, erroneously associated to an underlying isotropic SN flash: the optical re-brightening afterglow is in general the late re-crossing of the precessing Jet tail periphery toward the observer direction. The averaged integral optical signal (within thousand of seconds) hide the short-scale oscillatory behaviour of the precessing Jet. In particular the geometrical beaming offered in the rare GRB970508 a peculiar optical re-brightening and a manifest late radio oscillating afterglow of the cycloidal lighthouse Jet.

## 2 The geometrical multi-precessing Gamma Jet and GRB bursting signature

We imagine the GRB and SGR nature as the early and the late stages of jets fueled first by SN event and later by a disk or a companion (WD, NS) star. Their binary angular velocity $`\omega _b`$ reflects the beam evolution $`\theta _1(t)=\sqrt{\theta _{1m}^2+(\omega _bt)^2}`$ or more generally a multi-precessing angle $`\theta _1(t)`$ (Fargion & Salis 1996):

$$\theta _1(t)=\sqrt{\theta _x^2+\theta _y^2}$$
(1)
$$\theta _x(t)=\theta _bsin(\omega _bt+\phi _b)+\theta _{psr}sin(\omega _{psr}t)+\theta _Nsin(\omega _Nt+\phi _N)$$
(2)
$$\theta _y(t)=\theta _{1m}+\theta _bcos(\omega _bt+\phi _b)+\theta _{psr}cos(\omega _{psr}t)+\theta _Ncos(\omega _Nt+\phi _N)$$
(3)
where $`\theta _{1m}`$ is the minimal angle impact parameter of the jet toward the observer, $`\theta _b`$, $`\theta _{psr}`$, $`\theta _N`$ are, in the order, the maximal opening precessing angles due to the binary, spinning pulsar, nutation mode of the jet axis. For a 3D pattern and its projection along the vertical axis in an orthogonal 2D plane (Fig.1) The angular velocities combined in the multi-precession keep memory of the pulsar jet spin ($`\omega _{psr}`$), the precession by the binary $`\omega _b`$ and an additional nutation due to inertial momentum anisotropies or beam-accretion disk torques ($`\omega _N`$). On average, from eq.(5) the $`\gamma `$ flux and the $`X`$ optical afterglow decays, in first approximation, as $`t^2`$; the complicated spinning and precessing jet blazing is responsible for the wide morphology of GRBs and SGRs as well as their internal periodicity.(See Fig.2). The consequent $`\gamma `$ time evolution and spectra derived in this ideal models may be compared successfully with similar observed GRB data evolution as shown in Fig. 3, regarding X-ray precursors in $`GRB971210,GRB971212,GRB990518,GRB000131`$.

Similar descriptions with more parameters and with a rapid time evolution of the jet has been more recently also proposed by (Portegies Zwart et all 1999).

## 3  X Ray precursor as the fingerprint of Precessing Gamma Jets and the end of Fireballs

To choose for a model let us just consider with no prejudice the last reported (and most distant $`z=4.5`$) event: GRB000131 and its $`X`$ ray precursor:(fig 3,last). This event while being red-shifted and slowed down by a factor 5.5 exhibit on the contrary a short scale time fine structure not explicable by any fireball model , but well compatible with a thin, fast spinning precessing $`\gamma `$ jet. The extreme $`\gamma `$ energy budget, calling for a comparable $`\nu `$ one, exceeds few solar masses in its main emission even for ideal full energy conversion. Moreover one must notice the presence of a weak $`X`$-ray precursor pulse lasting 7 sec, 62 sec before the huge main structured $`\gamma `$ burst trigger. Its arrival direction (within 12 degree error) with main GRB is consistent only with the main pulse (a probability to occur by chance below $`3.610^3`$). The time clustering proximity (one minute over a day GRB rate average) has the probability to occur by chance below once over a thousand. The over all probability to observe this precursor by change is below 3.4 over a million making inseparable its association with the main GRB000131 event. This weak burst signal correspond to a power above a million Supernova and have left no trace or Optical/X transient just a minute before the real (peak power $`>billion`$ Supernova) energetic event. Similar X precursors occurred in a non negligible minor sample of GRBs (see for example Fig 2-4a) and also few SGRs event (Fig 4b). No isotropic GRB explosive progenitor could survive such a disruptive isotropic(million supernova output) precursor trigger nor any multi-exploding jet. Only a persistent precessing Gamma Jet crossing nearby the observer direction twice could naturally explain it.

## 4  Thin Precessing $`\gamma `$ Jet and the GRB-SN Connection

Statistical arguments (Fargion 1998, 1999) favor a unified GRBs model based on blazing, spinning and precessing thin jet. We assume that GRB jet arise in most SNe outbursts. The far GRBs are observables at their peak intensities (coincident to SN) while blazing in axis to us within the thin jet very rarely; consequently the hit of the target occurs only within a wide sample of sources found in a huge cosmic volume. In this frame work the GRB rate do not differ much from the SN rate. Assuming a SN-GRB event every 30 years in a galaxy and assuming a thin angular cone ($`\mathrm{\Omega }<`$(1/4)$`10^8`$) the probability to be within the cone jet in a ($`410^{10}`$)cosmic sample of galaxies (at limited Hubble R$``$ 28 magnitude ) within our main present observable Universe volume ($`z1`$-$`z4`$) during one day of observation at a nominal 10 sec GRB duration is quite small: ($`P<10^2`$). This value should be suppressed by nearly an order of magnitude because of the detector acceptance. However a precessing gamma jet whose decaying scale time is a thousand time longer than the GRB itself (decaying by a power law $`t^1`$) whose scale time is nearly ten or twenty thousand of seconds, may fit naturally the observed GRB rate. Moreover the jet pressure could also accumulate gas and form dense filaments. Such filaments fragment as well as gravitationally cluster leading to contemporaneous stellar arc formations.(Y.Efremov&D.Fargion, 2000). This phenomena might explain why GRBs seem to be associated with star formation regions.

## 5 GRBs Energy Budget in Fireball and the Neutrino $`\nu `$ and $`\gamma `$ ejection

Gamma Ray Bursts as recent $`GRB990123`$ and $`GRB990510`$ emit, for isotropic explosions, $`\gamma `$ energies as large as two solar masses annihilation. These energies are underestimated because of the neglected role of comparable ejected MeV (Comptel signal) neutrinos $`\nu `$ bursts and assume an unrealistic ideal energy conversion efficiency. Indeed, as often neglected, it is important to remind that the huge energy bath (for a fireball model) on GRB990123 imply also a corresponding neutrino burst. As in hot universe, if entropy conservation holds, the $`\nu `$ energy density factor to be added to the photon $`\gamma `$ budget is at least $`((21/8)\times (4/11)^{4/3})`$. If entropy conservation do not hold the energy needed is at least a factor $`[21/8]`$ larger than the gamma one. The consequent total energy-mass needed for the two cases are respectively 3.5 and 7.2 solar masses. Additional factors must be introduced for realistic energy conversion efficiency leading to energies as large as tens of solar masses. No fireball by NS may coexist with it. Jet could. These extreme power cannot be explained with any standard spherically symmetric Black Hole Fireball model. A too heavy Black Hole (hundred or thousands solar masses) or, worse, Star would be unable to coexist with the shortest millisecond time structure of Gamma Ray Burst. Cosmological and nearby Gravitational Red-shifts may only make the Fireball Model more inconsistent. Smaller size BH or NS do not offer enough mass reserve. Beaming of the gamma radiation may overcome the energy puzzle along with the short scale-time. However any mild ”explosive beam” event as some models (Wang & Wheeler 1998) $`(\mathrm{\Omega }>10^2)`$ would not solve the jet containment at the corresponding disruptive energies. Moreover such a small beaming would not solve the huge GRBs flux energy windows ($`10^{47}÷10^{54}`$ erg/sec), keeping GRB980425 and GRB990123 within the same GRB framework. Only extreme beaming $`(\mathrm{\Omega }10^8)`$, by a slow decaying, but long-lived precessing jet, may coexist with characteristic Supernova energies, apparent GRBs output and the puzzling GRB980425 statistics as well as the GRB connection with older, nearer and weaker SGRs relics. Therefore SGRs are very useful nearby astrophysical Laboratory where to study and test the far GRB process. SGRs are not associated with huge OT afterglow or explosive SN event : they have lost their primordial SNR shells elsewhere while escaping at high velocity far from the SN birth place.

## 6 Hard Gamma Jet by Inverse Compton Scattering by GeV Electron Pairs beam

A jet angle related by a relativistic kinematics would imply $`\theta \frac{1}{\gamma _e}`$, where $`\gamma _e`$ is found to reach $`\gamma _e10^3÷10^4`$ (Fargion 1994, 1998). At first approximation the gamma constrains is given by Inverse Compton relation: $`<ϵ_\gamma >\gamma _e^2kT`$ for $`kT10^310^1eV`$ and $`E_eGeVs`$ leading to characteristic X-$`\gamma `$ GRB spectra. The origin of $`GeVs`$ electron pairs are very probably decayed secondary related to primary inner muon pairs jets, able to cross dense stellar target. The consequent adimensional photon number rate (Fargion & Salis 1996) as a function of the observational angle $`\theta _1`$ responsible for peak luminosity becomes

$$\frac{\left(\frac{dN_1}{dt_1d\theta _1}\right)_{\theta _1(t)}}{\left(\frac{dN_1}{dt_1d\theta _1}\right)_{\theta _1=0}}\frac{1+\gamma ^4\theta _1^4(t)}{[1+\gamma ^2\theta _1^2(t)]^4}\theta _1\frac{1}{(\theta _1)^3}.$$
(4)
The total fluence at minimal impact angle $`\theta _{1m}`$ responsible for the average luminosity is

$$\frac{dN_1}{dt_1}(\theta _{1m})_{\theta _{1m}}^{\mathrm{}}\frac{1+\gamma ^4\theta _1^4}{[1+\gamma ^2\theta _1^2]^4}\theta _1𝑑\theta _1\frac{1}{(\theta _{1m})^2}.$$
(5)
These spectra fit GRBs observed ones (Fargion & Salis 1995). Assuming a beam jet intensity $`I_1`$ comparable with maximal SN luminosity, $`I_110^{45}ergs^1`$, and replacing this value in the above a-dimensional equation we find a maximal apparent GRB power for beaming angles $`10^3÷3\times 10^5`$, $`P4\pi I_1\theta ^210^{52}÷10^{55}ergs^1`$, just within observed ones. We also assumed a power law jet time decay as follows

$$I_{jet}=I_1\left(\frac{t}{t_0}\right)^\alpha 10^{45}\left(\frac{t}{310^4s}\right)^1ergs^1$$
(6)
where ($`\alpha 1`$) is a value able to reach, at 1000 years time scales, the present known galactic micro-jet (as SS433) intensities powers: $`I_{jet}10^{39}ergs^1`$. This offer a natural link between the GRB and the SGR out-put powers. We used the model to evaluate if April precessing jet might hit us once again. It should be noted that a steady angular velocity would imply an intensity variability ($`I\theta ^2t^2`$) corresponding to some of the earliest afterglow decay law.

## 7 Precessing Jet Relics in action

The Gamma Jet progenitor of the GRB is leaving a trace in the space: usually a nebulae where the nearby ISM may record the jet sweeping as a three dimensional screen. The outcomes maybe either a twin ring as recent SN1987A has shown, or helix traces as the Cat Eye Nebula or more structured shapes as plerions and hourglass nebulae. We imagine the jet as born by a binary system (or by an asymmetric disk accreting interaction) where the compact companion (BH or NS) is the source of the ultra relativistic electron pair jet (at tens GeV. Inverse Compton Scattering on IR thermal photons will produce a collinear gamma jet at MeV). The rarest case where the jet is spinning and nearly isolated would produce a jet train whose trace are star chains as the Herbig Haro ones (Fargion, Salis 1995). When the jet is modified by the magnetic field torque of the binary companion field the result may be a more rich cone shape. If the ecliptic lay on the same plane orthogonal to the jet in an ideal circular orbit than the bending will produce an ideal twin precessing cones which is reflected in an ideal twin rings (Fargion, Salis 1995). If the companion is in eccentric orbit the resultant conical jet will be more deflected at perihelion while remain nearly undeflected at a aphelion. The consequent off-axis cones will play the role of a mild ”rowing” acceleration able to move the system and speed it far from its original birth (explosive) place. Possible traces are the asymmetric external twin rings painted onto the spherical relic shell by SN1987a. Fast relics NS may be speeded by this processes (Fargion, Salis 1995a, 1995b, 1995c). Because of momentum conservation this asymmetric rowing is the source of a motion of the jet relic in the South-East direction. In extreme eccentric system the internal region of the ring are more powered by the nearby encounter leading to the apparent gas arcs. If the system is orbiting in a plane different from the one orthogonal to the jet the outcoming precessing jet may spread into a mobile twin cone whose filling may appear as a full cone or a twin hourglass by a common plerion shape. At late times there is also possible apparent spherical shapes sprayed and structured by a chaotic helix. External ISM distribution may also play a role enhancing some sides or regions of the arcs. The same presence of external relic shells may explain the possible new evidences of iron lines observed by Beppo-Sax satellite. We may explain these rare signal as the re-illuminated isotropic emission due to the under-line beamed inner light-house $`\gamma `$ Jet. The integral jet spray in long times may mimic even spherical envelopes but internal detailed inspection might soon reveal the thinner jet origin (as in recent Eta Carina string jets). Variable nebulae behaviours recently observed are confirming our present scenario.

## 8 Conclusions

Gamma Ray Burst call for extreme isotropic power or extreme collimation or both. Fore a decade the isotropic Fireball story hold the market. Now a Fireball-Jet compromise start to float on the community. Like a ghost.
We strongly believe that GRBs and SGRs are persistent blazing flashes from light-house thin $`\gamma `$ Jet spinning in multi-precessing (binary, precession, nutation) mode. These Jets are originated by NSs or BH in binary system or disk powered by infall matter: the Jet is not a single explosive event even in GRB, but they are powered at maximal output during SN event.The Jet power is comparable to SN at its peak; the $`\gamma `$ Jet has a chain of progenitor identities: it is born in most SN and or BH birth and it is very probably originated by very collimated primary muon pairs at GeVs-TeVs energies. These muons could cross the dense target matter around the SN explosions, nearly transparent to photon-photon opacities. These muon progenitors might be themselves secondary relics of pion decays or even by a more transparent beamed ultra-high energy neutrino Jet originated (by hadronic and pion showering) near the NS or BH. The high energy relativistic muons decay in flight in electron pairs is itself source of GeVs relativistic pairs whose Inverse Compton Scattering with nearby thermal photon is the final source of the observed hard $`X`$ \- $`\gamma `$ Jet. The relativistic morphology of the Jet and its multi-precession is the source of the puzzling complex $`X`$-$`\gamma `$ spectra signature of GRBs and SGRs. Its inner internal Jet contain, following the relativistic Inverse Compton Scattering, hardest and rarest beamed GeVs-MeVs photons (as the rarest and long life EGRET GRB940217 one) but its external Jet cones are dressed by softer and softer photons. This onion like multi Jets is not totally axis symmetric: it doesn’t appear in inner structure, on its front, as a concentric ring serial; while turning and spraying around it is deformed (often) into an elliptical off-axis concentric rings preceded by the internal harder center leading to a common hard to soft GRBs (and SGRs) train signal. In our present model and simulation this internal effect has been here neglected without any major consequence. The complex variability of GRBs and SGRs are simulated successfully by the equations and the consequent geometrical beamed Jet blazing leading also to the observed rare $`X\gamma `$ signatures. As shown in fig 2 the slightly different precessing configurations could easily mimic the wide morphology of GRBs as well as the surprising rare X-ray precursor shown in Fig. 3 above.