silico.biotoul.fr
 

M1 MABS BBS Math TD Modelisation

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Travail à réaliser)
m (Prise en main de la librairie deSolve)
(33 intermediate revisions not shown)
Line 1: Line 1:
-
= Système discret : plantes annuelles =
+
= Prise en main de la librairie deSolve =
 +
Nous allons commencer par modéliser des boucles de régulations de l'expression de gènes à partir d'exemples théoriques simplifiés à l'extrême.
 +
 
 +
Tout d'abord, vous avez ci-contre le schéma d'une auto-régulation positive ou négative avec la concentration (X/Xmax) ou quantité de produit de X en fonction du temps.
 +
 
 +
[[Image:autoregulation.png|frame|Boucles d'auto-régulation (positive ou négative). source: Shoval and Alon, Cell 2010]]
 +
 
 +
Pour le modèle d'autorégulation positive, nous aurons comme variable d'état du système <tt>X</tt> : la concentration de l'espèce moléculaire X, et donc comme début de script :
 +
<source lang='rsplus'>
 +
library(deSolve) # chargement de la librairie
 +
 
 +
state = c( X=0.1 )
 +
</source>
 +
 
 +
Ensuite la variation de concentration de X dépend de la différence entre la quantité produite et la quantité dégradée. Pour faire simple, les paramètres, sous forme de vecteur dont les éléments sont nommés, sont les suivants :
 +
* <tt>v</tt> la vitesse de synthèse
 +
* <tt>auto</tt> une constante d'auto-régulation
 +
* <tt>degradation</tt> la vitesse de dégradation
 +
 
 +
Ainsi, nous obtenons l'équation suivante :<br/>
 +
<math>\frac{dx}{dt} = v \times auto \times X - degradation \times X</math>
 +
 
 +
Les paramètres, là aussi sous forme de vecteur, sont définis de la manière suivante :
 +
<source lang='rsplus'>
 +
parameters = c(
 +
v = 0.04, degradation = 0.02, auto = 0.8
 +
)
 +
</source>
 +
 
 +
Et les équations sous forme de fonction R avec comme paramètres <tt>timelapse</tt> : le temps ; <tt>state</tt> : l'état du système ; et <tt>parameters</tt> : les paramètres du système
 +
<source lang='rsplus'>
 +
autoregulation = function(timelapse, state, parameters) {
 +
with(as.list(c(state, parameters)),{
 +
dX = v * auto * X - degradation * X
 +
##
 +
list(c(dX)) # liste renvoyée par la fonction
 +
})
 +
}
 +
</source>
 +
 
 +
Il n'y a plus qu'à effectuer une simulation du système et observer les résultats :
 +
<source lang='rsplus'>
 +
timelapse=1:200 # durée de la simulation
 +
sim1 <- ode(y = state, times = timelapse, func = autoregulation, parms = parameters)
 +
head(sim1)
 +
tail(sim1)
 +
plot(sim1, xlab = "times", ylab = "X")
 +
 
 +
plot(timelapse, sim1[,'X']/max(sim1[,'X']), type='l', xlab='time', ylab='X/Xmax') # plot avec en ordonnées X/Xmax
 +
</source>
 +
 
 +
Vous devriez obtenir le résultat suivant :<br/>
 +
[[Image:autoregulation_positive.png|400px]]
 +
 
 +
Effectuons à présent le même travail pour une boucle d'auto-régulation négative. Pour cela, considérons <tt>v</tt> comme la vitesse maximale de synthèse ; cette vitesse sera d'autant plus facilement atteinte que la concentration de X est faible, d'où :<br/>
 +
<math>\frac{dx}{dt} = v \times \frac{auto}{auto + X} - degradation \times X</math>
 +
 
 +
 
 +
Modifier la fonction préalable et effectuer la simulation. Vous devriez obtenir le graphique suivant (avec en bleu la courbe précédente et en noir la courbe correspondant à une régulation négative) :<br/>
 +
[[Image:autoregulation_pos_neg.png|400px]]
 +
 
 +
== Autre exemple : boucle de régulation négative ==
 +
Nous allons nous intéresser maintenant au circuit ci-contre (circuit de gauche) :<br/>
 +
[[Image:Negative-feedback loop.png|frame|Boucle de régulation négative. source : Shoval and Alon, Cell 2010]]
 +
 
 +
Dans ce circuit, X s'auto-régule positivement et active également Y, mais plus lentement (que lui-même X). Y, quant à elle, réprime X, de manière rapide, par exemple en stimulant la dégradation de X.
 +
 
 +
On pourra obtenir les équations (ultra simplifiées) suivantes :<br/>
 +
<math>
 +
\frac{dx}{dt} = v \times X - degradation \times (Y \times constante) \times X
 +
</math><br/>
 +
<math>
 +
\frac{dy}{dt} = v/2 \times X - degradation \times Y
 +
</math><br/>
 +
avec comme paramètres :
 +
* <tt>v</tt>: la vitesse de synthèse (et v/2 pour une activation lente de Y)
 +
* <tt>degradation</tt>: la vitesse de dégradation
 +
 
 +
On obtient donc le script suivant pour le scénario "damped oscillations" :
 +
<source lang='rsplus'>
 +
timelapse=1:2000
 +
parameters = c( v = 1, degradation = 0.01 )
 +
 
 +
state = c(X=0.10, Y = 0.10)
 +
 
 +
nfl = function(t, state, parameters) { with(as.list(c(state, parameters)),{
 +
dX = v * X - degradation*(Y*10) * X
 +
dY = v/2 * X - degradation * Y
 +
##
 +
list(c(dX, dY)) }) }
 +
 
 +
sim2 <- ode(y = state, times = timelapse, func = nfl, parms = parameters); X = sim2[,'X'] ; Y = sim2[,'Y']
 +
plot(timelapse, X, t='l', ylim= c(min(X,Y), max(X,Y)), col='red', ylab='')
 +
lines(timelapse, Y, col='blue')
 +
</source>
 +
qui produit le résultat suivant :
 +
[[Image:Damped oscillations.png|400px]]
 +
 
 +
= Système continu : Modèle de rythme circadien de cyanobactérie =
 +
Travail basé sur la publication [http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1557573/ A Model for the Circadian Rhythm of Cyanobacteria that Maintains Oscillation without Gene Expression] de Kurosawa ''et al.''.
 +
 
 +
Dans ce travail, ils proposent un modèle de régulation de la phosphorylation de la protéine KaiC capable de générer des oscillations circadiennes chez les cyanobactéries. Nous allons ici implémenter leur système et vérifier certains de leurs résultats.
 +
 
 +
L'état du système est le suivant :
 +
* la protéine KaiA se trouve sous 2 formes : KaiA (inactive) et Kai<sup>*</sup> (active)
 +
* la protéine KaiC également : KaiC et KaiC-P (phosphorylée)
 +
* la protéine KaiB aussi : KaiB (inactive) et KaiB<sup>*</sup> (active)
 +
 
 +
Le système évolue :
 +
* KaiA<sup>*</sup> stimule la phosphorylation de KaiC
 +
* KaiB<sup>*</sup> stimule la déphosphorylation de KaiC
 +
 
 +
* KaiC-P stimule l'inactivation de KaiB. L'activité de KaiB dépend de sa localisation et éventuellement de plusieurs formes chimiques de KaiB (comme pour KaiC).
 +
 
 +
Cela correspond à la figure 3C de la publication :<br/>
 +
[[Image:KaiABC schema.png]]
 +
 
 +
 
 +
Ceci se traduit par les équations 1,2 et 3 de la publication.<br/>
 +
Pour ''x'', la quantité de de KaiC-P, avec :
 +
* p et a : la vitesse de phosphorylation de KaiC
 +
* C<sub>0</sub> et s: la concentration initiale de KaiC et sa variation de concentration, donc C<sub>0</sub>s - x la quantité de KaiC non phosphorylée
 +
* b (y+f) : la vitesse de déphosphorylation de KaiC-P par KaiB<sup>*</sup>
 +
<math>
 +
\frac{dx}{dt} = pa (C_0s - x) - bx (y + f)
 +
</math>
 +
 
 +
Pour ''y'', la quantité de KaiB<sup>*</sup> (KaiB active), avec :
 +
* g : vitesse d'activation de KaiB
 +
* <math>k_1x^n / (q^n+x^n) </math> : la désactivation de KaiB<sup>*</sup> par kaiC-P
 +
<math>\frac{dy}{dt} = g(B_0s-y) - k_1 \frac{x^n}{q^n+x^n}y</math>
 +
 
 +
Les variations de concentration de KaiB totale et KaiC totale sont données dans l'équation 2, avec :
 +
* <math>\epsilon_1</math>: constante permettant d'ajuster les vitesse de synthèse et de dégradation
 +
* <math>\lambda</math>: la vitesse de transcription
 +
* <math>\mu</math>: la vitesse de décomposition
 +
* la transcription de kaiBC est inhibée par la concentration en KaiC-P, d'où la vitesse de transcription: <math>\frac{\lambda}{1+h_1x^m}</math>
 +
<math>\frac{dB}{dt}= \epsilon_1 (\frac{B_0 \lambda}{1+h_1x^m} - \mu B)</math>
 +
 
 +
<math>\frac{dC}{dt}= \epsilon_1 (\frac{C_0 \lambda}{1+h_1x^m} - \mu C)</math>
 +
 
 +
En synthétisant ces différentes équations et en remplaçant dans l'équation 1 C<sub>0</sub>s par C et B<sub>0</sub>s par B, on obtient le système suivant :<br/>
 +
<math>
 +
\frac{dx}{dt} = pa(C-x) - bx(y+f)
 +
</math>
 +
 
 +
<math>
 +
\frac{dy}{dt} = g(B-y) - k_1yx^n/(q^n+x^n)
 +
</math>
 +
 
 +
<math>\frac{dB}{dt}= \epsilon_1 (B_0 \lambda / (1+h_1x^m) - \mu B)</math>
 +
 
 +
<math>\frac{dC}{dt}= \epsilon_1 (C_0 \lambda / (1+h_1x^m) - \mu C)</math>
 +
 
 +
 
 +
Les valeurs des paramètres et concentration initiale pour obtenir la figure 4B sont fournis dans la légende des figures 2A et 4A, et 4B pour epsilon1.
 +
 
 +
Vous devriez obtenir les courbes suivantes :
 +
[[Image:LL-CBXY.png]] et [[Image:LL-kaiC-kaiB.png]]
 +
 
 +
= Exercice pour chez soi. Système discret : plantes annuelles =
Les [http://fr.wikipedia.org/wiki/Plante_annuelle plantes annuelles] tirent leur nom du fait qu'elles vivent moins d'un an. Elles effectuent leur cycle de vie en une seule année et passent l'hiver uniquement sous forme de graines. Nous allons nous intéresser à l'élaboration d'un modèle pour l'étude du développement et de la prolifération de ces plantes.
Les [http://fr.wikipedia.org/wiki/Plante_annuelle plantes annuelles] tirent leur nom du fait qu'elles vivent moins d'un an. Elles effectuent leur cycle de vie en une seule année et passent l'hiver uniquement sous forme de graines. Nous allons nous intéresser à l'élaboration d'un modèle pour l'étude du développement et de la prolifération de ces plantes.
Line 20: Line 180:
# Identifier sur ce schéma les variables d'états qui vont être utilisées
# Identifier sur ce schéma les variables d'états qui vont être utilisées
# Écrire les relations liant les variables d'une année <i>n</i> à l'année suivante <i>n + 1</i>.
# Écrire les relations liant les variables d'une année <i>n</i> à l'année suivante <i>n + 1</i>.
-
# Implémenter le modèle obtenu dans [http://www.r-project.org R] avec la librairie [http://cran.cict.fr/web/packages/deSolve/index.html deSolve]. Pour cela, on utilisera les valeurs suivantes pour les paramètres :
+
# Implémenter le modèle obtenu dans [http://www.r-project.org R] avec la librairie [http://cran.r-project.org/web/packages/deSolve/index.html deSolve]. Pour cela, on utilisera les valeurs suivantes pour les paramètres :
#* &gamma; = 2
#* &gamma; = 2
#* &alpha; = 0.5
#* &alpha; = 0.5
#* &beta; = 0.25
#* &beta; = 0.25
#* &sigma; = 0.8
#* &sigma; = 0.8
-
# Effectuer une simulation, sur 20 ans à partir de 100 plantes et d'un sol vierge.
+
# Effectuer une simulation, sur 20 ans à partir de 100 graines et d'un sol vierge
 +
# Au bout de combien d'années n'a-t-on plus de plante en partant de 50 graines ?
 +
# Déterminer le taux &gamma; minimal nécessaire au maintient des plantes pour les paramètres précédents (sauf &gamma;)
 +
'''Remarque :''' pour les simulations, vous utiliserez comme précédemment la fonction ode(...) mais il faudra bien préciser le paramètre <tt>method</tt> dans le cadre de données discrètes (on "saute" de l'année ''n'' à l'année ''n+1'' sans intermédiaire), donc : <tt>out = ode( ..., method='euler')</tt>
 +
<!--
'''Exemple d'utilisation de deSolve avec l'[http://bcev.nfrance.com/Lorenz/equations.htm équation de Lorenz] décrivant un comportement chaotique (1963) tirée de la vignette :'''
'''Exemple d'utilisation de deSolve avec l'[http://bcev.nfrance.com/Lorenz/equations.htm équation de Lorenz] décrivant un comportement chaotique (1963) tirée de la vignette :'''
Line 61: Line 225:
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)
</source>
</source>
 +
-->
 +
 +
<!--
 +
Dans la suite, les auteurs élaborent un modèle où les oscillations continuent même en l'absence de transcription, en supposant que l'abondance des protéines est constante. Ils supposent notamment que l'abondance de KaiA est constante et que la forme active oscille. Ils en déduisent l'équation 4.
 +
 +
'''Implémenter ce nouveau modèle avec les nouveaux paramètres indiqués.'''
 +
 +
Vous devriez obtenir quelque chose de la forme suivante :
 +
[[Image:LLDD-CBXYA.png]]
= Système continu : vinification =
= Système continu : vinification =
 +
Inspiré de la thèse d'Isabelle Queinnec
 +
 +
Pendant le processus de vinification, la fermentation alcoolique par les levures consiste à dégrader du glucose en anaérobie ce qui tend à favoriser la production d'alcool par rapport à une fermentation aérobie qui favoriserait la croissance de la biomasse. Dans certains cas, une faible aération facilite la synthèse de la biomasse et permet une augmentation de la productivité en alcool. La réaction théorique de Gay-Lussac est la suivante :
 +
 +
C6H12O6 -> C5H5OH + 2CO2 '''''(1)'''''<!--math>C_6H_{12}O_6 \leftarrow 2C_5h_5OH + 2CO_2</math-->
 +
<!--
 +
Son rendement théorique est de 0.511g d'éthanol produit par gramme de glucose consommé. Dans la pratique, ce rendement atteint 80% de cette valeur.
 +
 +
L'accroissement de la population de levures est proportionnel à la population elle-même, c'est-à-dire que la vitesse de croissance cellulaire r<sub>x</sub> est proportionnelle à la concentration en microorganismes X à chaque instant :
 +
r<sub>x</sub> = 1/V d(XV)/dt = µ X (en g/lh) '''''(2)'''''
 +
 +
On désigne par µ le taux de croissance, V<sub>S</sub> le taux de dégradation de glucose, et V<sub>P</sub> le taux de production d'alcool. Ces paramètres représentent les vitesses spécifiques de la réaction relatives à la biomasse. µ, V<sub>S</sub> et V<sub>P</sub> s'expriment selon les équations suivantes :
 +
µ = r<sub>x</sub>/X
 +
V<sub>S</sub> = r<sub>S</sub>/X '''''(3)'''''
 +
V<sub>P</sub> = r<sub>P</sub>/X
 +
 +
La croissance microbienne dépend d'un grand nombre de variables d'environnement comme la température, le pH, l'agitation, l'oxygénation, la dégradation du substrat, et l'accumulation de l'alcool. La difficulté de prendre en compte tous ces facteurs dans un modèle de croissance fait que l'approche usuelle consiste à fixer les variables d'environnement à des valeurs optimales déterminées expérimentalement et intervenant dans la valeur numérique des constantes. De ce fait, l'expression du taux de croissance ne représente plus que l'influence des concentrations de substrat, de produit et de biomasse.
 +
 +
L'influence de la teneur en substrat sur la croissance se traduit par une limitation lors d'un manque de substrat, et par une inhibition lors d'un excès de substrat. Le modèle d'Andrews propose une extension du celui de Monod permettant de prendre en compte l'effet limitant de l'excès de substrat :
 +
 +
µ(S) = µ<sub>max</sub> S / (K<sub>s</sub> + S + S<sup>2</sup>/K<sub>i</sub>) '''''(4)'''''
 +
 +
avec K<sub>s</sub> la constante de Monod et K<sub>i</sub> la constante d'inhibition.
 +
 +
Le substrat consommé est utilisé pour la croissance cellulaire et pour le maintien de la biomasse dans un état viable, ce qui produit de l'alcool. On exprime le taux de dégradation par :
 +
V<sub>S</sub> = µ/Y<sub>X</sub> + V<sub>P</sub>/Y<sub>P</sub> + m '''''(5)'''''
 +
avec
 +
* Y<sub>X</sub> le rendement de conversion de substrat en biomasse
 +
* Y<sub>P</sub> le rendement de conversion de substrat en produit
 +
* m le coefficient de maintenance
 +
Le coefficient de maintenance est généralement négligé devant les deux premiers termes de l'expression de V<sub>S</sub>. De plus, certains modèles expriment le taux de dégradation proportionnellement au seul taux de production.
 +
 +
Compte tenu des observations précédentes, pour simuler le système fermentaire, on retient un modèle dérivé de celui de Ghose-Tyagi, dans lequel l'inhibition par l'éthanol n'influence que la croissance des microorganismes :
 +
µ = µ<sub>max</sub>S/(K<sub>S</sub> + S + S<sup>2</sup>/K<sub>i</sub>)(1-P/P<sub>m</sub>)
 +
V<sub>P</sub> = V<sub>max</sub>S/(K<sub>S</sub> + S +S<sup>2</sup>/K<sub>i</sub>) '''''(6)'''''
 +
V<sub>S</sub> = V<sub>P</sub>/Y<sub>P</sub>
 +
avec µ<sub>max</sub> le taux de croissance optimal, V<sub>max</sub> le taux optimal de production d'alcool et P<sub>m</sub> le coefficient d'inhibition par l'alcool.
 +
 +
La première étape de la modélisation mathématique de la fermentation alcoolique consiste à déterminer les vitesses de réactions des variables macroscopiques de système, à savoir de la biomasse, du substrat et de l'alcool produit.
 +
 +
La deuxième étape permet de déterminer les équations d'état du système. Les variables d'état sont les concentrations en microorganisme X, en substrat S, en alcool P, et le volume réactionnel V. Les équations sont établies à partir de l'écriture des bilans matières relatifs aux composés de la réaction, c'est-à-dire en écrivant que la variation des quantités correspondantes pendant un intervalle de temps ''dt'' est égale à la somme de ce qui est apporté et créé, diminuée de ce qui est dégradé ou soutiré. Ceci nous permet d'écrire :
 +
d(VX) = r<sub>X</sub>V dt = µXV dt
 +
d(VS) = -r<sub>S</sub>V dt = -v<sub>S</sub>XV dt '''''(7)'''''
 +
d(VP) = r<sub>P</sub>V dt = v<sub>P</sub>XV dt
 +
 +
 +
# Implémenter le modèle à l'aide la bibliothèque deSolve
 +
# Au bout de combien de temps le jus dépasse 12° avec les paramètres et variables d'état suivants :
 +
#* K<sub>S</sub> = 5 g/l
 +
#* K<sub>i</sub> = 201 g/l
 +
#* µ<sub>max</sub> = 0.54 /h
 +
#* V<sub>max</sub> = 2.1 /h
 +
#* Y<sub>P</sub> = 0.43
 +
#* P<sub>m</sub> = 15
 +
#* X = 3 g/l
 +
#* S = 38 g/l
 +
#* P = 0 g/l
 +
#* V = 480 l
 +
 +
-->

Revision as of 08:45, 12 October 2018

Contents

Prise en main de la librairie deSolve

Nous allons commencer par modéliser des boucles de régulations de l'expression de gènes à partir d'exemples théoriques simplifiés à l'extrême.

Tout d'abord, vous avez ci-contre le schéma d'une auto-régulation positive ou négative avec la concentration (X/Xmax) ou quantité de produit de X en fonction du temps.

Boucles d'auto-régulation (positive ou négative). source: Shoval and Alon, Cell 2010

Pour le modèle d'autorégulation positive, nous aurons comme variable d'état du système X : la concentration de l'espèce moléculaire X, et donc comme début de script :

library(deSolve) # chargement de la librairie
 
state = c( X=0.1 )

Ensuite la variation de concentration de X dépend de la différence entre la quantité produite et la quantité dégradée. Pour faire simple, les paramètres, sous forme de vecteur dont les éléments sont nommés, sont les suivants :

  • v la vitesse de synthèse
  • auto une constante d'auto-régulation
  • degradation la vitesse de dégradation

Ainsi, nous obtenons l'équation suivante :
\frac{dx}{dt} = v \times auto \times X - degradation \times X

Les paramètres, là aussi sous forme de vecteur, sont définis de la manière suivante :

parameters = c( 
	v = 0.04, degradation = 0.02, auto = 0.8
)

Et les équations sous forme de fonction R avec comme paramètres timelapse : le temps ; state : l'état du système ; et parameters : les paramètres du système

autoregulation = function(timelapse, state, parameters) {
	with(as.list(c(state, parameters)),{
		dX = v * auto * X - degradation * X
		##
		list(c(dX)) # liste renvoyée par la fonction
	})
}

Il n'y a plus qu'à effectuer une simulation du système et observer les résultats :

timelapse=1:200 # durée de la simulation
sim1 <- ode(y = state, times = timelapse, func = autoregulation, parms = parameters)
head(sim1)
tail(sim1)
plot(sim1, xlab = "times", ylab = "X")
 
plot(timelapse, sim1[,'X']/max(sim1[,'X']), type='l', xlab='time', ylab='X/Xmax') # plot avec en ordonnées X/Xmax

Vous devriez obtenir le résultat suivant :

Effectuons à présent le même travail pour une boucle d'auto-régulation négative. Pour cela, considérons v comme la vitesse maximale de synthèse ; cette vitesse sera d'autant plus facilement atteinte que la concentration de X est faible, d'où :
\frac{dx}{dt} = v \times \frac{auto}{auto + X} - degradation \times X


Modifier la fonction préalable et effectuer la simulation. Vous devriez obtenir le graphique suivant (avec en bleu la courbe précédente et en noir la courbe correspondant à une régulation négative) :

Autre exemple : boucle de régulation négative

Nous allons nous intéresser maintenant au circuit ci-contre (circuit de gauche) :

Boucle de régulation négative. source : Shoval and Alon, Cell 2010

Dans ce circuit, X s'auto-régule positivement et active également Y, mais plus lentement (que lui-même X). Y, quant à elle, réprime X, de manière rapide, par exemple en stimulant la dégradation de X.

On pourra obtenir les équations (ultra simplifiées) suivantes :

\frac{dx}{dt} = v \times X - degradation \times (Y \times constante) \times X

\frac{dy}{dt} = v/2 \times X - degradation \times Y
avec comme paramètres :

  • v: la vitesse de synthèse (et v/2 pour une activation lente de Y)
  • degradation: la vitesse de dégradation

On obtient donc le script suivant pour le scénario "damped oscillations" :

timelapse=1:2000
parameters = c( v = 1, degradation = 0.01 )
 
state = c(X=0.10, Y = 0.10)
 
nfl = function(t, state, parameters) { 	with(as.list(c(state, parameters)),{
		dX = v * X - degradation*(Y*10) * X
		dY = v/2 * X - degradation * Y
		##
		list(c(dX, dY)) }) }
 
sim2 <- ode(y = state, times = timelapse, func = nfl, parms = parameters); X = sim2[,'X'] ; Y = sim2[,'Y']
plot(timelapse, X, t='l', ylim= c(min(X,Y), max(X,Y)), col='red', ylab='')
lines(timelapse, Y, col='blue')

qui produit le résultat suivant :

Système continu : Modèle de rythme circadien de cyanobactérie

Travail basé sur la publication A Model for the Circadian Rhythm of Cyanobacteria that Maintains Oscillation without Gene Expression de Kurosawa et al..

Dans ce travail, ils proposent un modèle de régulation de la phosphorylation de la protéine KaiC capable de générer des oscillations circadiennes chez les cyanobactéries. Nous allons ici implémenter leur système et vérifier certains de leurs résultats.

L'état du système est le suivant :

  • la protéine KaiA se trouve sous 2 formes : KaiA (inactive) et Kai* (active)
  • la protéine KaiC également : KaiC et KaiC-P (phosphorylée)
  • la protéine KaiB aussi : KaiB (inactive) et KaiB* (active)

Le système évolue :

  • KaiA* stimule la phosphorylation de KaiC
  • KaiB* stimule la déphosphorylation de KaiC
  • KaiC-P stimule l'inactivation de KaiB. L'activité de KaiB dépend de sa localisation et éventuellement de plusieurs formes chimiques de KaiB (comme pour KaiC).

Cela correspond à la figure 3C de la publication :
Image:KaiABC schema.png


Ceci se traduit par les équations 1,2 et 3 de la publication.
Pour x, la quantité de de KaiC-P, avec :

  • p et a : la vitesse de phosphorylation de KaiC
  • C0 et s: la concentration initiale de KaiC et sa variation de concentration, donc C0s - x la quantité de KaiC non phosphorylée
  • b (y+f) : la vitesse de déphosphorylation de KaiC-P par KaiB*


\frac{dx}{dt} = pa (C_0s - x) - bx (y + f)

Pour y, la quantité de KaiB* (KaiB active), avec :

  • g : vitesse d'activation de KaiB
  • k1xn / (qn + xn) : la désactivation de KaiB* par kaiC-P

\frac{dy}{dt} = g(B_0s-y) - k_1 \frac{x^n}{q^n+x^n}y

Les variations de concentration de KaiB totale et KaiC totale sont données dans l'équation 2, avec :

  • ε1: constante permettant d'ajuster les vitesse de synthèse et de dégradation
  • λ: la vitesse de transcription
  • μ: la vitesse de décomposition
  • la transcription de kaiBC est inhibée par la concentration en KaiC-P, d'où la vitesse de transcription: \frac{\lambda}{1+h_1x^m}

\frac{dB}{dt}= \epsilon_1 (\frac{B_0 \lambda}{1+h_1x^m} - \mu B)

\frac{dC}{dt}= \epsilon_1 (\frac{C_0 \lambda}{1+h_1x^m} - \mu C)

En synthétisant ces différentes équations et en remplaçant dans l'équation 1 C0s par C et B0s par B, on obtient le système suivant :

\frac{dx}{dt} = pa(C-x) - bx(y+f)


\frac{dy}{dt} = g(B-y) - k_1yx^n/(q^n+x^n)

\frac{dB}{dt}= \epsilon_1 (B_0 \lambda / (1+h_1x^m) - \mu B)

\frac{dC}{dt}= \epsilon_1 (C_0 \lambda / (1+h_1x^m) - \mu C)


Les valeurs des paramètres et concentration initiale pour obtenir la figure 4B sont fournis dans la légende des figures 2A et 4A, et 4B pour epsilon1.

Vous devriez obtenir les courbes suivantes : Image:LL-CBXY.png et Image:LL-kaiC-kaiB.png

Exercice pour chez soi. Système discret : plantes annuelles

Les plantes annuelles tirent leur nom du fait qu'elles vivent moins d'un an. Elles effectuent leur cycle de vie en une seule année et passent l'hiver uniquement sous forme de graines. Nous allons nous intéresser à l'élaboration d'un modèle pour l'étude du développement et de la prolifération de ces plantes.

Description du cycle

Les graines germent en mai pour donner des plantes. Ces plantes fleurissent et produisent des graines qui sont dispersées à partir d'août. Une partie des graines sont infectées par divers pathogènes et meurent ou bien sont consommées par différents animaux pendant l'hiver. L'année suivante, le cycle recommence : une partie des graines germent, ... Les graines n'ayant pas germé attendront l'année suivante, si elles ne germent pas l'année suivante, elles ne germeront plus.

Paramètres du modèle

  • γ : nombre de graines produites par une plante
  • α : proportion de graines produites l'année précédente qui germent en mai
  • β : proportion de graines d'il y a 2 ans qui germent en mai
  • σ : proportion de graines qui survivent à l'hiver

Etat du système

En début d'année, l'état du système est donné par le nombre de graines de l'année précédente et d'il y a 2 ans.


Travail à réaliser

  1. Faire un schéma du système/cycle
  2. Identifier sur ce schéma les variables d'états qui vont être utilisées
  3. Écrire les relations liant les variables d'une année n à l'année suivante n + 1.
  4. Implémenter le modèle obtenu dans R avec la librairie deSolve. Pour cela, on utilisera les valeurs suivantes pour les paramètres :
    • γ = 2
    • α = 0.5
    • β = 0.25
    • σ = 0.8
  5. Effectuer une simulation, sur 20 ans à partir de 100 graines et d'un sol vierge
  6. Au bout de combien d'années n'a-t-on plus de plante en partant de 50 graines ?
  7. Déterminer le taux γ minimal nécessaire au maintient des plantes pour les paramètres précédents (sauf γ)

Remarque : pour les simulations, vous utiliserez comme précédemment la fonction ode(...) mais il faudra bien préciser le paramètre method dans le cadre de données discrètes (on "saute" de l'année n à l'année n+1 sans intermédiaire), donc : out = ode( ..., method='euler')