silico.biotoul.fr
 

M1 MABS BBS Math TD Modelisation

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m
m (Prise en main de la librairie deSolve)
Line 6: Line 6:
[[Image:autoregulation.png]]
[[Image:autoregulation.png]]
 +
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 du delta entre la quantité produite et la quantité dégradée. Pour faire simple, on prendra comme paramètres :
 +
* <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>
 +
 +
On définit les paramètres 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 :
 +
<source lang='rsplus'>
 +
autoregulation = function(t, 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'>
 +
times=1:200 # durée de la simulation
 +
out <- ode(y = state, times = times, func = autoregulation, parms = parameters)
 +
pos=out
 +
head(out)
 +
plot(out, xlab = "time", ylab = "X")
 +
 +
plot(times, out[,'X']/max(out[,'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]]
= Système discret : plantes annuelles =
= Système discret : plantes annuelles =

Revision as of 14:19, 2 October 2014

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

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 du delta entre la quantité produite et la quantité dégradée. Pour faire simple, on prendra comme paramètres :

  • 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

On définit les paramètres de la manière suivante :

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

Et les équations sous forme de fonction R :

autoregulation = function(t, 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 :

times=1:200 # durée de la simulation
out <- ode(y = state, times = times, func = autoregulation, parms = parameters)
pos=out
head(out)
plot(out, xlab = "time", ylab = "X")
 
plot(times, out[,'X']/max(out[,'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) :

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 plantes et d'un sol vierge
  6. Au bout de combien d'années n'a-t-on plus de plante en partant de 50 plants ?
  7. Déterminer le taux γ minimal nécessaire au maintient des plantes pour les paramètres précédents (sauf γ)

Exemple d'utilisation de deSolve avec l'équation de Lorenz décrivant un comportement chaotique (1963) tirée de la vignette :

library(deSolve)
 
# Codage des paramètres du système (A, B et C) sous forme de liste 
parameters <- c(a = -8/3, b = -10, c = 28)
# Codage de l'état du système (X, Y et Z) également sous forme de liste
state <- c(X = 1, Y = 1, Z = 1)
 
# La fonction appliquant les équations d'état du système (calcul de dX, dY et dZ)
Lorenz<-function(t, state, parameters) {
	with(as.list(c(state, parameters)),{
		# rate of change
		dX <- a*X + Y*Z
		dY <- b * (Y-Z)
		dZ <- -X*Y + c*Y - Z
		# return the rate of change
		list(c(dX, dY, dZ))
	})
	# end with(as.list ...
}
 
# durée de la simulation (doit commencer par 0)
times <- seq(0, 100, by = 0.01)
# résolution du système et simulation
out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
 
head(out)
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")
plot(out[, "X"], out[, "Z"], pch = ".")
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)

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

Ceci se traduit par l'équation 1 de la publication.

L'abondance totale de KaiB et KaiC dépend de leur taux de transcription et de la vitesse de dégradation des ARNm. Les auteurs en déduisent l'équation 2.

Dans un premier temps, implémenter le modèle simple où l'abondance de KaiA active est constante. Vous utiliserez les paramètres spécifiés dans l'article.

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