Team:Oxford/Test/Modeling/Tutorial
Modelling Tutorials
Introduction
On this page we aim to explain the concepts and methodology behind the modelling we have done in our project. We hope this will make our modelling section more accessible for all members of the iGEM community.
Gene Expression Networks
Mass Action
Our goal is to model a system of reactions such as the reversible reaction
\[A \mathrel{\mathop{\rightleftharpoons}^{\mathrm{k_{+}}}_{\mathrm{k_{-}}}} B\]Where \(A\) and \(B\) are our two species and \(k_{+}\) and \(k_{-}\) will be determined shortly. To do this, we make the intuitive assumption that the rate at which something reacts is proportional to how much of it we've got. That is, the rate of change of \(B\) with time is
\[\dfrac{dB}{dt}=k_{+}A\]where \(k_{+}\) is simply a proportionality constant. This assumption is called the law of mass action. This of course is not the complete picture, as we forgot to include the fact that \(B\) can react back to form \(A\). So the complete system of equations to solve is
\[\dfrac{dA}{dt}=k_{-}B-k_{+}A\] \[\dfrac{dB}{dt}=k_{+}A-k_{-}B\]At this point, you can type these into your favourite computing program and solve them.
Michaelis-Menten Kinetics
It is possible to formulate a rate law that describes enzyme-catalysed reactions, this is known as Michaelis-Menten kinetics. In order to derive this rate law we will look at a generalised system of chemical events detailed below:
\[S + E\mathop{\rightleftharpoons}C_{1}\mathop{\rightleftharpoons}C_{2}\mathop{\rightleftharpoons}P + E\]In this system \(S\) represents our substrate, \(E\) our enzyme, \(C_{1}\) our enzyme-substrate complex, \(C_{2}\) our enzyme-product complex, and \(P\) our product.
Initially we will make two simplifications. Firstly, we will combine \(C_{1}\) and C2 into a single complex, \(C\), as we will assume that the time-scale of the conversion \(C_{1}\mathop{\rightleftharpoons}C_{2}\) is much faster than that of the association and dissociation events. Secondly, we will assume that the product never binds with the free enzyme. These two assumptions lead to the simplified network:
\[S + E\mathrel{\mathop{\rightleftharpoons}^{k_{1}}_{k_{-1}}}C\mathop{\rightarrow}^{k_{2}}P + E\]Using the laws of mass action (detailed above) we arrive at the following differential equation model:
\[\dfrac{d}{dt}s(t)=-k_{1}s(t)e(t) + k_{-1}c(t)\] \[\dfrac{d}{dt}e(t)=k_{-1}c(t) - k_{1}s(t)e(t) + k_{2}c(t)\] \[\dfrac{d}{dt}c(t)=-k_{-1}c(t) + k_{1}s(t)e(t) - k_{2}c(t)\] \[\dfrac{d}{dt}p(t) = k_{2}c(t)\]Concentrations are denoted as the lowercase letter eg the concentration of \(S\) is given by \(s\).
You may have spotted that the enzyme is not consumed in this series of reactions. Therefore the total amount of enzyme remains constant. We reflect this in the expression \(e_{T}=e+c\). We can now use this expression to eliminate \(e(t)\) from our model, leaving:
\[\dfrac{d}{dt}s(t)=-k_{1}s(t)(e_{T}-c(t)) + k_{-1}c(t)\] \[\dfrac{d}{dt}c(t)=-k_{-1}c(t) + k_{1}s(t)(e_{T}-c(t)) - k_{2}c(t)\] \[\dfrac{d}{dt}p(t) = k_{2}c(t)\]We have one further simplification to make. This is the rapid equilibrium approximation, by which we assume that equilibrium between \(s+e\) and \(c\) on a much faster time scale than the reaction of \(c\) to \(p\). With this approximation we can write the following equation:
\[0=-k_{-1}c(t) + k_{1}s(t)(e_{T}-c(t)) - k_{2}c(t)\]Leading to:
\[c^{*}(t)= \dfrac{k_{1}s(t)e_{T}}{k_{-1} + k_{2} + k_{1}s(t)}\]Now we can see that \(c\) is no longer an independent variable, but instead tracks the other variables. This can be referred to as \(c\) being in quasi-steady state. By including this new expression into our model we are left with:
\[\dfrac{d}{dt}s(t) = -\dfrac{k_{2}k_{1}s(t)e_{T}}{k_{-1} + k_{2} + k_{1}s(t)}\] \[\dfrac{d}{dt}p(t) = \dfrac{k_{2}k_{1}s(t)e_{T}}{k_{-1} + k_{2} + k_{1}s(t)}\]With this we have an expression describing our enzyme catalysed system in the form of a single reaction. The rate of this reaction is known as a Michaelis-Menten rate law. By defining \(K_{max}=k_{2}e_{T}\) and \(K_{half}=\dfrac{k_{-1}+k_{2}}{k_{1}}\) we can express in the more familiar form:
\[rate \: of \: S \mathop{\rightarrow} P = K_{max}\dfrac{s}{K_{half} + s}\]It is worth noting that there are a number of ways to derive this form, each involving different approximations. These methods may lead to the constants \(K_{max}\) and \(K_{half}\) being defined differently.
Stochastic
The simplest way to add an element of 'randomness' into an equation is to add a random term into the differential equation. This turns it into a stochastic differential equation, such as:
\[\dfrac{dx}{dt}=c+r\]
which determines the evolution of a particle travelling at a speed \(c\), with a random number \(r\) added into the mix just for fun. We choose \(r\) to be normally distributed with a mean value of 0 and a variance of, say, 1. We evaluate the position of the particle at multiple timesteps (with some small time \(dt\) separating them) and plot our solution in the Figure below.
We can do better than this when dealing with systems such as our example in the previous section. We use the Gillespie algorithm to find:
- The time taken between reactions
- The reaction that took place in that time.
as neither of these could be determined by our initial method [2, 4].
1) To find the time taken, we make an assumption about the form of the probability distribution that the time between reactions holds. It makes sense to think of it as a falling exponential, such that reactions more often happen quickly rather than slowly, and that the more probable a reaction is, the shorter the time between subsequent reactions. This naturally leads to the form
\[P(t)=e^{-\alpha_{0}t}\]where \(P(t)\) is the probability of any reaction occuring in time \(t\), and \(\alpha_{0}\) tells us something about the likelihood of the reaction. We plot this in the Figure. In fact we call \(\alpha_{0}\) the sum of the propensities \(\alpha_{i}\) (where the dummy variable \(i\) runs from 1 to the number of possible reactions, \(N\), that could occur). In the case of our example system, the propensity for the reaction \(A\) to \(B\) would be \(k_{+}A\) and the propensity for the reaction \(B\) to \(A\) is \(k_{-}B\) so in this case \(\alpha_{0}=k_{+}A+k_{-}B\).
To decide the time it took for a reaction to occur, we sample a random number \(r_{1}\) in the range [0,1] and stick that number on the y axis of the graph. We read off from our graph a corresponding time and we're done - we now know when the next reaction took place.
2) To decide which reaction took place, we need the probabilities of each reaction occuring being proportional to the ratio of their propensities. We split up the domain of 0 to 1 into the normalised ratios of the propensities of our problem, recalling the propensities we calculated. We could have, in this case, selected the backwards reaction B to A. We can generalise this to the form given in the Gillespie Algorithm above.
