Team:BCCS-Bristol/BSim/Case studies/Repressilators

Repressilators coupled by quorum sensing - an agent based approach
Recent theoretical work has demonstrated that by coupling the individual dynamics of a number of synthetic genetic regulatory networks (GRNs) such as the repressilator [1], correlated population level dynamics, such as synchronisation can emerge [2]. Studying synthetic GRNs on an individual level can provide us with a deeper understanding of how natural GRNs function, however by expanding our perspective to a population level, the full functional capacity of GRNs can be uncovered.

What is a repressilator?
A repressilator is a synthetic GRN formed of one or more genes, where each gene represses its successor in the cycle. Depending on the number of genes, the network can exhibit sustained oscillations. An example of an analogue in the world of electronics would be a ring oscillator. The first in vivo construction of a repressilator was a cycle of three genes. By connecting the output of the repressilator to a green fluorescent protein (GFP) reporter, oscillations could be observed visually.



The individual repressilators can be coupled by extending the network to include a 'quorum-coupling module' (fig. 2) [2]. The coupling module contains an autoinducer (AI) component that allows individual repressilators to communicate via diffusion through the extracellular space.



Previous work
The quorum-coupled repressilator is an excellent initial model for investigating population level GRN dynamics. As a result there have already been a number of theoretical models produced to investigate the behaviour of such a system, demonstrating that population level synchronisation is achievable, although is is dependent on a number of factors, especially cell density [2].

However, all the previous studies of quorum coupled repressilators have assumed that the AI is well-mixed across the extra-cellular space. In reality chemical concentrations will vary across a large space, therefore it is important to consider the effects of a nonuniform extracellular chemical field on GRN dynamics.

Why we used BSim
By bringing together an agent based approach with the standard ordinary differential equation methods typically used for modelling GRNs, we hope to be able to extensively study spatial factors affecting the behaviour of coupled repressilators across a population.

BSim coupled repressilators
When modelling the population of repressilator-containing bacteria in BSim, we used many of the same features as the model outlined in [2]. The repressilators themselves are modelled as a system of 7 ODEs; 3 ODEs representing the 3 different mRNA levels, 3 ODEs for the 3 corresponding proteins, and one ODE for the internal level of autoinducer. These are in turn coupled to an external spatially varying chemical field via the autoinducer term, which incorporates physically correct diffusion and degradation characteristics. For the modelling of the autoinducer behaviour, we made the assumption that the autoinducer in question was AHL as this is a common quorum signalling molecule. The parameters for diffusion (in space and through the cell wall) and decay were then set accordingly. The parameter governing the ratio between mRNA and protein degradation was chosen from a random distribution as in the first part of [2].

Overview
Due to the complex nature of the complete model, there are a large number of possible directions for an investigation to take. As a result of time constraints of the iGEM project we were only able to look into the effects of a few of the parameters of the system. However we were able to verify the validity of our model by reproducing some of the key results in [2], which we summarise below.

However, the project continues, and we have outlined some possible directions in which we plan to take our investigation in the final part of this summary.

Results
So far, we have managed to investigate some of the effects of diffusion on the level of synchronisation in the system. Using realistic parameters for AHL, when considering the very small size of the simulation region, results in almost instantaneous diffusion. This is good as it allows us to simulate closely a 'well-mixed' environment as an initial assumtion and to evaluate the validity of our model. The following series of images shows the external chemical field oscillating as a result of the rapid diffusion speed. This can also be observed in the video at the end of this section.

The following video shows 200 bacteria swimming in a 100x100x100 micron volume. Each bacterium has a system of ODE's inside which model the essential dynamics of a repressilator. In this case, the individual repressilators are coupled via a diffusing autoinducer signal (in this case AHL). The colour of a bacterium represents the level of lacI mRNA in that bacterium (the lacI gene is one of the genes present in the repressilator); the bacterium will change from yellow to red as the internal level of lacI mRNA increases. In the example shown here, all 200 of the individual repressilators are initialised with random conditions, but quickly synchronise due to the effect of the AHL communication.

   

Due to the long time scale of repressilator oscillations the simulation was run for 10 hours in total; the video has therefore been sped up a fair amount and as a result individual bacterial motion may be difficult to see. However, a number of key features can be observed, namely the oscillation of individual repressilators and the pulsng oscillation in the external chemical field resulting from AHL diffusing into and out of the bacteria.

Of the two diffusivity parameters present in our system, it was found that cell wall diffusivity had the most noticeable effect on synchronisation. The effect of external diffusivity was negligible due to the small size of our simulation region. In the following sections we summarise some of the results that were obtained when cell wall diffusivity was varied.

Phase locking
The following images show the time series of a bacterium vs. 3 other bacteria for the final ~6 hours of the simulation. It can be clearly seen that at this point synchronisation has occurred, though phase locking is not always perfect! This is due to the randomness in our degradation ratio parameter. (Left to right: increasing degrees of phase locking):

Frequency coupling
The following three images illustrate the logarithmic power spectral density of repressilators coupled across a population of 200 bacteria. The x-axis indicates the bacterium and the y-axis indicates the frequency of oscillation. Colour variation corresponds to log amplitude of oscillations at a given frequency for a given bacterium.

From left to right, the images show increasing cell wall diffusivity parameter from 0.01 to 1.0 and finally to 4.0. This was to try to get an initial idea of roughly how the synchronisation level varies across the population as this parameter is changed.

As can be clearly seen, for a very low diffusivity parameter there is a lot of variation in the frequency distributions across the population, given that the simulation was performed over a long time period. As the cell wall diffusion coefficient is increased, we observe that the degree of frequency based synchronisation increases. This demonstrates that over a long period of time the effect of the external chemical field does actually cause the population to oscillate at the same frequency as a whole, even if phase locking may not be perfect.

Synchronisation transition
Finally, we analysed population level synchronisation across a wide range of values of cell wall diffusion coefficient. The parameter R used to measure synchronisation is the same as in [2]. R represents a general order parameter defined as the ratio of the standard deviation of the time series of the average signal to the standard deviation of each individual time series averaged over the range of bacteria. An R value of 1 would indicate a perfect level of synchronisation according to this measure, likewise R=0 indicates no synchronisation.

Clearly an increase in cell wall diffusion rate results in an increase in population level synchronisation according to a sigmoidal distribution. The left graph shows synchronisation transition for random initial conditions and the right hand graph shows synchronisation for uniform initial conditions. As expected, synchronisation is greater when starting from uniform initial conditions due to a small amount of synchronisation already being present at the start. It is expected that performing a range of calculations while varying our other parameters such as cell density throughout the simulation region would produce similar graphs.

Further work
As was mentioned in the overview, there are a large number of factors present in our system that it would be prudent to investigate to obtain a full picture of the repressilator's behaviour. With the full spatial chemical field model there are relatively few simplifications that can be made, therefore it is important that the parameters present in our system and their effect on population level synchronisation are fully profiled before any further changes are made.

Subsequently we would like to investigate a number of areas, the main ones being:
 * A further investigation into the effects of cell density when simulated with a full spatial model.
 * The effect of proper stochasticity (stochastic ODEs) in the system. At the moment our mRNA/protein degradation ratio is simply chosen from a suitable distribution, however it would also be possible to model the ratio using a stochastic ODE. As we have already implemented ODE solvers in BSim, these could be adapted and extended to perform integration of full stochastic ODE systems. It would then also be possible to model other aspects of the system in a more realistic manner, incorporating stochastic effects, for example stochastic diffusion through the cell membrane.
 * External control of the system. This could be performed at its simplest level by injecting a pulse of autoinducer into the extracellular field. This injected signal could be varied to investigate whether entrainment on a spatial level is possible. There are two main ways in which the system could be controlled:
 * Open-loop entrainment of uncoupled repressilators to some external periodic fluctuations of the autoinducer concentration obtained by additional inflow (or inflows) localized within the medium. Here the effects of diffusion are expected to be extremely relevant.
 * Closed-loop control of cell population via additional controlled inflow of autoinducer. This could potentially include the use of a virtual bsim population to decide the amount of autoinducer to be added/removed to a real population of bacteria in the wet-lab.
 * A better measure of synchronisation (integral order parameter perhaps) which incorporates spatial characteristics. This would also allow for much better analysis of effects such as localised synchronisation (clustering), which may arise as a result of spatial terms or non-uniformity in the cell population. One would also perhaps expect to see some form of travelling waves arising across the extracellular field and the population itself as a result of nonuniform cell populations and clustering.