Team:KU Leuven/Project/Oscillator/Modelling

We transfered our oscillation model into a preliminary ODE system composing a rate equation for each of the six components of this circuit. We expanded upon this system by adding in the fact that the production of quorum-sensing molecules goes through an enzymatic step and by using an approximation of the OR gate that is easier to realize in the wet lab.

Then we wanted to define a parameter space for which these ODE systems produce oscillations. This was done by random probing a realistic range of the different parameters to test if sustained oscillations take place. The next test we performed was testing the synchronization and rapid synchronization of the most complete system. Each of the two tests were conducted with a multitude of oscillating parameter sets. We tested the synchronized oscillation when cell-to-cell variability is existent by assuming a colony of ten different cells of which each cell has different parameters. The rapid resynchronization was tested by using again a colony of ten different cells, however a number of the cells began with different initial conditions.

Our modeling showed that we had designed a very robust oscillator:

1. Our system produces sustained oscillations for a very broad range of parameters.
   
2. The oscillator can withstand a high amount of cell-to-cell variability.
   
3. Rapid resynchronization occurs even when 40% of the cells start out of sync.

To conduct our computations we made use of the Flemish Super Computer Centre (VSC), which is a partnership of the five Flemish university associations. To use this we first had to rewrite the tests so they could be run in parallel on the supercomputer cluster at our university.

For the transformation into a biological network, we refer you to the wetlab page.

1.1 Exploring our possibilities

Firstly, it is important to clarify what we exactly mean with ‘parameter’. A biochemical system has a very high variability in ranges for transcription rates, translation rates, degradation rates, etc. On top of that these parameters are not perfectly quantified and are subject to changes in the conditions. We checked what part of the parameter space (Box 1: Explanation of parameter space) creates oscillations. There are possibilities for doing this in a purely mathematical manner. Tyson (2002) gives a good example of how to study systems that produce biochemical oscillations. Polynikis, Hogan and Bernardo (2009) described modelling approaches for gene regulatory networks more generally. This is typically done by investigating the eigenvalues of the Jacobian matrix of the system. This becomes increasingly more difficult when the number of parameters and variables increases. A high level of non-linearity complicates the study of the behaviour even further. In order to see what is possible we contacted Professor Dirk Roose, an expert in non-linear systems analysis. We explained him we want to investigate what parameter values create a synchronized oscillation. However, our parameter space consists of about 20 parameters that can each vary with more than a factor ten. On top of that we will have highly non-linear equations. Professor Roose told us this amount of variability would make a clean-cut mathematical examination of our model impossible . Since it is not possible to reduce our parameter space, without diverting from our goal of fully studying our system, we decided to use another strategy. We will study this enormous parameter space by generating random sets of parameters throughout this space. This offers a less theoretical, nonetheless effective means of assuring this model robustly produces oscillations.

1.2.2 Enzymatic production of quorum-sensing molecules

As mentioned on the design page the quorum-sensing molecules are produced by an enzymatic step. Equation 2 shows the equations that result in that behaviour. The transcription factors control the production of the enzymes rather than A and B directly. The production of A and B is modelled by using Michaelis-Menten kinetics (S_A/(K_A+S_A)). To not overly complicate our model we will assume that the substrates S_A and S_B are at saturation levels. This means they are a lot higher than K_A, K_B respectively and the Michaelis-Menten consequently returns 1. This resulting system of equations is shown in Equation 3.

Equation 2: ODE system exhibiting the behaviour of the oscillator, including an enzymatic production step for the QS molecules
Equation 3: ODE system exhibiting the behaviour of the oscillator, including an enzymatic production step for the QS molecules with the assumption of saturation of those enzymes

1.2.3 Alternative OR gate

In the part before we used an actual OR gate, however, in the wetlab we use an approximation of a real OR gate. We implement the gene for production of C (D) twice, once with a promoter that is induced by A (B) and once with one that is induced by X (Y). This is easier to implement although the behaviour is somewhat different. There is now a higher expression when both the inducers are present. There is no difference in the fact that this displays a asymmetric time delay, of which the importance is discussed in the design page. This altered circuit is displayed as Figure 3.

This is translated into Equation 4 by using the format instead of . It is clear that this will return 1 in the case that either A or X are present, 0 when both are absent and 2 when both are present.

1.3.1 Selection of the parameter range

For the degradation rates of the proteins we assume they follow the N-end rule and are not in one of the unstable categories. This means the dilution rate controls the effective degradation rate and thus we use values that are appropriate for that. More information on how these values have been found can be found on the methyl salicylate modelling page. The ranges can be found in Table 1.
For the colony-wide molecules this approach for the degradation rate cannot be used, since they are present inside as well as outside of the cells. Those molecules thus become diluted when the volume of the colony increases, but not per se when the amount of cells increases as is the case for stable proteins. When we think about using quorum sensing molecules, the colony-wide molecules are metabolites. A potential problem might then be the slow degradation rate as encountered by Danino et al. (2010). On the wetlab page, we describe the use of an actively degrading enzyme as a solution for this issue of slow degradation. Since our model does not require a specific chemical species for oscillating, we will allow a broader spectrum of possibilities than through exactly that degrading enzyme. We will simplify the model by assuming a uniform degradation throughout the colony, but by allowing two decades of variance, we will still allow for a broad spectrum of possibilities. We will use the degradation (dilution) rate of the proteins as a mark for the degradation rate of the colony-wide molecules, since the volume of a colony also increases when the amount of cells augments. The two decades of variance in the degradation rate of colony-wide molecules will range from one tenth of the degradation rate used for proteins to ten times that degradation rate. Allowing values lower than the degradation rate used for proteins is based on the fact that the growth of the volume of a colony is slower than the growth of the number of cells. Allowing values higher than the degradation rate used for proteins is based on the possibility that the colony-wide molecules are for instance actively degraded by some produced enzymes or a fast spontaneous degradation. These values can also be seen in Table 1.

The maximal production rates (β) of the proteins can vary enormously, for instance due to differences in copy number. However, the goal of this oscillator is to have the components in the chromosomal DNA of the cells. We thus assume values that are in the rage of production rates of chromosomal genes. For those ranges we use the assumption that the steady-state concentrations of proteins is equal to β/α (see Box 2), we used the α’s as determined above and we use protein abundances as determined by Ishihama et al. (2008). They determined a steady-state concentration for enzymes to be around 65 per cell and for components of the biosynthetic machinery to be around 105. An important remark is that the assumption that steady-state concentrations are equal to β/α is only correct for as long as the production is fully on and not when there were some inhibitory effects active. This means that the lowest protein abundances that are determined by Ishihama et al. can be a lot lower than β/α due to for instance regulation of the production rate, which was consequently not at its maximum. When the production rate is only half of its maximum for instance, the steady-state concentration of the protein is only equal to 0.5∙β/α instead of β/α, which means the lowest values determined by Ishihama et al. be too low. On the other hand the highest values might also be too high, since we are not dealing with components of the biosynthetic machinery. It thus seems wiser to us to discard the extremes and only probe between 6.5∙10² and 10⁴. We will use these values and transform them to production rates by using the formula from Box 2 and the degradation rates from Table 1. The calculated values can be found in Table 1.

The maximum production rates due to the enzymes cannot be determined this way, but also span an enormous range. Shaefer et al. (1996) found proof, however, that enzymes that produce quorum-sensing molecules have a low maximum production rate, and they found specifically for the Vibrio fischeri LuxI protein that β equals 1.1 molecule per minute. Since our model is not confined to one enzyme, we will again test a bigger scope. We will again use a range of two decades for this parameter, as can be seen in Table 1.

The final parameters we will look into is the threshold value. Our method will consist of randomly picking production rates degradation rates of metabolites. As mentioned in Box 2 the steady-state concentration at full production capacity equals β/α for proteins and , which means there is no possibility of having a higher level of that protein or metabolite. For the proteins we will use this as an upper bound for setting the possible threshold, since the threshold otherwise would never be reached. We will assume the threshold is reached somewhere in the process, but this is not necessary to have production as long as the Hill coefficient is smaller than infinity. As a lower bound we will allow one decade of variance. The formulas used are displayed in Table 1. The steady-state concentration at full production capacity of metabolites can however vary significantly from the thresholds of the proteins that act as transcription factors. Since the threshold is in reality a function of the interaction between promoter and transcription factor and not of the steady-state concentration of the transcription factor, it would be unfair to use such different values for thresholds for proteins on the one hand and thresholds for metabolites on the other. Instead we will use an the median of all β/α of the proteins for the thresholds for metabolites. Since we use a uniform probability distribution for both α and β the median will be the . The resulting value is shown in Table 1, with a lower bound that is one decade lower.

1.3.2 Setup of the in silico experiments

We have already skimmed the fact that we will use random probing of the parameter space in order to establish the fact that sustained oscillations take place. Here we will discuss how we exactly probed the parameter space, what we took as initial conditions for the ODE solver and how we tested for oscillating behaviour. Since we cannot assume any probability distribution, we will stick to a uniform distribution for each of the parameters. Each of the parameters will be allowed to vary independently from one another, except for the protein degradation rates and the metabolite degradation rates. Those will all be chosen to be the same since the value is based on the growth rate of the colony, which is of course the same from the perspective of different proteins within that same colony.

The initial conditions determine quite a lot of the behaviour of the network when there are non-linear functions involved. In order to test whether the system behaves as predicted we will use initial conditions that correspond to the explanation on the design page. This means we start with a high concentration of A and zero of each of the others. It is now only necessary to determine what qualifies as a high enough concentration. For this we use the maximum of K_AX and K_AC as a point of reference because that can be seen as the threshold for A to have an impact. We simply chose twice that threshold, which should suffice to spark the system.

Detection of oscillating behaviour was done differently for tests involving ODEs and PDEs. The latter had to be checked manually, because of their often very different appearance. The ODE tests however were checked on stable oscillations by a self-written script. This program takes the Fourier transform of the signal after the transient part using the Matlab command for Fast Fourier Transform (FFT). The the components of this FFT are split into the zero frequency component and all other components. Then, the sum of the absolute values of the components with a nonzero frequency is compared to the absolute value of the zero frequency component. If this ratio is greater than one, the signal is marked as oscillating. Practically, this means that the signal has to have some periodic components that are of the same order of magnitude as the zero frequency component, which is the average of the signal. Signals with negative values were automatically discarded due to their physical impossibility. If any of the signals in a test run was oscillating, the whole test run was marked as positive.

As an output we will determine what the oscillating fraction of the tested parameter space is. This is a good measure to find out how easily the system reaches oscillations. Since we randomly probe throughout the parameter space and find a binary result of either oscillating or either not. As this is a Bernoulli trial process, we can calculate the 95% confidence interval of the fraction of the parameter space that oscillates. The formula displayed as Equation 5 gives that 95% confidence interval and according to Brown Cai and DasGupta (2001) this formula should function well as long as both n∙p > 5 and n∙(1-p) > 5. This means that only in the case of a fraction close to 0 or close to 1 there can be a problem. However, when it is that close to either of the extreme values, the results are already clear enough because it simply means the system hardly ever oscillates, or oscillates almost all the time. When p is as close to 1 that n∙p > 5, a maximum 4 out of 100 tests produced a non-oscillatory system and we can definitely say that the system easily produces oscillations in that parameter space.

1.3.4 Testing the other ODE systems constructed earlier

We made confidence intervals for the fraction of the parameter space that has oscillations on the other representations of our model as well. We again used the parameter ranges from Table 3 and let each of the parameters vary independently from one another.

The first of the systems we test here is displayed as Equation 3. Here we add enzA and enzB which are responsible for the production of A and B (that are now metabolites instead of proteins). We do this in order to come closer to the proposed implementation of our model. Another thing that is important to mention here is the fact that this extra step means an extra delay in the system. This is useful because in an oscillating system the different peaks should not have too much overlap. This perhaps counters the fact that we did not resort to DDE’s. This should benefit the occurrence of oscillations. The entirely different ranges of production and degradation of the metabolites (here A and B) in respect to the proteins that act as transcription factors (here X, C, Y and D) may put a strain on the occurrence of oscillations, because this makes their ranges of concentrations also differ enormously. This means that problems may occur when we use the same order of magnitude as thresholds, because we determined earlier that they should at least be of reachable height. On the other hand, we do not allow pre-determined differences in parameters (e.g. using entirely different threshold ranges for proteins than for metabolites), since it has been our purpose all along to create a system that oscillates for a broad range of different in vivo components (e.g. promoters). It would thus be unfair to already determine that thresholds for metabolites have to be 100 times as high as for proteins for instance, because the mechanism they work with does not differ per se.

We did the same test as earlier in order to know the fraction of the parameter space that produces oscillations and the results are shown in Table 5. These results show a significant increase with respect to those in Table 4, with fractions that are close to 1, indicating that this system is an oscillator. This proves that the extra delay benefits the oscillating behaviour of the system.

In the next step we added an alternative version of the OR gate, of which the resulting system is shown as Equation 4. The resulting confidence intervals are shown in Table 6 and exhibit no significant differences with the results in Table 5.

Each of the models show a very frequent occurrence of oscillations, proving that our proposed system inherently produces oscillations. A recurrent phenomenon was the big difference in amplitude of the different components. However, we chose not to go in depth on this topic and which component’s oscillations had the highest amplitude was not always the same. A second noteworthy feature is the type of oscillations that occur, namely relaxation oscillations. This is a type of oscillations in which the production shifts when thresholds are reached. This is exactly what Hill functions represent and thus circadian rhythms, a widespread biological oscillator, often exhibit this type of oscillations (Maheshri and O’Shea, 2007).

1.4 Conclusion

As required our system produces sustained oscillations for a very broad range of parameters. We simply used broad ranges for each of the parameters and the oscillations followed easily. There was only one parameter our used ranges were not backed up that strongly, namely the threshold value K. However, the system produced oscillations as long as the value was reachable, so its value is not crucial for the production of oscillations. Of course the production of oscillations in its own is not an impressive feature, since many known networks do that, as for instance the important example proposed by Elowitz and Leibler (2000).

The features of the proposed system that result in synchronization and rapid synchronization are undoubtedly the most important. An explanation of how this exactly happens can be found on the design page. Here we will show how we tested these properties. When we tested the parameter space for oscillations, we could stick to one cell only. Now, however, we will need spatial heterogeneity, which expresses the fact that there is ‘empty’ space in between cells.

In order to test whether synchronized oscillations occur, we will use a colony of ten different cells of which each has different parameters. This allows us to test whether synchronized oscillations occur even though there is cell-to-cell variability. This phenotypic variability can for instance arise due to differences in their position in the cell cycle and thus in overall protein production rates. A second feature we will test here is the occurrence of rapid resynchronization. This will be tested by using again a colony of ten different cells, however the difference will now lay in the level of each of the chemical species (A, X, C, B, Y and D) inside the cell. This expresses the differences that arise due to stochastic effects (Kærn et al., 2005). This is relevant because at any time a cell can express each of the proteins at an entirely different level than the others. We require that even then they rapidly resynchronize, which is an inherent feature of our model.

Each of the two tests will be conducted with a multitude of oscillating parameter sets which we received from the procedure above. As for the reaction rates we will stick to Equation 4 since this is the most complete system and it is this one we are eventually interested in. To test whether oscillations occur we use the same strategy as before.

2.3 Rapid resynchronization

The last test we conducted is testing whether rapid resynchronization occurs. We do this by having completely different initial conditions for a number of the cells and test whether oscillations still occur. We will each time let X and C be four times as high in some cells and Y and D be one fourth compared to the other cells. We chose to let 40% of the cells have these different values since this is a very extreme test and we want to seek the borders of what’s possible. In Figure 8 and Figure 9 we show two examples in which a different group of cells started out of sync. These examples show that there is no real difference in speed of resynchronization and we will simply take the setup of Figure 9 in our tests. Important to note is that the resynchronization happens incredibly fast, which is very beneficial for an in vivo oscillator.

We’ve tested this for 100 randomly probed parameter sets in each of the threshold ranges discussed earlier. The results of this test are shown in Table 8 and except for the second threshold range, the confidence intervals overlap with those of Table 6, which shows the oscillating fractions when only one cell is studied. This, together with the very high fractions of oscillations shows that the studied system exhibits rapid resynchronization, as is crucial for any colony wide oscillator.

Within the range of reasonable parameters, a high fraction resulted in oscillatory behaviour. This is required for three reasons. First of all this means that a broad range of components can be used to create our model in vivo, which was our goal all along. Secondly it is important to have oscillations occur in a number of different cell densities, since this influences the metabolism of the cells in many ways. Example given: the production of the colony-wide molecules is spread out over a bigger volume when cell densities are low than when they are high. Thirdly, this means that oscillations will still occur even if parameters change due to an extra load on the cell’s system, which is for instance the case if we couple a production module to our oscillator.

The oscillator clearly can withstand a high amount of cell-to-cell variability, which is an aspect that is obviously very important in synchronized oscillations. The rapid resynchronization is also a feature our model has and this is important in the context of the stochastic effects desynchronise the cells. Our results also showed that the resynchronization occurs incredibly fast, which is very beneficial in order to have sustained oscillations.

The several tests conducted on our model gave satisfactory results and indicate our system creates synchronized oscillations. However, extrapolating this to saying it will also work in vivo is a bridge too far. This is because there is an enormous amount of complexity we did not take into account. We did use the most widespread technique for studying genetic networks, namely through ODE’s (Polynikis, Hogan and di Bernardo, 2009) so our results show it is definitely worth testing it in vivo. For this test the first steps have been taken in our wetlab page about the subject.

Brown, L. D., Cai, T. T., DasGupta, A. (2001). Interval Estimation for a Binomial Proportion. Statistical Science, 16(2):101-133.
Danino, T. et al. (2010). A synchronized quorum of genetic clocks. Nature, 463:326-330.
Ishihama, Y. et al. (2008). Protein abundance profiling of the Escherichia coli cytosol. BMC Genomics, 9:102.
Kærn, M., Elston, T. C., Blake, W. J., and Collins, J. J. (2005). Stochasticity in gene expression: from theories to phenotypes. Nature Reviews Genetics, 6:451-464.
Elowitz, M. B., and Leibler, S. (2000). A synthetic oscillatory network of transcriptional regulators. Nature, 403(6767):335-338.
Maheshri, N. and O’Shea, E. K. (2007). Living with noisy genes: how cells function reliably with inherent variability in gene expression. Annu. Rev. Biophys. Biomol. Struct., 36:413-434.
Park, S. et al. (2003). Motion to form a quorum. Science, 301:188.
Polynikis, A., Hogan, S. J., and di Bernardo, M. (2009). Comparing different ODE modeling approaches for gene regulatory networks. Journal of Theoretical Biology, 261:511-530.
Rosenfeld, N. et al. (2005). Gene Regulation at the Single-Cell Level. Science, 307(5717):1962-1965.
Schaefer, A. L., Val, D. L., Hanzelka, B. L., Cronan, J. E. JR., and Greenberg E. P. (1996). Generation of cell-to-cell signals in quorum sensing: Acyl homoserine lactone synthase activity of a purified Vibrio Fischeri LuxI protein. PNAS, 93:9505-9509.
Stewart, P. S. (2003). Diffusion in biofilms. Journal of Bacteriology, 185:1485-1491.
Tyson, J. J. in Computational Cell Biology (eds. Fall, C. P., Marland, E. S., Wagner, J. M. and Tyson, J. J., 2002), 230-260 (Springer, New York).