Kinetic model

In our project different transistor elements are put together in order to create logic gates and the idea is to chain these gates in order to create complex logic circuits within cells. Because of this chainability of different elements, we have a cascade of reactions which will eventually reach a stationary state. In order to study in depth the behavior of our system we decided to use a kinetic model, where time dependency of the concentration of different species is taken into account explicitly.

Assumptions

A model which is too complex and too detailed can surely describe experimental results with precision, but it is often difficult to interpret and to understand in details. It is therefore useful to come out with the simplest model that can reproduce the system behavior, which can then be understood in details and used to predict new results.

In order to find a simple model which can reproduce experimental results we have to do many simplifications and assumptions over reality. In this section we will try to keep track of the assumptions we made in order to clarify our model. Every assumption is analyzed in details: assumptions are generally not uniques, but we will try to justify our choices.

The most important assumption underlying the kinetic model is the fact that the concentration of a given specie does not depends on spatial coordinates, i.e. it is the same in every region of the cell. This assumption is quite strong but it's a common approach in the literature and usually gives good results. Note however that within a cell the validity of concentration idea itself can be doubtful, since the number of molecules can be small and since these molecules can be localized to membranes or particular organels; in these cases it is necessary to consider the stochastic behavior of individual trajectories rather than global averages [1].

A challenging problem we faced to build our model concerns kinetic constants. Since dCas9 is a newly discovered gene regulation technology [2], gRNA/dCas9 and gRNA+dCas9/DNA binding/unbinding kinetic is still under investigation. It is easy to imagine that binding/unbinding constants depend (at least sligthly) on the gRNA sequence because of the different chemical properties of nucleotides. However, we will consider that these binding/unbinding constants are 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 generated at random.

dCas9 degradation is a fundamental process in our system: high levels of dCas9 within the cell are toxic [3], thus a continuous change in dCas9 population is needed to propagate the signal from one gate to another (remember, the output of a gate is a gRNA which will bound to a free dCas9 in order to propagate the signal to the next gate). We can imagine three different ways of gRNA/dCas9 complex degradation: the degradation of the whole complex, degradation of the dCas9 leaving the gRNA and the degradation of the targeting sequence of the gRNA leaving a non-functional but occupied dCas9. Since the unbinding probability of gRNAs from dCas9 proteins is extremely low [4], we consider only the degradation of the whole complex neglecting the ubinding; in addition we assume that this degradation rate is the same of the dCas9 not bound to a gRNA.

Our dCas9/gRNA complex is used as a transcription factor, in order to activate or inhibit a promotor. An important assumption of our model is that the number of transistors (i.e. the number of promotors) does not change. This means that the total number of transistors is constant. We have six possible states for the transistor: activated (\(a\)), inhibited (\(i\)), basal (\(b\)), activated/inhibited (\(ai\)), double inhibited (\(ii\)) and activated/double inhibited (\(aii\)). This traduces mathematically to \[ [Ta] + [Tb] + [Ti] + [Tai] + [Tii] + [Taii] = \text{cst.} \] where for a given experimental setup not all possible states are available. 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. \]

Summary

• Concentrations depends only on time
• dCas9 binding/unbinging is gRNA-independent
• gRNAs does not unbind from dCas9
• gRNA/dCas9 is degraded as a complex, whit the same rate as dCas9 alone
• The total namber of transistors is constant

Equations

Writing down the large set of ordinary differential equations (ODEs) governing our system is a nontrivial and error prone process. When the number of species present is low, doing it by hand and double check the equations can be sufficient. However, when the number of transistors composing our system increases (many chained gates), keeping track of all gRNAs and their promoting/inhibiting interaction with DNA (when bound to dCas9) become almost impossible. For this reason we created a Python program which does the dirty job for us (see the Software section): it is sufficient to input the circuit structure (gates with input/outputs) and the program automatically generate the kinetic model, writing it in LaTeX format and creating a Python function representing our system of ODEs, which can be solved numerically in a second time.

Activation

Our activation model consists in a simple transistor taking a gRNA/dCas9 complex as an input (A) and producing an output gRNA (C). The gRNA/dCas9 complex enhances the production of the output C, which is otherwise produced basally. Experimentally, the output C is used to enhance the production of GFP, in order to quantitatively measure the activation (with respect to the basal expression).

dCas9 is constitutively produced from a low copy plasmid and it's degraded proportionally to its concentration. Other changes in free dCas9 population are given by the binding between dCas9 and a gRNA; the unbinding of this two species is neglected. \[ \frac{d}{dt}[\mathit{dCas9}] = K_{dCas9} -\Gamma_{dCas9}[\mathit{dCas9}] -R_{dCas9}[\mathit{dCas9}][\mathit{gRNA_{A}}] \]

The gRNA is produced by an inducible promoter. gRNAs are degraded proportionally to their concentration; another way of diminishing their population is the binding with a free dCas9 protein. As stated previously, the dCas9/gRNA complex is degraded as a whole (at the same rate of dCas9 degradation) and therefore there is no gRNA production from complex 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 target site of a promoter, these complex can enhance or inhibit gene transcription. To test activation, gRNA sequences target only activating sites: the binding of dCas9/gRNA to the promotor enhance 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_{A}}] -R_{Ta}[\mathit{Tb}][\mathit{dCas9\text{-}gRNA}] \]

Our simple transistor needed to test 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 asses the functionality of our transistor and the targeting of the dCas9/gRNA complex we need a measurable output, which in our case is a green fluorescent protein (GFP): \[ \frac{d}{dt}[\mathit{GFP}] = L_{GFP}[\mathit{Tb}] + K_{GFP}[\mathit{Ta}] -\Gamma_{GFP}[\mathit{GFP}] \] The model consider also the leakiness of GFP, i.e. the production of GFP when the transistor is found in a basal state.

Inhibition

Activation and inhibition

Double inhibition

Activation and double inhibition

\[ f(a) = \frac{1}{2\pi i}\oint_\gamma \frac{f(z)}{z-a}dz \]

Constants

gRNA production/degradation rates

For the gRNA degradation rate \(\Gamma_{gRNA}\) we considered the mean mRNA half-life of Refs. [5-6] and we computed the degradation rate subsequently. 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. Remember 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.