"
Team:ULB-Brussels/Project/Mathematical
From 2009.igem.org
Contents |
Mathematical modeling
Introduction
Synthetic biology is a recent development in biology which aims at producing useful material via biological agents. In this context, a biological system can be seen as a complex network composed of different functional parts (\cite{uri}). Mathematical tools allow one to make prediction on the dynamical behaviour of a given bio-system. We aim at studying the system which has been built in the second section of this report. We already know from our experiment that our system is able to produce glue in the presence of IPTG inductor. We now address the following question: what is the influence of the different parameters on the global dynamics (degradation rate, production rate, level of initial quantity for the C2P22 repressor,...) with and without IPTG. In order to achieve this task, we will consider three different models: the first one shows the basal property we are looking for: a system able to produce glue. The two next models can be seen as two different steps of improvement. The purpose of these models is to improve the experimental control of the glue production. Parameters values result from the literature and previous iGEM team wiki's. We know these remain qualitative and that lab work should be carried out to specify them more precisely.
We start from the global biological system established in the first section (cf. Fig.\ref{fig:without199}(a)). In order to be able to obtain the dynamical equations for the global behaviour of our system, we have to make some assumptions and simplifications. Firstly, we completely neglect the detailed composition of the biobricks. Each of them is considered as a simple block. The next step of our mathematical design is the modeling of the interactions between all the different blocks. The regulation network of our system is modelized using Hill functions for a repressor and an activator. \cite{uri}
- Hill function for activator :
- Hill function for repressor :
In these expressions, [x] is the concentration of activated genes, p is called the Hill coefficient, k is the activation coefficient and α is the maximum expression level of the promotor.
For each block of the system, we can obtain a dynamical equation by considering its interactions with the other blocks of the system. For each block we can build an equation of the following form:
Where, [x]i is the concentration of the gene i, R(Ha([x]),Hp([x])) is the regulating function which is a combination of Hill functions. The second term of the right side is the destruction term, where γi is the maximum destruction rate of the gene i.
We must bear in mind, however that the robustness of a given operational regime with respect to external perturbations strongly depends on the value of the Hill coefficients. \cite{arn} In particular, the robustness is expected to increase with the value of the hill coefficient. The cooperativity behaviour is also a function of the Hill coefficient. For these reasons we will consider a situation for which p = 2.
A first simplified model
In this first step, we present a simplified model of our system. We make the following assumptions: we consider that the LuxR+HSL complex is formed quickly at the beginning of the dynamics. This assumption, allows us to modelize the quorum sensing system by considering the complex LuxR+HSL only. The simplified schema is shown on Figure \ref{fig:without199}(b): the effect of the LuxI and LuxR is represented by the autoregulation arrow on the box of the complex LuxR+HSL. We also neglect the effect of the block parE in this first approach.
For the system shown in Figure \ref{fig:without199}(b) we can obtain the following equations (writing here [L] for the LuxR+HSL concentration:
- Equation for the c1 repressor block (designed by c1 in the equation):
In this equation, the parameter β has the following explicit form:
In our case, [LacI] can be considered as a constant, then we have β = β([IPTG]). In this first model we consider only the situation with IPTG inside the system: β ≠ 0. - Equation for the c2 repressor block (designed by c2 in the equation):
- Equation for LuxR+HSL block (designed by L in the equation):
- Equation for the glue production (designed by Gl in the equation, or by hfsGH in the text)
We choose a set of values for the different parameters:
Stationary state
The first step of our analysis is the study of the stationary point. In order to do that we consider the following algebraic system:
By solving this system of equations for the previous set of parameters we can find different sets of stationary states. So we have (considering only the states with real values):
Linear stability analysis
To check the stability of the stationary states, we perform a linear stability analysis. To this end, we assume that the configuration of the system is perturbed by:
ω is the growth rate of the perturbation.
is the vector representing the perturbation of each concentration. Linearizing the equations, we can rewrite them in the following form \cite{nico}:
(11)
where
and
is the Jacobian matrix. Taking into account the form of the perturbation (\ref{a}), we can write the equation (11) as follows:
and the eigenvalues ω satisfy the characteristic equation:
Given the roots ωi of the above equation, we can easily deduce the steady state. From (\ref{a}) it is easy to see that a given state is unstable if any of the ωi has a positive real part. In our case we obtain the following expression for 
:
Equation (14) is a polynomial of 4th order and its first root is of the form:
The others are the roots of a polynomial expression of the third order. We could obtain these roots analytically but their expressions are too complicated to be useful. So, we solve numerically the equation for the previous sets of parameters and we put into the equation (15) the corresponding values for the different sets of stationary states. We have:
- set one =
- set two =
- set three =
Two states are stable and the last one is unstable. The system can stay near its initial configuration, which corresponds to the second set of values or it can also go to another state for which we have glue production (first set). To learn more about the dynamics of this system, we need to study its complete behaviour.
Global dynamics
In order to have an overview of the dynamics of the system we perform a numerical integration of the dynamical equations with the set of previous values for the different parameters. When we solve the system, we see that the system travels from its initial configuration to the nearest stationary stable state which is the second set, corresponding to a state without glue production. We have to stress that we do not find any different situation by varying the set of parameters. We have always a stationary state with no glue production: the system chooses this final configuration. In order to force the system to go to the other state, we add a little basal production term (κ) in the dynamical equation for the complex LuxR+HSL:
noindent In this case, as we can see on the Figure below, the system reaches the first stationary state, for which we have a glue production.
 Dynamical evolution of the concentrations |
 Dynamical evolution of the glue concentration |
Discussion
In this first attempt to obtain a mathematical model of our biological system, we have been able to build a model for which we can obtain stable stationary states with or without glue production. The problem is the following one: the system shows a sensitiveness to the values of the different parameters and the only way to obtain the glue is to introduce a basal production term. The presence of this term is not unrealistic but it is clear that an improvement of this model is needed in order to obtain a dynamical system for which the switch from a state without glue to a state with glue is regulated by the value of β, meaning by the concentration of the inductor in the system (which is the control parameter from the experimental point of view). This is not the case with this simplified model because with κ ≠ 0 we have a glue production independently of the value of β.
Detailed model for the Quorum Sensing system
The aim of this section is to present a more realistic system in order to solve the problem pointed out in the last section. The main difference in this new model is the following one: as we can see on the figure below, we consider a more complete description of the quorum sensing system by taking into account the presence of the LuxR and LuxI blocks. We also consider the effect of the parE block.
For the system shown on previous figure we have the following equations (writing here [Li] for the LuxI concentration and [Lr] for the LuxR concentration):
- Equation for the c1 repressor block, which is the same as in the previous model:
- Equation for the c2 repressor block, which is the same as in the previous model:
- Equations for the LuxI block and LuxR block: this is the mean improvement in this model.
- Equation for parE block, this is a new equation:
- Equation for the HfsGH block (Gl in the equation), it is almost the same form as before, but there is the effect of the parE block:
- Equation for LuxR+HSL block:
There is no explicit auto regulation in this model.
In the two next sections of this study we will consider the following set of values for the different parameters which appear in our model:
Considering this set of values, we will study the dynamical behaviour of our system in absence of IPTG and when the inductor is added. We will focus mainly on the influence of the initial value of the [c2] concentration and on the quantity of IPTG which is added in the system. These two concentrations are the main values which can be controlled from an experimental point of view. In the last section of this part, we will discuss the behaviour of the system when some of the intrinsic parameters are modified.
Dynamics of the system without inductor
Steady state
In absence of an inductor (β = 0), the system is expected to stay at rest and not to produce a large amount of glue. We then consider the following initial conditions:
Due to the complexity of the system we do not expect that it will stay at rest with these initial values but that at least this state will be stable. Indeed, if we try to obtain the steady state by considering the set of algebraic equations obtained by assuming that the time derivatives of all concentrations are set to zero:
It can be seen (for the set of values which has been chosen) that the system will reach a steady state near the initial configuration (It has to be noticed that the glue (hfsGH or Gl in the model) continues to oscillate near zero then it is numerically impossible to obtain a true stationary state for which d[Gl]/dt is strictly zero.) ( [c2] increases to a very high value and the other variables go to almost zero).
Linear stability analysis
As previously explained, we assume that the configuration of the system is perturbed by :
with
and
The eigenvalues ω satisfy the characteristic equation:
In this case we obtain the expression:
This is a polynomial of 7th order and all its roots are of the form:
As we can see, the expression is much more simpler that the one we obtained in the first model. From the linear analysis point of view, the system can reach any stationary state.
 Bifurcation diagram in function of the initial value of [c2] |
 Time evolution of the glue concentration, for values of [c2](t= 0) = 10,20,30,40,50, with β = 0 |
On previous figure (right) we show the time evolution of [hfsGH] for different initial values of [c2], we see that a higher value of [c2] tends toward a diminution of the maximum value of [hfsGH] concentration. As we can see its value is very low, there is just a small increase at the beginning. During the same time the [c2] concentration reaches its steady state. Its concentration is very high; this is mainly due to the fact that the other variables as [L] or [c1] (which can repress [c2]) go to a insignificant value.
Existence of a critical value for [c2]
We have to stress that there is a minimum value for [c2] which has to be found in the system at the beginning. Indeed, as it can be seen on previous figure on the left, when the initial concentration [c2] is too small, the production of glue can start without inductor. The initial concentration [c2] has to be sufficiently high so that the c2 repressor can really represses the other blocks in the system. The critical value for [c2] depends on the parameters of the components of the system so that αi, γi and ki: [c2]c = [c2]c(αi,&gammai,ki).
Dynamics of the system with inductor
 (a):Diagram in function of the amount of inductor for different initial values of [c2] |
 (b): Evolution of the different concentrations for [c1],[c2],[hfsGH] and [L] inside the system with IPTG (model without toggle switch), for initial value of [c2]: 10,20,30,40,50 (quantities have been rescaled) |
When IPTG is added, the system leaves its initial configuration (25) and starts producing glue. Eventually, the system settles down in a steady glue producing state. On Figure (b) we show the dynamical evolution of the system in the presence of IPTG. As we can see the production of the glue begins when the repressor [c2] reaches its final value. During the decrease in [c2], [c1] symmetrically increases. When the production of the glue begins, [c1] rapidly increases to reach its steady value. In this system only a few amount of IPTG is enough to start the process (independently of the initial value of [c2]), as we see on Figure (a). We also have to notice that the concentration of [Li] and [Lr] follows the same temporal evolution as the [L] complex,this is due to the symmetry of the system. The initial value of [c2] concentration does not have any influence on the final amount of produced glue , but the initial value of [c2] has an influence on the rate of the global process: if we increase [c2], the system will take more time to begin the glue production.
In this new model, we see that the glue production can be controled by adding IPTG. But, the system is now too sensitive to the presence of the inductor. In the next section we show how this problem can be solved.
Detailed Quorum Sensing system with integrated Toggle Switch
Sensibility to the initial condition
In the previous section, we have seen that the production process can start with only a few concentration of IPTG inside the system. This characteristics can be a problem from a practical point of view. If we look at the situation from a commercial or industrial point of view it is obvious that it is mandatory to implement a control system in order to avoid a non desired production of glue. In this section we present a theoretical improvement of the actual system which can solve this problem.
An improvement proposal
A natural way to improve our system is to use the well known toggle switch system in order
to reduce the sensitiveness of our system. The toggle switch system is described in [15].
The schema for this new model is shown on Figure 19. The only difference with the model in
the Figure 16 is the presence of a negative regulation of the c2 repressor on the c1 repressor.
The equation 18 is then modified as follows:
Let focus on the two first equations of the system. Assuming that we start from an initial state
in which [L] (t = 0) = 0, the two first equations reduce to (see the corresponding diagram on
Figure 20):
This is similar to the dynamical system described in [15]. The difference is that in the first equation the initial value of ¯ appears. If we choose the parameters α1 and α2 so that the system is initially in the bistability region, we can put the system in a state for which [c2] dominates the global dynamics. In that way, we are sure that the system will stay in its initial configuration.
We can see that in this case, we have a strong connection between the initial concentration of [c2] and the concentration of IPTG we have to add in order to switch to the production mode. It means that the system will not start to produce glue if a too small amount of inductor is added inside the system (by accident for example).
Dynamics of the system without inductor
When IPTG is lacking (β= 0), we find the same behaviour as in the previous model ( figure 21), the results from the previous linear stability analysis remain valid. On the figure 21(a) we see that as previously we need to put a minimal value of [c2] to stay in the initial configuration without inductor. However, it can be seen on the figure 21(b) that we recover the same influence of the initial value of the [c2] concentration on the amplitude of the initial perturbation.
 (a): Diagram in function of the initial value of [c2] in the system with integrated toggle switch |
 (b): time evolution of the glue concentration in absence of IPTG ( integrated toggle switch model) |
Dynamics of the system with inductor
In this case we focus on the correlation which has been established between the minimal amount of IPTG which is needed to start the glue production and the initial value of the [c2] concentration. The bifurcation diagram is shown on Figure 22. If we compare with the analog diagram of the previous model ( Figure 18(a)), we see that in this new model, there is an obvious correlation between the initial value of [c2] and the quantity of IPTG which has to be added in order to produce the glue. We also notice that the needed IPTG values are higher. On Figure 23 we observe the time evolution of the different concentrations when the glue production is activated. About the glue production we have the same behaviour as for the model without toggle switch. The main difference lies in the dynamics of [c1]: it increases rapidly to a very high value. This is due to the very high initial concentration of [IPTG] which is added in this case. When [c2] varies we observe the same type of influence on the time evolution as for the case without toggle switch. We have to notice that, like in the previous model, when the glue production process starts, the final amount of glue which is produced is not influenced by the value of IPTG concentration or [c2] concentration.
 Bifurcation diagram for different initial values of [c2](from left to right: [c2] = 10; 20; 40; ; 50; 70; 100; 300; 500) |
 Time evolution of the concentrations for β = 600 for [c2] = 10, 20, 30, 40, 50 |
Discussion and conclusion
In this part, we described our biological model from a dynamical point of view. Our aim was to identify the function of the main biobricks components present in the model. We showed that the system is able to switch from a stage of zero glue production to a stage where a steady state of glue production is reached. In our last model, this transition is entirely regulated by the quantity of IPTG inductor added in the system. Lacking such an inductor we observed that the system reaches a steady state for which there is no glue. But we saw that if the initial concentration [c2] is too low, the glue production can start without inductor. In the second configuration of our model, the glue production can start for a very small quantity of IPTG. This can be a problem from a practical point of view. Because of this high sensitiveness, the glue production could start just by accident. In order to increase the robustness of our system in regard to the IPTG concentration, it is useful to improve the current system by adding a toggle switch system between the c1 repressor and the c2 repressor biobricks. Indeed, with this last improvement, there is a minimal quantity of IPTG which is needed to start the glue production, This minimal quantity is strongly correlated to the concentration of [c2] initially found in the system. Moreover an increase in the initial quantity of [c2] leads to a diminution of the disturbance amplitude in the glue concentration if there is no IPTG. In the presence of inductor, the decrease in the initial concentration [c2] leads to an increase in the growth rate of the glue production.
From a more theoretical point of view, we can also address the following question: for this global system, which genes should we choose in order to obtain an optimal equilibrium between robustness and glue production? To answer that we studied the behaviour of our theoretical model with integrated toggle switch when the parameters values (such as activation concentration, degradation rates which are intrinsic to the genes) are modified. Firstly, an increase in the [c2] degradation rate (γ2) leads to an increase in the glue production rate when IPTG is added. Secondly, when we lower the value of the kLr and kLi constant, an increase in the glue production rate is also observed. Thirdly, a decrease in the k1 constant value leads to a crash of glue production even in presence of IPTG. As we said in the previous section, the total amount of produced glue depends also on the intrinsic characteristics of the biobricks: a decrease in the γparE coefficient leads to a lower quantity of glue when the stationary state is reached, α7has also a similar influence.
