15 - Introducing gene models#

In the next part of the module, we will look at models for the regulation of gene expression in cells, and how simple feedback loops underlie some of the basic dynamics that cells exhibit: regulation, switching, and oscillation. Before we do this, it is important to have at least a basic idea of what it is we are trying to model.

A quick guide to cells and genetic networks#

All living things are made of cells - self-contained structures bounded by a cell membrane (and sometimes also a cell wall). Many organisms exist predominantly as single cells (unicellular, e.g. bacteria, amoebae), while others exist as patterned collections of cells (multicellular, e.g. animals, most plants). While different cells possess and maintain a well-defined identity (which is how we can name them), they are far from static structures, and depend on balanced dynamical processes. Mathematical modelling plays an important role in understanding these dynamical processes and how they are regulated in cells.

Unicellular organisms encounter variable environments, and therefore need to be able to adapt their behaviour. For example, a bacterium living in your gut will have to adapt its metabolism to whatever food you present it with. Adaptability often requires a cell to switch between different behaviours (e.g. between metabolising different types of sugar). Multicellular organisms often develop from a single cell through sequential rounds of cell division (the process in which a cell duplicates its contents and then splits to make two cells). However, the cells in multicellular organisms do not remain identical to each other, but take on stable and well-defined characteristics, so we can talk about skin cells, liver cells,muscle cell, blood cells, etc. in a meaningful way. This adoption of distinct well-defined characters is called differentiation, and again requires cells to be able to switch their behaviour.

How do cells manage to change or switch their state in a coherent way? They use a combination of information from their environment and internal mechanisms. Central to the internal mechanisms are what we call genetic regulatory networks (or gene networks). Almost all cells contain large structures centred on DNA (deoxyribonucleic acid). For our purposes, DNA can be viewed as a set of long, linear sequences of letters (A, C, G and T). This sequence is highly stable, is accurately copied during cell division, and is identical in all cells of a multicellular organism. The full DNA sequence in a cell is referred to as the genome.

How is the genome involved in the regulation of the state of a cell (either stably maintaining it or switching it)? A key concept is that of the gene, and the basic unit of dynamics that we will concern ourselves with in this module is the expression of a gene. For our purposes, a gene is a defined subset of the DNA sequence that can code for (i.e. be copied, or transcribed to) an mRNA molecule. In turn the mRNA produced from the gene codes for (i.e. can be translated into) a particular protein molecule. We will refer to the regulated production and degradation of the mRNA and protein corresponding to a gene as the expression of that gene. The amount of the gene itself does not change, and so the dependent variables of our models will be the concentrations of mRNA and protein in a cell.

In the next section, we will consider the regulated expression of a single gene. Later, we will consider how genes can interact by regulating each others’ expression. We will again focus exclusively on differential equation models, in which the concentrations/amounts of mRNA and protein are represented by continuous variables. These are essentially population models, just like the ones in the module so far.

A first model#

In this section, we consider the dynamics of gene expression in two scenarios:

  1. The expression of the gene is regulated by some external transcription factor (whose concentration may or may not be a function of time).

  2. The expression of the gene is regulated by its own protein product (which is therefore a transcription factor).

The models have two dependent variables:

\(M(t)\) : the concentration of mRNA.

\(P(t)\) : the concentration of protein.

The dynamics of \(M(t)\) and \(P(t)\) are regulated by production and degradation. In this context, the molecular processes underlying the production of mRNA and protein are called transcription and translation, respectively (see above).

We will make the following assumptions:

  1. The rates of mRNA and protein degradation are proportional to their concentration (i.e. degradation is linear ).

  2. The rate of translation is proportional to the mRNA concentration (i.e. translation is linear ).

Then the general form of the models, written as a pair of ordinary differential equations, is

(66)#\[\begin{align} \frac{dM}{dt} &= R(t) - \mu M\\ \frac{dP}{dt} &= kM - \nu P. \end{align}\]

where \(\mu\) is the degradation rate of mRNA, \(\nu\) is the degradation rate of protein, \(k\) is the translation rate, and \(R(t)\) is the transcription rate. Note that for the model to make sense biologically, \(\mu\), \(\nu\) and \(k\) must be positive constants, and \(R(t) \geq 0\). Note also that if \(R(t)\) is not constant, then the model will not go to an equilibrium.

Time-courses#

Starting with the dynamics of mRNA expression, we see that it is independent of \(P\), so we can solve it in isolation by using an integrating factor. We have

(67)#\[\begin{equation} \frac{dM}{dt} +\mu M = R(t) \end{equation}\]

Have a go

Using an integrating factor, show that the solution is,

(68)#\[\begin{equation} M(t) = e^{-\mu t}\int e^{\mu t}R(t)dt. \end{equation}\]
Click for solution

The integrating factor will be \( e^{\mu t}\). We now multiply both sides by this integrating factor to find,

\[\begin{align*} & e^{\mu t}\frac{dM}{dt} +e^{\mu t}\mu M = e^{\mu t}R(t)\\ \implies & \frac{d}{dt}e^{\mu t}M=e^{\mu t}R(t), \end{align*}\]

since the left-hand side is the inverse of the product rule for differentiation. The form of the left-hand side is one that we can readily integrate (that is the whole purpose of using an integrating factor after all). However, because \( R(t)\) is some (as yet) unknown function of time we cannot yet compute that integral. This means we can write down the solution as,

\[\begin{align*} &e^{\mu t}M=\int e^{\mu t}R(t)dt\\ \implies &M(t)=e^{-\mu t} \int e^{\mu t}R(t)dt. \end{align*}\]

Note that we didn’t account for the constant of integration when we integrated the left-hand side; we will deal with this if/when we integrate the term on the right-hand side.

It looks like the function \(R(t)\) is going to be very important for how the time course plays out. Let us initially assume that \(R(t)=r\), that is, it is constant. As such we can simplify the equation to,

(69)#\[\begin{equation} M(t) = re^{-\mu t}\int e^{\mu t}dt. \end{equation}\]

Recalling how we integrate exponenitals, we can then find,

(70)#\[\begin{align} M(t) &= \frac{r}{\mu}e^{-\mu t}\left[e^{\mu t}+C\right],\\ &= \frac{r}{\mu}\left[1+Ce^{-\mu t}\right], \end{align}\]

where \(C\) is a constant of integration. Now let us assume we know the initial concentration of mRNA at time \(t=0\) is \(M_0\). We can use this to find our constant of integration as,

(71)#\[\begin{equation} C=\frac{M_0\mu}{r}-1. \end{equation}\]

Substituting back in allows us to reach the final expression of,

(72)#\[\begin{equation} M(t) = \frac{r}{\mu}+\left[M_0 - \frac{r}{\mu}\right]e^{-\mu t}. \end{equation}\]

From this we can see that as \(t\to\infty\) the concentration will tend towards an equilibrium value of \(M=r/\mu\). Also, as the expression shows we have exponential decay or growth, we can get an idea of the speed at which the concentration approaches the equilibrium by finding its half-life. In particular, assume that \(M_0=0\) (i.e. there is no mRNA until time 0 when we “switch the gene on”). This means,

(73)#\[\begin{equation} M(t) = \frac{r}{\mu}(1-e^{-\mu t}). \end{equation}\]

Half the equilibrium value will be \(r/(2\mu)\). If we call the time at which this value is reached \(\tau\), then we have,

(74)#\[\begin{equation} \frac{r}{2\mu} = \frac{r}{\mu}(1-e^{-\mu \tau}). \end{equation}\]

Some inspection finds that this requires \(e^{-\mu\tau}=1/2\), meaning \(\tau=\ln(2)/\mu\). Therefore the faster the degradation rate of the mRNA, the more slowly it reaches its equilibrium.

What about the time-course of \(P(t)\)? Now that we have an expression for \(M(t)\) we can also turn the equation for \(dP(t)/dt\) into a linear ODE, and can again use an integrating factor.

(75)#\[\begin{align} &\frac{dP}{dt}+\nu P = k M(t)\\ \implies &\frac{d}{dt}[Pe^{\nu t}]=kM(t)e^{\nu t}\\ \implies &P(t)=ke^{-\nu t}\int M(t)e^{\nu t}dt. \end{align}\]

Let us again assume \(M(0)=M_0=0\) and also that \(P(0)=P_0=0\). We just found that in this case \(M(t)=r/\mu(1-e^{-\mu t})\). Therefore,

(76)#\[\begin{align} P(t)&=\frac{kr}{\mu}e^{-\nu t}\int (1-e^{-\mu t})e^{\nu t}dt,\\ &=\frac{kr}{\mu}e^{-\nu t}\int (e^{\nu t}-e^{(\nu-\mu) t})dt. \end{align}\]

We need to consider the cases where \(\nu=\mu\) and \(\nu\neq\mu\) separately.

\(\nu=\mu\)

(77)#\[\begin{align} P(t)&=\frac{kr}{\mu}e^{-\mu t}\int (e^{\mu t}-1)dt,\\ &=\frac{kr}{\mu}e^{-\mu t}\left(\frac{1}{\mu}(e^{\mu t}-1)-t+C\right) \end{align}\]

With the initial condition \(P(0)=0\) we can find that \(C=0\), and we can re-arrange slightly to give,

(78)#\[\begin{equation} P(t)=\frac{kr}{\mu^2}\left(1-(1+\mu t)e^{-\mu t}\right) \end{equation}\]

\(\nu\neq\mu\)

(79)#\[\begin{align} P(t)&=\frac{kr}{\mu}e^{-\nu t}\int (e^{\nu t}-e^{(\nu-\mu) t})dt,\\ &=\frac{kr}{\mu}e^{-\nu t}\left(\frac{e^{\nu t}}{\nu}-\frac{e^{(\nu-\mu) t}}{\nu-\mu}+C\right). \end{align}\]

Hopefully you can see from this expression why we needed to take the case \(\nu=\mu\) separately. Recalling we set \(P(0)=0\) we can calculate the constant of integration \(C=\mu/(\nu(\nu-\mu))\).

Have a go

Substitute in the expression for \(C\) and rearrange to find the final solution,

(80)#\[\begin{equation} P(t)=\frac{kr}{\mu\nu}\left(1+\frac{\mu}{\nu-\mu}e^{-nu t}-\frac{\nu}{\nu-\mu}e^{-\mu t}\right). \end{equation}\]

In either case we see that as \(t\to\infty\), \(P\to kr/(\mu\nu)\). The figure below shows the time-course of the system for some example parameter values.

Time-course of single gene model

Figure: Time-course of single gene model, with \(r=1\), \(\mu=0.5\), \(\nu=0.25\), \(k=0.1\) and \(M(0)=P(0)=0\).

So, in this model, even when we started with 0 concentrations of both mRNA and protein, both eventually tend towards equilibrium values, with the speed of that approach largely dictated by the mRNA degradation rate.

3 key points#

  1. We can model the dynamics of cells, genes and proteins in much the same way as we modelled larger populations.

  2. In a simple model of gene expression, the equations are linear and can be solved with integrating factors.

  3. Even when both concentrations start at 0, they will approach a non-zero equilibrium.