Template:CLSB-UK Model MAK

Mass action kinetics


Mass action kinetics is the assumption that the rate of reaction is proportional to concentration of the reactants. We chose to create a mass action kinetics model as we could use relatively few parameters. Reducing the number of parameters, whilst keeping the model representative of the system, made it easier to understand the effect of varying each parameter. We built the model in MATLAB first, and later ported it to JavaScript to make it a better design tool for the biologists and better communicate how the model works.

Graph interactivity

We spent a long time reviewing previous teams models. A problem we found was that most models weren’t communicated in a way that was engaging and easy to understand. We’ve tried to explain our model in a way that someone with no experience in modeling could understand.

We’ve used interactive graphs on this page - when you get to them you can drag the sliders to change parameters and see how they affect the system.

We initially ported our model to a webpage so the design team could try out new stuff when they were designing our toehold switches. This turned out to be very useful for both of us. We decided to make it look a lot nicer and improve its compatibility with different browsers so we could put it on the wiki.

As a webpage it’s much more accessible - anyone (even on a phone) could use it without having to install any custom programs which would take both time to install and might require a license.

As far as we’re aware, we’re the first iGEM team to ever put up something like this. The closest was Pretoria UP 2016, but our graphs also contained our entire mass action kinetics model, instead of just a visualization.

It was built in JavaScript and therefore not limited by proprietary license agreements, and doesn’t require Java applets to be enabled (which modern browsers no longer support due to stability and security issues).[1]

JavaScript follows ECMA’s open standard so it can run on purely open-source tools. This helps acheive “iGEM’s goal of making everything in the competition open source”.[2] There are also a lot more open source libraries available for JavaScript which means other teams can build upon our work easier. These include Chart.js, the library we are using to draw the graphs.

Finally, our JavaScript code ran faster than our MATLAB code, which we attribute to both our greater experience in JavaScript and the aggressive optimisation that’s been done by vendors trying to compete for the fastest browsers.[3]

Initial mass action kinetics model

We created our first model following a meeting with Thomas Ouldridge at Imperial, who recently coauthored a paper about mathematically modeling toehold-mediated strand displacement. In the meeting we agreed that we’d need to model the different states of the toehold switch. We also discussed modeling the miRNA binding and unbinding to and from the toehold base-by-base. We later decided against this, as whilst it may have been useful for learning about the kinetics of toehold switches it would not have given us useful information to put into our design process and would have overcomplicated our model.

We first modelled a simple system of 4 substances with reversible reactions as follows:



Where CTS, PBTS and OTS represent the toehold switch in closed, partially bound and open states respectively. At this point we didn’t know the value of the rate constants so we used placeholders.

This was useful as it gave us experience programming in MATLAB, which we hadn’t done before, as well as forced us to learn the foundations for the maths we’d need later; mainly ODEs. Khan Academy was very helpful with this.

With all three reactions being reversible, MATLABs dsolve function was not able to solve our differential equations. We tried them by hand but were not able to solve them analytically. We had to solve numerically and chose the Euler method. This involves iterating with short time intervals (here we used dt ≈ 1s) between each calculation using a loop to keep a record of how the substances changed over time. We then graphed these values as follows:




We’ve recreated this model so it’s actually running in the browser - try dragging the sliders above the graph to modify the rate constants.

Mass action kinetics model with transcription and translation

The first simple model neglected transcription and translation completely, which meant it was an excessive simplification and consequently was not that useful. To address these shortcomings, we created a second model, derived from parts of Tobias Stögbauer’s PhD, with many more species.[4]



This model was very slow to run (so slow we’re not going to try run it in the browser), as there were many species and thus more calculations needed to be done per loop. However, it primarily was slow as we needed a small dt value to keep the model accurate so there were more loops in total.

We had to use a very small dt value or errors could easily add up with so many species. It ended up taking about a minute to run each simulation which wasn’t acceptable for us, as we wanted to be able to quickly change parameters and ideally test a range of parameters. Speeding up the model would make it more useful to the wetlab team. They then could vary each of our parameters and see how they affected GFP output, for example seeing the effects of different switch designs and experimental setups.

Having lots of species also made it harder to find accurate rate constants for all the steps in the reaction. Tobias Stögbauer’s PhD was very useful for this, although many of his values weren’t that applicable to us as we used a different cell-free system leading to inaccuracies - which were amplified when dt was increased.

It was also difficult to get useful data from the model as there were so many parameters that changing just one was unlikely to have a significant effect on the final result. This made it harder to draw conclusions from the model. For similar reasons it would also have made it hard to adjust the parameters to fit to our lab data when we had it. More parameters also increased the total uncertainty, decreasing the precision of our model.

Streamlined mass action model

We met with Thomas Ouldridge again after the previous model and he showed us how we could combine many steps of the reaction to slim our model down to just 4 parameters and 5 species, as well as giving us some more accurate parameters for our system.

With fewer parameters it was much easier to see the effect that varying each parameter has, and thus draw meaningful conclusions from the relationship between the parameters and GFP output.

Given the knowledge of the system that we had, it was reasonable to incorporate:

  • the GFP maturation into ktranslation
  • the opening of the toehold switch into the binding of the miRNA

We used mass action kinetics to model the decay of the RNAs and the opening of the toehold switch, however we modelled transcription and translation such that DNA and the open switches weren’t used up. We assumed the pool of resources was constant, as our construct would not significantly deplete the cell free extract’s resources.[5]

We also made the assumption that our inducible promoter could be modelled as a constitutive one. Given arabinose concentration is constant and the assumption that the time taken to induce the promoter was negligible, no changes to the way transcription was modelled were necessary, other than neglecting resource degradation. This was confirmed experimentally when we characterized part BBa_K808000.

In this model we included a term for miRNA decay, assuming that it decayed at the same rate as the other RNAs in the system. This was on the recommendation of Thomas Ouldridge as the degradation of miRNA is not well documented in our cell free system.

As we were not confident in this assumption we tried running the model with different miRNA degradation rates and found that between 5×10-4 and 5×10-3 the amount of GFP produced in the positive case increased by a factor of 3. Outside of this range the difference was small. As our value is within this range the uncertainty introduced by the miRNA degradation rate was high. However, this did not affect the model's usefulness in understanding the qualitative effect of varying the parameters as the miRNA degradation rate stretches the GFP against time graph horizontally. This changes the time our test takes, but not the final difference in GFP output.

As with our previous model we assumed protein degradation was negligible as our cell free system does not contain any proteases.



Following correspondence with Oxford’s modeling team, we changed the values of some of our parameters. Oxford recommended we used Karzburn’s 2011 paper, ‘Coarse grained dynamics of a cell free system’ for the translation constant. They also pointed out that we had forgotten to change some of the parameters’ units so our model was not dimensionally consistent, and we updated the model to fix those.


Parameter name Where we got it Value
ktranscription Used Tobias Stoegbauer’s PhD for this value, as it is equivalent to his kts value of 2.2NtPs-1 1.1×10-3s-1
ktranslation Karzburn’s 2011 paper gave 4 amino acids s-1[6]. As the GFP sequence we used was 239 amino acids long,[7] we simply did 4aa s-1 ÷ 239 aa. 1.7×10-2s-1
kopen When we met Thomas Ouldridge at Imperial he suggested a range of 105 to 106 for this value, given his experience modeling toehold switches. 6×105 M-1s-1
kdecay Estimate of messenger RNA decay rate in our cell free extract.[8] 1.28×10-3


One of the biggest potential problems we identified was that if enough translation of closed switches occurred, it would be difficult to see any meaningful results. The kleakage parameter represents the ratio of translation rate on the closed and the open switch. We concluded that would not get any observable difference if translation on the closed switch was above a ten thousandth of the rate on the open switch, something we would never have been able to quantify easily without the model.

Without the model we wouldn’t have realized quite how crucial reducing leakage was. We spent a lot of time designing our switches to minimize leakage.

The lab team also asked us to find the optimal concentration of DNA for both output variation and overall output. By varying ktranslation and the concentration of DNA we found that a higher concentration was better in general, however it would plateau eventually as miRNA consumption became the limiting factor. In the specific case of miRNA degradation being 2 orders of magnitude lower than the rate of degradation of the other RNAs the GFP concentration plateaued around 5×10-9 mol dm-3.

Try out the model for yourself and see if you can see our conclusions by changing the parameters to see the final output of GFP here:




Results fitting and analysis

Having collected our experimental results we reevaluated our model. We scaled the results to fit them on the graph, as they were in arbitrary units of fluorescence. We can plot this fluorescence is proportional to GFP concentration.[9]


00.511.522.533.54×104 time / smiRNA = 9×10-9 MGFP concentration / M00.20.40.60.811.21.41.61.8miRNA = 9×10-8 MmiRNA = 9×10-7 MInitial model and results ×10-11


The major differences between the results and our theoretical values was the later maximum GFP concentration and the difference between the measured maximum and the theoretical maximum for 9×10-9moles/dm-3 of miRNA. In order to delay the peak in GFP concentration the rate of RNA decay had to be decreased. We also introduced a term for GFP degradation, to account for the peak in fluorescence. We then tried to fit the data by changing the parameters with these ideas in mind. The new values for these parameters were:

  • kdecay ≈ 3×10-4
  • GFPdecay ≈ 1×10-5


×10-10 00.511.522.533.54×104 00.511.522.5miRNA = 9×10-9 MGFP concentration / Mtime / smiRNA = 9×10-8 MmiRNA = 9×10-7 MRevised model and results


We varied the other parameters within our uncertainties. To make the model more accurate, especially for the lower miRNA concentrations we decreased kopen to 105, the lowest possible value in our range.

Another way to improve the lower miRNA concentrations was increasing the decay rate of just the miRNA. We had assumed it decayed at the same rate as the other RNAs in our system. Increasing just the miRNA decay rate had similar effect to decreasing kopen.

With our initial parameters our model was a poor representation of our system. Once we updated our parameters, it was able to represent our results relatively well.
  1. ↑ Chrome Developers (2013, September 23). Saying Goodbye to Our Old Friend NPAPI. Chromium Blog. Retrieved October 9, 2017, from https://blog.chromium.org/2013/09/saying-goodbye-to-our-old-friend-npapi.html
  2. ↑ (2017). Competition/Deliverables/Wiki - 2017.igem.org. Retrieved October 1, 2017, from https://2017.igem.org/Competition/Deliverables/Wiki
  3. ↑ CreativeJS (2013, June 3). The race for speed part 1: The JavaScript engine family tree. Retrieved October 8, 2017, from http://creativejs.com/2013/06/the-race-for-speed-part-1-the-javascript-engine-family-tree/index.html
  4. ↑ Stögbauer, S. (2012) Experiment and quantitative modeling of cell-free gene expression dynamics. Ludwig Maximilian University of Munich, Germany
  5. ↑ Shin, J., & Noireaux, V. (2012). An E. coli cell-free expression toolbox: application to synthetic gene circuits and artificial cells. ACS synthetic biology, 1(1), 29-41.
  6. ↑ Karzbrun, E., Shin, J., Bar-Ziv, R. H., & Noireaux, V. (2011). Coarse-grained dynamics of protein synthesis in a cell-free system. Physical review letters, 106(4), 048104.
  7. ↑ SnapGene (n.d.). eGFP Sequence and Map. Retrieved October 8, 2017, from http://www.snapgene.com/resources/plasmid_files/fluorescent_protein_genes_and_plasmids/EGFP/
  8. ↑ Shin, J., & Noireaux, V. (2010). Study of messenger RNA inactivation and protein degradation in an Escherichia coli cell-free expression system. Journal of biological engineering, 4(1), 9.
  9. ↑ Furtado, A., & Henry, R. (2002). Measurement of green fluorescent protein concentration in single cells by image analysis. Analytical biochemistry, 310(1), 84-92.