From OpenWetWare

Jump to: navigation, search

EchiDNA 2014

Experiment 5 - Modelling of Cooperative Biosensor


The mathematical model of our cooperative biosensor builds on the model proposed by Bai et al. to explain the cooperative sensing and switching mechanism of the bacterial flagellar motor. Their work, in turn, follows on from the seminal paper by Duke et al. in the more general area of conformational spreading in cooperative multi-protein complexes. The code described in this section was written by Team EchiDNA in MatLab R2014a for Windows 7 64-bit on a i7-4790k (4GHz) system with 16GB system memory.


To simulate, in silico, the cooperative nature of our multi-subunit biosensor, and feed the results back into the design in order to iteratively optimise its behaviour.

Methods and Materials


The model relies on the state of each subunit of the assembled switch complex. Each subunit has two parameters giving rise to energetic differences that are taken into account:

  • Binding of the signal DNA to a subunit
  • Switching the conformation of the subunit from inactive (closed & non-fluorescent) to active (open & fluorescent)
The thermodynamically-favoured state of an unbound subunit is inactive.
The thermodynamically-favoured state of a bound subunit is active.

We will represent a as an inactive subunit and A as an active subunit. We will represent b indicating unbound signal DNA and B representing bound signal DNA.

The first thermodynamic value we consider is the free energy change of binding the signal DNA to a subunit, which we call EA, as per Bai et al. We have decided to use a symmetric model, in which EA is both the energy cost of AB to aB, and ab to Ab. In the limit of large EA, coupling between ligand binding and subunit conformational change is absolute.

Fig 1. Diagram based on Ma et al. Explanation of EA. Lower free energy is favourable. In the perfectly symmetric case, EC, not discussed here, is also zero. The above figure demonstrates that the preferred conformation for unbound sensor is closed, whereas a bound sensor prefers the open conformation.

The second thermodynamic value we consider is the energy cost of conformational change taking into account the conformation of the neighbours, which we call EJ, as per Bai et al. By design, our switch prefers neighbouring subunits to be of the same conformation. Therefore we introduce an energy cost of EJ for each nearest-neighbour of different conformation, and an energy cost of 0 for each nearest-neighbour of the same conformation. In the limit of large EJ, each subunit is forced to behave the same as its neighbours, and a simultaneous concerted switch is forced.

Fig 2. Diagram based on Ma et al. Explanation of EA. Illustration of EJ reflecting the energy changes upon switching conformations of a single subunit. Positive energy changes are unfavourable, and EJ is a positive quantity.

This means that for each conformational switch of a subunit of the switch, we have several possible energies depending on the local environment of the switch. These energies will be calculated throughout the simulation, and can hold six possible values: ±EA, ±(EA + 2EJ) or ±(EA – 2EJ)


In order to investigate how the model behaves as a potential cooperative DNA-sensing switch, the system was simulated using a kinetic Monte Carlo method.

In the standard Metropolis Monte Carlo method, no temporal information is gathered and the simulation is only interested in finding the final, equilibrium state of the system. The principle of Metropolis Monte Carlo is that subunits of a system are randomly examined and their states potentially changed. Each change of a subunit of a system has an energy cost associated with it. If the energy cost is favourable, the change is accepted. If the energy cost is unfavourable, the change is accepted with a certain probability dependent on the energy cost. This process happens many, many times before equilibrium is reached.

Kinetic Monte Carlo is similar, in that it randomly examines subunits of a system, but instead generates a random time in the future at which it expects the state of the subunit to switch or signal to bind. In our simulation, the random time generated is weighted according to first order kinetics, according to the equation t-t0 = -loge(rand)/k. k is the rate constant (calculated as shown below), t0 is the time at which the calculation is made, and rand is a random number generated using Matlab’s inbuilt rand function (using the multiplicative congruential method). The simulation steps forward in time, executing the next subunit switch or signal binding each step. In addition, each subunit switch is accompanied by a recalculation of times for its nearest neighbours as the energetics of the local environment are changed.

To do this we calculated rate constants using the Arrhenius equation for the conformational switch of the subunits and the binding of signal as shown below. In future models, we will seek to investigate the appropriate

Eq 1. Based on Bai et al. Explanation of ka. Expressions based on the Arrhenius rate equation for energy for subunit conformational change. In the symmetric model, we choose λa to be 0.5 as an intermediate value. Like Bai et al., we lack sufficient information about the thermodyanmics of our system to make a more informed choice. ωa is chosen to be 104s-1 in the model of Bai et al. which is a typical rate for conformational change for proteins.

Eq 2. Based on Bai et al. Explanation of kb. Expressions based on the Arrhenius rate equation for energy energy of binding of signal DNA. This rate is dependent on concentration for binding (but not for unbinding). However, the rate is not dependent on the conformation of the subunit in the model of Bai et al., and λb is chosen to be zero. ωb is chosen to be 10 s-1, identical to the model of Bai et al. because we lack information about the kinetics of DNA binding in our situation. c0.5 is the concentration for neutral bias of the switch (neither on nor off).


The base code used for Experiment 5 is available here.

  • n is the number of subunits in the switch.
  • c is the concentration of signal DNA as a multiple of c0.5 (which itself is the concentration for neutral bias for the switch - neither on nor off)
  • EA and EJ are as described above, referred to below in units of kBT

As an initial investigation, simulations were run with n = 30, c = 2 for several values of EA and EJ. Animations showing the typical behaviour of the switch ring were recorded and analysed. Below is an image demonstrating the typical features of a cooperative switch, with a nucleated domain. These domains have the potential to either expand and hence activate the entire switch, or to collapse, and hence return the switch to an inactive state.

Fig 3. A frame showing the switch from the MatLab output demonstrating the salient features of the visualisation.

The animations can be summarised in images that show the state of the switch over time to allow more efficient analysis and to be able to grasp the behaviour of the switch under different conditions as a whole. The images shown below show the fraction of the 30 subunits that are active in blue, and the fraction of subunits with ligand bound in red.

Click on the images below to see their associated simulation videos

Fig 4. EA = 1, EJ = 1, n = 30, c = 2. Notice the relatively chaotic behaviour, even though an overall trend is observed following the number of subunits with bound signal.

Fig 5. EA = 1, EJ = 4, n = 30, c = 2. With a large EJ, we can see that each subunit is strongly held in the same state as its neighbour. Hence, we see well-defined domains appearing and disappearing and a clear switching event.

Fig 6. EA = 4, EJ = 1, n = 30, c = 2. With a large EA, we instead see that binding of signal is strongly linked to activation of the subunit. Hence, the number of active subunits follows closely the number of signal-bound subunits.

Fig 7. EA = 4, EJ = 4, n = 30, c = 2. Whilst this may seem similar to Fig. 6, notice how early on, the red line is above the blue, but later, the red is below the blue. This is due to coupling of signal binding to activation, which is disfavoured due to the the high cost of creating a new active subunit in the middle of an inactive domain. Therefore red is above blue. In the later stages, enough signal has bound so that the system wants to "fill in the gaps" between active subunits and we get more active subunits than signal-bound subunits.
Based qualitatively on these results, we can appreciate that for our cooperative sensor, we have a small window of ideal cooperative behaviour as shown in Fig. 5. To accomplish this, it is apparent that we need a larger EJ than EA. To further investigate the effects of EA and EJ on the behaviour of the switch, we ran simulations to examine how long it would take switches with n = 30, c = 1 to reach a state of all subunits active as a measure of how efficient and sensitive the switch was. Each value is an average of 10 runs. The result is shown below:

Fig 8. EA from 0.25 to 5 in 0.25 increments, EJ from 0.25 to 5 in 0.25 increments, n = 30, c = 2. Each block is an average of 10 runs.

The results reflect the previous simulations. Switching is achievable for certain values of Ea and EJ. In particular, notice the only scenario where a switch is achieved in reasonable time is for low EA and a moderate value of EJ. Because this particular set of simulations is so vital to understanding the time-scale of the switch, we re-ran the simulation with more energetic resolution. This revealed a clearly defined region as before where the switch appears to operate on a reasonable time scale. The result of 56 hours of calculation is shown below.

Fig 9. EA from 0 to 12.5 in 0.1 increments, EJ from 0 to 25 in 0.1 increments, n = 60, c = 10. Each block is an average of 10 runs.

As before, we see that low Ea and moderate EJ values are required. In particular, we note that the general shape of the allowable area is the same as for the previous images despite using different n and c values. This supports the universality of effect of both EJ and EA.

This graph is informative, but unfortunately does not give us much information about the specificity. Does the switch activate randomly without the presence of signal? In order to determine the behaviour of our switch, it is necessary to characterise how the switch responds to changing concentrations. Therefore, we ran a series of simulations with n = 12 in which EA was held constant at 1, whilst concentration was scanned from 0.2 to 5 in increments of 0.2 for a range of EJ values ranging from 0 to 5 in increments of 0.1. The information obtained was the fraction of active subunits after 1 second of simulation time. Each data point obtained was computed as the average value of 10 runs. The results of these simulations are shown below:

Fig 10. n = 12, EA = 1, EJ from 0 to 5 in 0.1 increments, c from 0.2 to 5 in 0.2 increments. Each data point is an average of 10 runs.

From this set of simulations, we can see that when EJ is 0, poor cooperativity is observed (Purple line, Fig. 10). However, increasing EJ to about 1.2, we get the onset of cooperativity which strengthens as EJ further increases. To quantify this, we attempted to curve-fit each of the response-concentration to the Hill equation:

Fig 11. n = 12, EA = 1, EJ from 0 to 5 in 0.1 increments, c from 0.2 to 5 in 0.2 increments. Each data point is an average of 10 runs.

θ is the fraction of active subunits, c is the concentration of signal DNA, Vmax is the maximum value, k is the dissociation constant and n is the Hill coefficient. We did not force Vmax to be 1 as it is possible that for low EA or EJ the switch may statistically never approach a completely active state over finitely many trials. This method resulted in the following Hill coefficients being observed.

Fig 12. Data was calculated under the conditions for Fig. 11. Non-linear curve-fitting was performed in OriginPro 8 using the inbuilt Hill Equation specification.
Note that the Hill coefficients follow a generally increasing trend. There is larger spread of the hill coefficients at higher EJ since for such tight curves more data points and more runs to average over are required for smooth fitting. Small anomalies in the inherently random simulations will result in noise being added. Almost all of the curves with EJ > 1 were fit with an adjusted R2 > 0.95 indicating reasonable fit. The other two parameters were not of interest to us.

Finally, we decided to run simulations to assess the overall effect of all four of n, c, EA and EJ to observe trends in how the switch behaves. From these graphs, trends can be observed that tell us what we should do in order to increase the sensitivity of the switch. Or the specificity of the switch. What is the optimum size of switch? These are all pertinent questions in our design. The result of this monstrous calculation is shown below.

Fig 13. EA and EJ values are indicated on the graph, c from 0.2 to 5 in 0.2 increments, n from 2 to 50 in increments of 2. Each data point is an average of 10 runs.

These simulations provided sufficient information for us to qualitatively tune the behaviour of our cooperative switch in an iterative fashion. Experimental data from each design are able to be fed into the simulations and allow us to refine the initially unknown coefficients in order to create a more accurate model. This in turn allows us to devise a better combination of energies and kinetic parameters that can be used to rationally design the sequence of our DNA to approach as closely as possible the values computed in silico.


After much modelling of a cooperative switch in silico, a lot of valuable information was gained about how the behaviour of the switch depends on several parameters that we can rationally control and incorporate into the design of our cooperative biosensor:
  • EJ was shown to regulate neighbour interactions and EA was shown to affect correlation of signal-binding with subunit activation, as intended in the model.
  • We identified that a low non-zero EA and a moderate EJ is optimal for the intended purpose of the switch.
  • Finally, we confirmed that cooperativity emerges in our switch under the correct conditions as shown by analysis of the Hill coefficients.
This forms the basis of our proposed iterative approach to switch design. We have successfully constructed a model that is both intuitive and reproduces both previously known information, as well as providing new insight into the anticipated behaviour of our biosensor.

Looking ahead at further work with the model, we propose that single molecule fluorescence could be of great importance in obtaining experimentally verified kinetic parameters for the rate constants underlying the mathematical model. This would allow us to better model our specific system instead of a generalised cooperative switch.

In addition, many important questions remain unexplored! We would like to understand the scaling of computational time with increasing size of the cooperative ring. We believe that experimenting with a kinetic Monte Carlo algorithm would further drive our discoveries about this model, our biosensor, and the phenomenon of cooperativity in nature.

Click here to go back to the Lab Book Overview

Personal tools