*Repressilator*, that consists of a network of three proteins, each of which represses the transcription of the next (akin to the rock-paper-scissors game). Specifically, we examine a simple deterministic model of the dynamics of the Repressilator (phrased as a set of coupled ordinary differential equations), and study some properties of that system. In a subsequent exercise, Stochastic and Deterministic Repressilator, we example a more complex model of the Repressilator, and compare deterministic models of its dynamics with stochastic models.

### Learning Goals

*Science*: You will learn about a simple engineered gene
regulatory network, the Repressilator, and explore a simple model of
its dynamics. You will examine the phase diagram of the Repressilator
for different parameter combinations.

*Computation*: ODE integration and root finding

### Procedure

- Read Elowitz and Leibler, "A synthetic oscillatory network of transcriptional regulators", Nature 403, 335-338, 2000.
- Start from the following equations of motion for the concentrations of
the mRNAs m and the proteins p, and derive the reduced equations found
under the heading "Deterministic, continuous approximation" found on p. 337
of the Repressilator paper:
The scaling is done as described in the paper, or alternatively, as described in Philips Physical Biology of the Cell: time is rescaled in units of the mRNA lifetime; protein concentrations are written in units of K_b**(1/n) (or K_M in the Repressilator paper); and mRNA concentrations are rescaled by their translation efficiency, the average number of proteins produced per mRNA molecule.

- Download the SimpleRepressilator Hints file and rename it as SimpleRepressilator.py.
- Flesh out the definition of the
`HillRepressilator`class. (We call it this since the interaction of each protein with the promoter it represses is approximated with a simple Hill function. This will involve steps described below. - Implement the
`__init__`method for the`HillRepressilator`. This mostly involves just setting input parameter values to instance variable (e.g.,`self.alpha = alpha`). Also set the dynamical state of the system (`self.y`) to the specified initial condition`y0`.- Note: The default parameter values given in the
`__init__`method for`alpha`,`alpha0`, and`beta`are chosen such that the Repressilator should spontaneously oscillate, rather than converge to a steady state. See Fig. 2(b) of the Repressilator paper: the X in parameter space is roughly at`alpha=216`,`beta=5`, with`alpha0/alpha=0.001`, and`n=2`. (Note that beta as defined on the y-axis of Fig. 2(b) is the inverse of that described on p. 337 of the paper, as lifetimes are proportional to the inverse of decay rates.)

- Note: The default parameter values given in the
- Implement the method
`dydt(y,t,alpha,n,alpha0,beta)`that returns the right hand side of the scaled dynamical equations. This will be used with`odeint`for integrating the dynamics in time. You'll notice in the hints file that this is declared as a`staticmethod`, which is a way of defining a method without having the`self`passed as the implicit first parameter to the method. We need to do this since`odeint`needs a right-hand-side function that takes`y0`and`t`as argument, not an instance of a class. See the Python docs to learn more about staticmethods. (Alternatively, if you're feeling very adventurous, define instead a method`make_dydt`that returns a function for use with`odeint`, with the current values of the model parameters`alpha`,`alpha0`, and`beta`hardwired in, rather than passed as`args`.) - Implement the
`run`method. This should integrate the dynamical equations for a specified time interval, and store the trajectory in an instance variable`self.traj`. Since you might want to call`run`multiple times in succession, you will want to append the result of the integration to`self.traj`rather than overwriting it with each function call.- NOTE:
`scipy.append(array1, array2)`can be used to append one array to another, but its behavior is different than the`append`method for Python lists. Whereas the`list.append`method modifies the list in place (and returns`None`), the`scipy.append`method*concatenates*two arrays and returns the result. (I find this choice unfortunate.)

- NOTE:
- Implement a method
`plot`that will plot the dynamical trajectory of the Repressilator. You might want to include options for showing just the mRNA or protein concentrations, as well as an option for converting back from scaled units to absolute units to compare with Fig. 1(c) of the Repressilator paper. - Implement a method
`find_steady_state`to find the steady-state concentrations of mRNAs and proteins, based on the values of the model parameters. This will involve using`scipy.optimize.fsolve`to find a set of concentrations`cstar`for which the concentrations do not change in time. - The steady state is stable only for certain parameters. Where it
is unstable, it gives rise to the oscillations of interest. Fig. 1(b) of
the Repressilator paper shows a phase diagram delineating the stable
and unstable regimes in parameter space. Write a method
`plot_phase_boundary`that sweeps over parameters`alpha`, and generates the corresponding values of`beta`, using`fsolve`, that delineate the stability boundary. Make a plot of the phase boundary for the parameter values specified in Fig. 1(b).- If you're interested in a derivation of this stability result,
look at
section 19.3.3 of
*Physical Biology of the Cell*

- If you're interested in a derivation of this stability result,
look at
section 19.3.3 of
- Implement a method
`in_steady_state`that returns`True`if the system is in a steady state, or`False`if it is not. Compute the instantaneous velocity (i.e., that returned by`dydt`) and see if its vector norm is small (e.g., below 1.0e-3). - For various sets of parameters in the alpha-beta plane, compute whether the Repressilator has a stable steady state. Plot the sets of parameters for which the steady state is stable, and overlay this with a plot of the computed stability boundary.