1.1 Boucle d’auto-régulation positive

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.

Fig. 1 : boucles d’auto-régulation. 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.autopos = c( X=0.1 ) # variable(s) d'état du système

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 qui correspondent à des constantes, là aussi sous forme de vecteur, sont définis de la manière suivante :

params.autopos = c( 
  v = 0.04, degradation = 0.02, auto = 0.8
)

Et les équations sont écrites 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.

Ecriture de l’équation sous forme de fonction R

func.autopos = 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.autopos, times = timelapse, func = func.autopos, parms = params.autopos)
head(sim1) ; tail(sim1)

##      time         X
## [1,]    1 0.1000000
## [2,]    2 0.1012072
## [3,]    3 0.1024291
## [4,]    4 0.1036656
## [5,]    5 0.1049171
## [6,]    6 0.1061837

##        time        X
## [195,]  195 1.025741
## [196,]  196 1.038124
## [197,]  197 1.050657
## [198,]  198 1.063340
## [199,]  199 1.076177
## [200,]  200 1.089169

plot(sim1, xlab = "time", ylab = "X", main='Auto-régulation positive [X]')

plot(sim1[,'time'], sim1[,'X']/max(sim1[,'X']), type='l', xlab='time', ylab='X/Xmax', main='Auto-régulation positive [X]/[Xmax]') # plot avec en ordonnées X/Xmax

state.autoneg1 = c( X=0.1 ) 
params.autoneg1 = c( v = 0.04, degradation = 0.02, auto = 0.8 )
autoregulation.negative = function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
  dX = v * (auto/ (auto + X)) - degradation * X  # autoregulation negative
  list(c(dX))
})}

sim2 <- ode(y = state.autoneg1, times = timelapse, func = autoregulation.negative, parms = params.autoneg1)

head(sim2); tail(sim2)

##      time         X
## [1,]    1 0.1000000
## [2,]    2 0.1325927
## [3,]    3 0.1633825
## [4,]    4 0.1925396
## [5,]    5 0.2202054
## [6,]    6 0.2465019

##        time         X
## [195,]  195 0.9249078
## [196,]  196 0.9249606
## [197,]  197 0.9250117
## [198,]  198 0.9250613
## [199,]  199 0.9251094
## [200,]  200 0.9251561

plot(sim2, xlab = "time", ylab = "-", ylog=TRUE, main='Auto-régulation négative')

plot(sim2[,'time'], sim2[,'X']/max(sim2[,'X']), type='l', ylim=c(0,1), xlab='time', ylab='X/Xmax', main='Auto-régulation positive (bleue) et négative (noir)')
lines(sim1[,'time'], sim1[,'X']/max(sim1[,'X']), type='l', col='blue')

1.1.2 Auto-régulation négative, exemple à 2 composants

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

Fig. 2 : boucles d’auto-régulations négatives. 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} = \frac{1}{2} \times v \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 code suivant pour le scénario damped oscillations :

timelapse=1:2000
params.nfl = c( v = 1, degradation = 0.01 )

state.nfl = 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)) }) 
}

sim.nfl <- ode(y = state.nfl, times = timelapse, func = nfl, parms = params.nfl)

X = sim.nfl[,'X'] ; Y = sim.nfl[,'Y']
plot(sim.nfl[,'time'], X, t='l', ylim= c(min(X,Y), max(X,Y)), col='red', xlab='time', ylab='')
lines(sim.nfl[,'time'], Y, col='blue')

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

3.2 Paramètres du modèle

• \(\gamma\) : nombre de graines produites par une plante
• \(\alpha\) : proportion de graines produites l’année précédente qui germent en mai
• \(\beta\) : proportion de graines d’il y a 2 ans qui germent en mai
• \(\sigma\) : proportion de graines qui survivent à l’hiver

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

3.4 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 :

• \(\gamma = 2\)
• \(\alpha = 0.5\)
• \(\beta = 0.25\)
• \(\sigma = 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 \(\gamma\))

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

Réalisation

3.4.1 Schéma du système

Schéma

3.4.2 Variables d’état

• \(P_n\) : Le nombre de plantes à l’année \(n\)
• \(S^1_n\) : Le nombre de graines de l’année précédente à l’année \(n\)
• \(S^2_n\) : Le nombre de graines ayant 2 ans à l’année \(n\)

3.4.3 Equations

\[ \begin{cases} P_n = \alpha S^1_n + \beta S^2_n \\ S^1_{n+1} = \sigma \gamma P_n \\ S^2_{n+1} = \sigma (1 - \alpha) S^1_n \end{cases} \]

3.4.4 Implémentation

parameters = c( a=0.5, b=0.25, g=2, s=0.8)

syst=function(t, state, parameters) {
    with(as.list(c(state, parameters)),{
        P = S1 * a + S2 * b
        dS1 <- s*g*P   # P plantes * gamma * sigma 
        dS2 <- s*(1-a)*S1 # S1 * (1-a) * sigma

        dS1 = dS1 - S1
        dS2 = dS2 - S2

        list(c(dS1, dS2))
    })
}

3.4.5 Simulations

A partir de 100 graines

state = c( S1=100, S2=0)

times=0:20
out <- ode(y = state, times = times, func = syst, parms = parameters, method="euler")
# method euler car le pas de temps est discret cf. man
head(out) ; tail(out)

##      time     S1     S2
## [1,]    0 100.00  0.000
## [2,]    1  80.00 40.000
## [3,]    2  80.00 32.000
## [4,]    3  76.80 32.000
## [5,]    4  74.24 30.720
## [6,]    5  71.68 29.696

##       time       S1       S2
## [16,]   15 50.55520 20.94065
## [17,]   16 48.82042 20.22208
## [18,]   17 47.14517 19.52817
## [19,]   18 45.52740 18.85807
## [20,]   19 43.96515 18.21096
## [21,]   20 42.45650 17.58606

par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-", ylog=TRUE)
mtext(outer = TRUE, side = 3, "plantes annuelles", cex = 1.5)

A partir de 50 graines

state = c( S1=100, S2=0)

times=0:150
out <- ode(y = state, times = times, func = syst, parms = parameters, method="euler")
# method euler car le pas de temps est discret cf. man
min(which(out[,'S1']<1))

## [1] 129

par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-", ylog=TRUE)
mtext(outer = TRUE, side = 3, "plantes annuelles", cex = 1.5)

3.4.6 \(\gamma\) minimal pour le maintient des pantes

On cherche \(\gamma\) pour \(\frac{P_{n+1}}{P_n} \ge 1\)

\[\frac{\alpha S^1_{n+1}+\beta S^2_{n+1}}{P_n} = \frac{ \alpha (\sigma\gamma P_n) + \beta (\sigma (1 - \alpha) S^1_n) }{P_n}\]

\[= \frac{ \alpha \sigma \gamma P_n }{P_n} + \frac{ \beta \sigma (1-\alpha)S^1_n}{P_n}\]

\[= \alpha \sigma \gamma + \beta \sigma (1-\alpha) \frac{S^1_n}{P_n}\]

\[= \alpha \sigma \gamma + \beta \sigma (1-\alpha) \frac{\sigma \gamma P_{n-1}}{P_n}\]

\[= \alpha \sigma \gamma + \beta \sigma^2 (1-\alpha) \gamma \frac{P_{n-1}}{P_n}\]

\[= \gamma (\alpha \sigma + \beta \sigma^2 (1-\alpha) ) \frac{P_{n-1}}{P_n} \ge 1\]

d’où

\[\gamma \ge \frac{1}{ (\alpha \sigma + \beta \sigma^2 (1-\alpha) ) \frac{P_{n-1}}{P_n} }\]

avec \(\frac{P_{n-1}}{P_n} =1\) on a

\[\gamma \ge \frac{1}{ (\alpha \sigma + \beta \sigma^2 (1-\alpha) )}\]

soit :

with( as.list(parameters), {1 / (a*s + b*s**2*(1-a))})

## [1] 2.083333

a=0.5 ; b=0.25 ; s=0.8
1 / (a*s + b*s**2*(1-a))

## [1] 2.083333