Gene Expression

Modeling results

The main process we want to model in this section is the expression and translation of the genes involved in our circuits. Qualitatively, this whole process is triggered by the attachment of an enzyme that synthesizes RNA, known as RNA polymerase (RNAP), to the DNA near a gene. Promoters contain specific DNA sequences such as response elements that provide a secure initial binding site for RNA polymerase and for proteins called transcription factors that recruit RNAP. Then, mRNA associated with that gene is generated which can then be translated into a protein by ribosomes. This second process (translation) is triggered when a ribosome binds at a special region, called Ribosomal Binding Site (RBS). At the end of the day, a protein is produced.

At this point, the protein might be in its function conformal structure or not. We will stop at this point for now and model this whole process from rather coarser level of details: we will simulate populations of mRNAs and proteins and predict their concentration over time. As a sneaky-peek of what is to come, this is the basic prediction of our model:

Notice that when Rhamnose concentration (black line) drops, TetR also drops and HokD starts to be produced, killing the cell.

Main objective

Before we start dwelling with our mathematical modeling, let's state what are our main objectives and goals with this model. We want to ...

predict protein/enzyme concentration
take into account the inducer concentration
...

Deterministic modeling

As mentioned above, we want to model the gene expression. Fundamentally, RNA polymerase (RNAP) binds to a transcription factor and triggers the transcription of a gene into mRNA; subsequently, ribosomes starts the translation of that mRNA into a protein/enzyme. The schematics below summarizes it.

Each process is based on the combination of two elements: first, a RNAP and the gene; then, ribosome and mRNA. The simplest way to model such combinations and their products is to use a Dynamical Systems approach [1,2], similar to techniques employed in Population Dynamics [3]. For a introductory text on this subject, click here.

RoxA and Lcp Production

In particular, we assume an inducible promoter such as pLac (BBa_R0010) or pRha (BBa_K914003). Let $P_{p}(t)$ be the concentration of promoter inducer over time, which will control the level of expression of our circuit. The concentration of RoxA and Lcp proteins are $P_{r}(t)$ and $P_{\ell}(t)$, respectively. Finally, TetR concentration is denoted by $P_{t}(t)$, which blocks the Kill Switch mechanism, hindering cell growth by the production of HokD with concentration $P_{h}(t)$.

Let's start defining the differential equation for the production of RoxA and Lcp. Their concentration depend only on the concentration of promoter inducer and their own degradation rate. Let's define $\delta_r$ and $\delta_{\ell}$ as the degradation rate of RoxA and Lcp respectively. Additionally, to set the promotion levels, we define the $\beta$ parameter. Then, the differential equations become:

\[ \frac{dP_{r}(t)}{dt} = - \delta_r P_{r}(t) + \beta \frac{P^n_{p}(t)}{K^n + P^n_{p}(t)} \] \[ \frac{dP_{\ell}(t)}{dt} = - \delta_{\ell} P_{\ell}(t) + \beta \frac{P^n_{p}(t)}{K^n + P^n_{p}(t)} \]

Notice that we have modeled the promoter activity as a Hill function. Jointly, tetR will be also expressed and its differential equation can be written as follows:

\[ \frac{dP_{t}(t)}{dt} = - \delta_{t} P_{t}(t) + \beta \frac{P^n_{p}(t)}{K^n + P^n_{p}(t)} \]

Kill Switch: hokD

The Kill Switch mechanism induces a self-destruction mechanism, based on a toxic enzyme, HokD. It is regulated by pTet (BBa_R0040). This assures cells won't leaeve. expression will then depend on $P_t (t)$,

\[ \frac{dP_{h}(t)}{dt} = - \delta_{h} P_{h}(t) + \beta_0 + \beta_t \frac{P^n_{t}(t)}{K^n + P^n_{t}(t)} , \]

where $\beta_0$ is the leakeage term [2] and $\beta_t$ defines the level of promotion of the tetR promoter.

Fitting the parameters

Precisa passar do caderno

Fluctuation analysis

Deterministic models within a context of population dynamics (i.e., populations of chemicals evolving and interacting in time) are well suited for very large concentrations, where fluctuations are indeed suppressed and a clear average effect tends to dominate. If that is not the case, fluctuations are rather important and fundamental to describe the phenomenon. To describe some level of randomness in the system (such as the possibility of RNAp not binding to a promoter) we use Stochastic models. In general, chemical reactions are somewhat at the edge of both cases: fluctuations are expected and may give rise to interesting effects (such as stochastic resonance), but averages is a very good first step to describe our quantities of interest. However, there is a price to using Stochastic Processes for modeling: they are considerably more expensive and complex than their deterministic counterparts. Thus, beyond doubt, it would be essential for us to have some kind of control or estimation on fluctuations over the deterministic approach presented above. To achieve this, we will next introduce a classical technique in Stochastic Processes and show some analytical and numerical procedures to embed and investigate fluctuations in our modeling.

Interestingly enough, there is a way to closely related Stochastic and Deterministic models: it is not very hard to construct a Stochastic Process whose average dynamics matches quite well with a given Deterministic model. To better understand this bridge, suppose you have a stochastic process describing the evolution of a variable $X(t)$. You want the average evolution of $X(t)$ to match the dynamics a given Deterministic model. You can calculate the $n$-th order momentum $\langle X^n \rangle (t) = E\big[X(t)^n\big]$, where $E[\cdot]$ denotes a probability average (or expectation), as function of time. For instance, supposing you have a master equation (as we will below), then you can manipulate it into a differential equation for $\langle X \rangle (t)$ (i.e., average of $X$) that depends on $\langle X^2 \rangle (t)$. If you write another differential equation, from you master equation, for $\langle X^2 \rangle (t)$ it will depend on $\langle X^3 \rangle (t)$. Guess what? The equations will never close, always depending on higher and higher order momenta. This is a problem, since you would have an infinite number of equations to solve!!

The simples possible way to solve this issue is as follows. Thus, as a first approximation to the equation for $\langle X \rangle (t)$ you can neglect $\langle X^2 \rangle$ term as simply $\langle X \rangle^2$, having then a equation that depends only on the first order momentum $\langle X \rangle$. As a second approximation, you can go up to the equation for $\langle X^2 \rangle (t)$ and regard $\langle X^3 \rangle (t)$ as $\langle X \rangle^3 (t)$. Then you will end up with two equations that depend only on $\langle X \rangle (t)$ and $\langle X^2 \rangle$. Techniques based on this principle are called momentum-closure techniques.

Back to top

Team:Brasil-USP/Modeling/GeneExpression