Team:KU Leuven/Project/Oscillator/Modelling

In this part of the wiki we will describe how we performed an analysis of our proposed oscillating system and the results. This text starts with a small introduction on what we want to achieve with this oscillator, a topic that is more thoroughly elaborated on the design page. Before we started the analysis that is stated here, we looked up how similar networks have been analyzed before in order to see what direction we will take. A full-scale analysis would go beyond the scope of the project so we will stick to an elaborate indicative study. The first step is to translate our network into ODE’s (ordinary differential equations), which we will make more realistic step by step. We will use these to see how easily sustained oscillations form. This is of course not the most impressive feature since there are many known networks that easily produce oscillations. We chose not to include the effect on amplitude and frequency, since that would make the scope of this study explode and the many assumptions we have to make, render it unrealistic. The important feature of our network is its synchronization features. In order to check whether our model achieves rapid resynchronization, we will solve systems of PDE’s (partial differential equations). Since those are a lot more computationally intense, we used the Flemish Super Computer Centre (VSC) in order to do our computations time-efficiently. The explanation of how the network functions and how it attains synchronized oscillations can be found on the page explaining the design. 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) 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.

Box1 | Parameter Space

We will explain this concept with a simple example. Let’s say we want an iGEM team (our system) to develop an oscillator. Two parameters of our system will determine whether this endeavour is successful. These parameters are the number of engineers in a team and the amount of coffee they consume per engineer. As to determine the parameter space we have to find relevant ranges for each of the parameters. We assume more than 10 engineers in the iGEM team is unlikely and a coffee consumption of more than 1.2 liter per day per engineer can also be discarded. The resulting parameter space for this two-parameter model can be easily visualised by a rectangle. This resulting parameter space can be probed and the parameter values we obtained were then used in a simulation which gave as an output whether an oscillator is produced or not. In the figure below we show the result of such a random probing of the parameter space. The blue dots represent unsuccessful endeavours and the orange dots successful ones. The results allow us to hypothesise on what is necessary to produce an oscillator. Apparently a higher amount of engineers still makes success possible, but only when the coffee consumption is not too high. We hypothesise here that this might be because the engineers get too much energy and start to fight amongst themselves. We used a similar approach in the testing of our model, with some distinctions though. We deme parameter sets positive when synchronized oscillations take place. But we have 26 parameters and that makes it impossible to visualise our positive parameter space. So we chose another way to look at the parameter space and look at the total fraction of the parameter space that exhibits the wanted behaviour.

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.

Box2 | Steady-state concentrations

A steady-state concentration of proteins can be easily determined by setting each of the rates to zero. In the case of maximal production rate, the rate equation has the following form: and when , the equation can be easily converted to , which gives the steady-state concentration of Z, or , when all but the production rates are known.
In the case of metabolites, that are produced by enzymes the rate at full production capacity equals which becomes at steady-state. When the enzyme is at its own steady-state concentration during maximal production rate this becomes .

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 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.