Team:EPF Lausanne/Modeling
Kinetic model
In our project, multiple transistor elements are assembled to create logic gates. We envision the chaining of such gates in order to create complex logic circuits within cells. Predicting the behaviour of these complex cascades of reactions and way they reach a stationary state can be challenging. Modeling here represents the best way to understand our system and quantify the influence of various biochemical parameters. On the long run, modeling may also help fine-tune the design of biological systems and wet lab experiments.
We attempted to model the dynamics and interactions our system's components to predict the temporal response of our biological circuits. In this kinetic model, time dependency of the concentration of different species is taken into account explicitly.
Assumptions
Modelers of biophysical system often have a difficult mission. They must strike a compromise between developing a very detailed and accurate model and producing a simplistic one that may be easily analysed and compared to experimental data.
Developing model implies making a certain number of simplifications and assumptions to enable it to conform to reality. In this section we will keep track of the assumptions we made in order to clarify our model. Since the same system may be modeled using different approaches, we will justify our choices and set of assumptions.
The most important assumption underlying our kinetic model is the fact that the concentration of a given species does not depends on spatial coordinates, i.e. it is the same in every region of the cell. This assumption is quite bold but is a common approach in the literature and usually produces accurate results. The validity of this assumption is limited by the small number of molecules present in the cell and localization in membranes or particular organels. In some cases it is necessary, it may be necessary to consider the stochastic behavior of individual trajectories rather than global averages [1].
One of the challenges we faced while building the model was the lack of readily available kinetic constants. Since dCas9 is a newly discovered gene regulation technology [2], gRNA/dCas9 and gRNA+dCas9/DNA binding/unbinding kinetics are still under investigation. One may imagine that the binding/unbinding constants are somewhat related to the gRNA sequence given the different chemical properties of nucleotide sequences. However, we will consider binding/unbinding constants to be gRNA-independent. For the gRNA/dCas9 interaction, this assumption is justified by the fact that the gRNA scaffold is always the same. For the gRNA+dCas9/DNA interaction, we can justify this hypothesis by thinking of an average nucleotide composition, as our synthetic sequences are randomly generated
dCas9 degradation is a fundamental process in our system. Since high levels of dCas9 are toxic for the cell [3], a constantly varying dCas9 population is required for a signal to propagate from one gate to another (remember, the output of a gate is a gRNA which will bind to a free dCas9 in order to propagate the signal to the next gate, cf Project Description). We can imagine three different types of degradation for the gRNA/dCas9 complex : degradation of the whole complex, dCas9 detaching from its gRNA and the degradation of the gRNA's targeting sequence (resulting in a non-functional but occupied dCas9). Since the probability of a gRNA unbinding from a dCas9 protein is extremely low [4], we only considered the degradation of the whole complex, thus neglecting the unbinding reaction. In addition, we assumed that this degradation rate is the same of gRNA-free dCas9.
In our project, the dCas9/gRNA complex is used as a transcription factor. It may alternatively activate or inhibit a targeted promoter. An important assumption of our model is that the number of transistors (i.e. the number of promoters) does not change. This implies that the total number of transistors is constant. We defined six different states for a transistor: activated (\(a\)), inhibited (\(i\)), basal (\(b\)), activated/inhibited (\(ai\)), double inhibited (\(ii\)) and activated/double inhibited (\(aii\)). Mathematically, this translates into: \[ [Ta] + [Tb] + [Ti] + [Tai] + [Tii] + [Taii] = \text{cst.} \] where the simultaneous existence of all states was not tested experimentally. This reduces the degrees of freedom of our system of ODEs: \[ \dfrac{d[Ta]}{dt} + \dfrac{d[Tb]}{dt} + \dfrac{d[Ti]}{dt} + \dfrac{d[Tai]}{dt} + \dfrac{d[Tii]}{dt} + \dfrac{d[Taii]}{dt} = 0. \]
To have only first-order reactions we assume that we can reach transistors states with two dCas9/gRNA complexes bound only through simple states. THis imply that is impossible to go from \(Tai\) to \(Tb\) in a single time step: fist a dCas9/gRNA complex detach (leaving the transistor in an active or inhibited state, depending on which complex detached), then the remaining complex can leave as well. This assumption relies on the fact that we neglect the possible interaction of two dCas9/gRNA complexes when bound to the same promotor on different position; in this case the simultaneous detachment of the two dCas9/gRNA complexes has a low probability. The same apply for the states \(Tii\) and \(Taii\).
Summary
Equations
Expliciting the large set of ordinary differential equations (ODEs) governing our system is a nontrivial and error prone process. When few species are present, doing it manually and double-checking the equations can be sufficient. However, when the number of transistors composing our system increases (due to the chaining of gates), keeping track of all gRNAs and their activating/inhibiting effect on specific promoters (when bound to dCas9) become almost impossible. For this reason we created a Python program which does the tedious tasks for us (see the Software section). With this script, inputting the circuit structure (i.e the gates with input/outputs) will generate of the entire kinetic model, write down the equations in LaTeX format and create a Python function representing the apposite system of ODEs. This system may then be solved numerically.
Activation
Our activation model consists of a simple transistor that takes a gRNA/dCas9 complex as an input (A) and produces a gRNA (C). The gRNA/dCas9 complex enhances the production of the output C, which is otherwise produced basally. In our experimental setup, the output C is used to enhance the production of GFP to enable a quantitative measurement of the activation (with respect to the basal expression).
dCas9 is constitutively produced from a low copy plasmid and its degradation is proportional to its concentration. The population of unbound dCas9 is also affected by the binding between dCas9 and a gRNA. The unbinding of this two species is neglected (see justification above). \[ \frac{d}{dt}[\mathit{dCas9}] = K_{dCas9} -\Gamma_{dCas9}[\mathit{dCas9}] -R_{dCas9}[\mathit{dCas9}][\mathit{gRNA}] \]
gRNA is produced by an inducible promoter. gRNAs are degraded proportionally to their concentration. Unbound gRNAs may also bind to a free dCas9 protein. As stated previously, the dCas9/gRNA complex is degraded as a whole (at the same rate as unbound dCas9 degradation). Therefore, there is no gRNA production that results from the complex's dissociation. \[ \frac{d}{dt}[\mathit{gRNA}] = K_{IP}[\mathit{IP}] -\Gamma_{gRNA}[\mathit{gRNA}] - R_{dCas9}[\mathit{dCas9}][\mathit{gRNA}] \]
dCas9/gRNA complexes are our gene regulatory units: depending on the targeted site of a promoter, these complex can enhance or inhibit gene transcription. To test activation, gRNA sequences were generated to exclusively target activating sites. In this case, the binding of dCas9/gRNA to the promotor enhances transcription, which is otherwise in a basal expression. The binding of a dCas9/gRNA to a promoter creates an activated transistor \(Ta\), which is otherwise in a basal state \(Tb\). \[ \frac{d}{dt}[\mathit{dCas9w\text{-}gRNA}] = R_{dCas9}[\mathit{dCas9w}][\mathit{gRNA}] +D_{Ta}[\mathit{Ta}] -\Gamma_{dCas9} [\mathit{dCas9\text{-}gRNA}] -R_{Ta}[\mathit{Tb}][\mathit{dCas9\text{-}gRNA}] \]
Our simple transistor (used to simulate activation) can be found in two states: activated (\(Ta\)) or basal (\(Tb\)). The switch between basal and activated states is obtained by the binding/unbinding of the dCas9/gRNA complex: \[ \frac{d}{dt}[\mathit{Tb}] = D_{Ta}[\mathit{Ta}] -R_{Ta}[\mathit{Tb}][\mathit{dCas9\text{-}gRNA}] \] \[ \frac{d}{dt}[\mathit{Ta}] = R_{Ta}[\mathit{Tb}][\mathit{dCas9\text{-}gRNA}] -D_{Ta}[\mathit{Ta}] \]
In order to assess the functionality of our transistor and the targeting of the dCas9/gRNA complex, we used a measurable output: green fluorescent protein (GFP): \[ \frac{d}{dt}[\mathit{GFP}] = L_{b}[\mathit{Tb}] + K_{GFP}[\mathit{Ta}] -\Gamma_{GFP}[\mathit{GFP}] \] The model also considers the leakiness of GFP, i.e. the production of GFP when the transistor is found in a basal state.
Inhibition
Inhibition works similarly to activation, with the difference of the targeting site of the dCas9/gRNA complex. In the case of activation, the dCas9/gRNA complex targets a sequence on DNA such that the RNA polymerase (RNAP) recruiting unit (\(\omega\) unit in E. Coli and VP64 in yeast) is close to the RNAP binding site: the recruiting unit enhance RNAP binding and therefore gene transcription. In the case of inhibition, the dCas9/gRNA targeting sequence on DNA is so close to the RNAP binding site that the RNAP is not able to bind anymore; this steric inhibition prevents gene transcription.
Differential equation describing the activation remains substantially unchanged. The only change is the possible states of our transistor, which can be in its basal state (nothing bound to it) or inhibited (dCas9/gRNA complex bound close to the RNAP binding site). Thus the equation representing the variation of \([Ta]\) is now substituted by the following equation: \[ \frac{d}{dt}[\mathit{Ti}] = R_{Tb}[\mathit{Tb}][\mathit{dCas9\text{-}gRNA}] -D_{Ti}[\mathit{Ti}] \]
The production of GFP, which is final measurable output, is proportional to the concentration of inhibited transistors (instead of the activated one): \[ \frac{d}{dt}[\mathit{GFP}] = L_{b}[\mathit{Tb}] + L_{i}[\mathit{Ti}] -\Gamma_{GFP}[\mathit{GFP}] \] The is caused by the fact that a biological transistor is by no mean perfect and we can always have some kind of leakiness. Despite this, we expect to fulfill the condition \[ L_{i} < L_{b} \] and therefore the GFP expression will be lowered in the inhibited case with respect to the basal expression because the presence of inhibited transistors diminish the population of transistors in the basal state (remember, we assumed a constant total numbe of transistors).
Activation and inhibition
The possibility of activating and inhibiting a transistor at the same time, i.e. have dCas9/gRNA complexes bound to activating and inhibiting sites at the same time, is a key point of our project. In fact the design of complex gates relies strongly on the idea that when dCas9/gRNA complexes bind to the same transistor in inhibiting and activating positions at the same time the transistor is inhibited (see LINK for details on complex gates design). This is true in yeast [] but to our knowledge has not been shown in E. Coli yet (see LINK to our experimental results).
In the case of activation and inhibition we have to combine the previous equations. In this case we have a gRNA (A) targeting the activation site, while another gRNA (B) targets the inhibition site. \[ \frac{d}{dt}[\mathit{gRNA}_A] = K_{IP}[\mathit{IP}] -\Gamma_{gRNA}[\mathit{gRNA}_A] - R_{dCas9}[\mathit{dCas9}][\mathit{gRNA}_A] \] \[ \frac{d}{dt}[\mathit{gRNA}_B] = K_{IP}[\mathit{IP}] -\Gamma_{gRNA}[\mathit{gRNA}_B] - R_{dCas9}[\mathit{dCas9}][\mathit{gRNA}_B] \] The presence of the two gRNA modifies the kinetic of dCas9, which can now bind two different species: \[ \frac{d}{dt}[\mathit{dCas9}] = K_{dCas9} -\Gamma_{dCas9}[\mathit{dCas9}] -R_{dCas9}[\mathit{dCas9}][\mathit{gRNA}_A]-R_{dCas9}[\mathit{dCas9}][\mathit{gRNA}_B] \] \[ \frac{d}{dt}[\mathit{dCas9w\text{-}gRNA}_A] = R_{dCas9}[\mathit{dCas9w}][\mathit{gRNA}_A] +D_{Ta}[\mathit{Ta}] -\Gamma_{dCas9} [\mathit{dCas9\text{-}gRNA}_A] -R_{Ta}[\mathit{Tb}][\mathit{dCas9\text{-}gRNA}_A] \] \[ \frac{d}{dt}[\mathit{dCas9w\text{-}gRNA}_B] = R_{dCas9}[\mathit{dCas9w}][\mathit{gRNA}_A] +D_{Ti}[\mathit{Ti}] -\Gamma_{dCas9} [\mathit{dCas9\text{-}gRNA}_B] -R_{Ta}[\mathit{Tb}][\mathit{dCas9\text{-}gRNA}_B] \] Note that because of the assumption that the recruitment of gRNA from the dCas9 complex is gRNA-independent (and the same is assumed true for the dCas9/gRNA complex binding to the promotor) we can obtain a system of ODEs similar to simple activation or simple inhibition for the total concentration of gRNAs \[ [gRNA] = [gRNA_A] + [gRNA_B] \]
Because we cannot assure that each transistor is in the \(Tai\) state and because we assumed that the activated/inhibited states cannot be reached in a single step from the basal state, we have four possible states for our transistor: basal, activated, inhibited, activated/inhibited. This leads to the following equations:
GFP production ...
Double inhibition
Activation and double inhibition
\[ f(a) = \frac{1}{2\pi i}\oint_\gamma \frac{f(z)}{z-a}dz \]
Constants
| Name | Description | Value | Source |
|---|---|---|---|
| \(K_{dCas9}\) | dCas9 production | 50 | Freiburg |
| \(\Gamma_{dCas9}\) | dCas9 degradation | 50 | Freiburg |
| \(R_{dCas9}\) | gRNA recruitment from dCas9 | 50 | Freiburg |
| \(\Gamma_{gRNA}\) (E. Coli) | gRNA degradation | 0.2 (1/min) | Bernstein et al. [5] |
| \(\Gamma_{gRNA}\) (Yeast) | gRNA degradation | 0.03 (1/min) | Wang et al. [6] |
| \(\Gamma_{dCas9-gRNA}\) | Complex degradation | Freiburg | |
| \(R_{Ta}\) | dCas9 binding to promotor | 50 | Freiburg |
| \(D_{Ta}\) | dCas9 detachment from promotor | 50 | Freiburg |
gRNA production/degradation rates
For the gRNA degradation rate \(\Gamma_{gRNA}\), we considered the mean mRNA half-life of Refs. [5-6] and subsequently computed the degradation rate. In Ref. [6] the mean half-life is explicitly stated, while we computed ourself the mean half-life of Ref. [5] for different E. Coli strains. Note that for an exponential decay, the half-life \(t_{1/2}\) is linked to the lifetime \(\tau\) by \[ t_{1/2} = \tau \log(2) \] and the lifetime \(\tau\) is in turn linked to the decay constant \(\Gamma\) by \[ \Gamma = \frac{1}{\tau}. \]
Simulation
References
[1] R. Phillips et al., Physical Biology of the Cell, Second Edition, Garland Science, 2013.
[2] L. S. Qi et al., Repurposing CRISPR as an RNA-Guided Platform for Sequence-Specific Control of Gene Expression, Cell 152, 1173–1183, 2013.
[3]
[4] Team Duke iGEM 2014, Wiki
[5] Bernstein et al., Global analysis of Escherichia coli RNA degradosome function using DNA microarrays, PNAS, vol. 101, 2758 –2763, 2004.
[6] Wang et al., Precision and functional specificity in mRNA decay, PNAS, vol. 99, 5860 –5865, 2002.
