Team:Paris/Analysis/Math+Sim

From 2008.igem.org

Revision as of 08:09, 28 October 2008 by YannLC (Talk | contribs)


Mathematical analysis and simulations


Contents

FIFO

  • In order to go further into the details of mathematical analysis of the FIFO, let us remind you the set of equations used and the corresponding network (note that you can have an explanation on the way they have been built on the Model Construction page):


Subsystem1.jpg
FliA dynamics.jpg
CFP.jpg
YFP.jpg
RFP.jpg

where CFP, YFP, and RFP will be denoted below as respectively Z1,Z2 and Z3. We have implemented this system using Matlab (see the corresponding codes)









  • Then, in order to obtain simulations of the expected FIFO behaviour, we set ideal conditions, by imposing flhDC as a step function, equal to one then to zero.


FlhDC Test FIFO cresc.jpg
FlhDC Test FIFO decresc.jpg
  • We saw during the overall description of the core system that without FliA, the FIFO would presumably not work. We then simulated a first system, where [FliA] stays to zero value to confirm that qualitative conclusion :


Essai without fliA.jpg

Indeed, we may see that there is a LIFO behavior rather than the FIFO we expect...

  • Then, we simulated the entire system (that means with FliA), to check if we had
    • the lasting burst due to FliA (more important for Z3 than for Z2, and more important for Z2 than Z1) in the increasing phase.
    • the effect of fliA which maintained the concentrations to their maximum (more important for Z3 than for Z2, and more important for Z2 than Z1) in the decreasing phase. We have obtained the following curves, which are in perfect agreement with the theoretical definition of a 'FIFO' process :


Essai with fliA.jpg

  • FliA enables the curves to cross, and adds a delay on the genes that are most affected, with gives a better observability of the FIFO behavior.

Oscillations

  • In addition to the previous FIFO system, we will add a negative feedback on FlhDC promoter to observe oscillations. Thus, we obtained the following network and the corresponding equations (notice that we decide to model in a single arrowthe process through which that EnvZ repress FlhDC):
Core system.jpg


Eqn flhDC.jpg
FliA dynamics.jpg
CFP.jpg
YFP.jpg
Eqn EnvZ-RFP.jpg

If you want more details about the way we have built them, you can go to the model construction page.

If one tries to draw simulations of this system, he will obtain this kind of results :

Steady state.jpg

In a nutshell : It does not work !
Nonetheless, there are two questions we could focus on:

  • Is it possible to prove theoretically that it will not produce oscillations ? For detailled answers :
↓ mathematical analysis ↑


Understanding the dynamics

  • Theory introduction

The classical approach consists in finding the equilibrium state, by setting

Math0bis.jpg

Then, we evaluate the jacobian matrix, so as to put the system under its linearized form:

Matrix truebis.jpg


Where [EnvZ]* represents the equilibirum state of [EnvZ] and :

Def fbis.jpg

which gives :

Math1.jpg
Math2bis.jpg

Then, we want to find the eigenvalues of the jacobian matrix, because they make us understand the behavior of the system. Here is the theoretical explanation :

DE1.jpg
DE2.jpg
DE3.jpg
DE4.jpg
DE5.jpg

The λ and γ coefficients make the convergence. In our case they are negative terms, which explains the fact that we have a quick convergence for the system. Ideally, it could be convenient to find which parameter influence this coefficient, so as to play with it and be able to propose a better control of the convergence.
The μ coefficient is responsible for the oscillations. In fact, in another base, with only real matrix, this would beget cosinuses and sinuses terms. Furthermore, it is strongly linked to the period of the oscillations.

  • What about the periods
We can evaluate a value of the period we would have obtained without damping. Indeed, we have:
Period.jpg
With the parameters presented in the array before, we have evaluated the equilibrium values by solving
Math0bis.jpg

We obtained three solutions :

-0.2382 - 0.6139i

-0.2382 + 0.6139i

0.4764

We are only interested in positive real solutions. We can note that 0.4764 corresponds well to the experimental equilibrium value. Then we evaluated the eigenvalues for the jacobian :

-1.0159 + 0.5042i

-1.0159 - 0.5042i

-0.1111

We can note that the λ and γ coefficient are negative, which corroborates the fact that we obtain a convergence.

Finally for
Period.jpg
since μ=0.5042 we find a theoretical period of Ttheoretical=12.4617.
  • Understanding the attenuation

Furthermore, this theoretical study helps us understand that reducing the attenuation is tough. In fact, whatever the base in which you consider a matrix, the trace is conserved. Here, we obtain :

Trace matrix.jpg

Since λ and γ are assumed (thanks to simulations) to be negative, we cannot avoid to have their sum to be (with absolute value) superior to 2. The only chance we have would be that the weight goes in the λ parameter because this parameter plays a role in the third dimension, in which we are not interested. However, for every simulation we made, we found the opposite. Trying to evaluate with non-numerical calculus does not provide relevant material.

However, we now know that βFliA influences the attenuation. We ran a simulation with βFliA=0.02

Damped BliA.jpg

This confirms what we saw previously. In this second case, the system is therefore less attenuated, and some very tiny damped oscillations can be observed.

Last but not least, one should not forget that we chose a normalized interpretation of the system. Thus, even though the attenuation phenomenon remains unchanged, the tiny oscillations we see may indeed be observable.

  • Does this result depend on the integration methods used ? For more detailled considerations :
↓ Impact of integration methods ↑


Impact of simulation methods

This section wil just be made of some simulations to show the influence of the solver used. Indeed, when it comes to detect oscillations, some weird phenomenon can just be produced by simulation algorithms' aterfact. We have just obtained the following curves for the same initial system of ODEs.

  • Numerical simulations with Euler algorithm
Simueuler.jpg
Simueulerzoom.jpg

(The second picture is just a zoom, as defined by the red box)

  • Numerical simulations with variable-step methods
Simuode45.jpg
Simuode45zoom.jpg

(The second picture is just a zoom, as defined by the red box)

One can notice that some oscillations seem to occur in the second case, which is not the case in simulations made with constant step. Therefore, the precise tuning of the algorithm used have to be undertaken in order to avoid false detection of oscillations.


Finally, the logical continuation of the process will be to try to improve the biological system in order to find oscillations. This will be described in the part below, as well as in the section System Improvements.

System analysis

Simulations of the core system displayed above reveal that it does not exhibit an oscillatory behavior. In this section we use the model of the core system to try to figure out the contribution of some key characteristics of the network topology on the dynamics of the system. This analysis is done by successively simulating altered forms of the system.

The following table lists the different forms of the system considered and their respective outcome :

Description Simulation Evaluation Interpretation
Core system0.png
Core system (without modifications).
Core0.png
Orgsad.gif Non oscillating system.
Core system1.png
Removing the activation of envZ via FlhDC.
Core1.png
Orgstar.gif Non oscillating system. Activation of envZ via FlhDC has a negative effect on the presence of oscillations. This interaction creates a 'shortcut' in the negative loop, diminishing the delay between FlhDC and EnvZ. Notice that without this shortcut EnvZ increases after Flia. This interaction is already low in the core system, its deletion thus has a moderate positive effect.
Core system2.png
Removing the auto-activation of flia
Core2.png
Orgstar.gifOrgstar.gif Damped oscillations. Auto-activation of Flia slows down its dynamics. Without this auto-activation the slope of Flia is similar to the slope of FlhDC. For oscillations to appear, different species need to have comparable time scales : oscillations are much more easier without Flia auto-activation.
Core system3.png
Enhancing the inhibition of FlhDC via envZ
Core5.png
Orgstar.gif Non oscillating system. A stronger inhibition slightly improves the behavior of the system.
Core system4.png
Combining three modifications all together.
Core4.png
Orgstar.gifOrgstar.gifOrgstar.gif Damped oscillations with greater amplitude. Combining the three previous variants provides the best possible behavior.


The core system cannot provide oscillations. In particular, auto-activation of flia and activation of envZ by FlhDC both have serious negative effects on oscillations. But these actions are required for the FIFO behavior of the system.

In the next section, we examine significant topology improvements of the core system that are potentially better for oscillations and that preserve FIFO required interactions.

Navigator

Back to the overall presentation of our system and model
Top of the page
Mathematical analysis Effects of integration methods