Compte rendu des TP de l’UE MDET: Dynamique adaptative

Importation des packages utiles

rm(list=ls()) # Nettoyage de l'environnement
# Liste des packages requis
packages <- c("deSolve", "tidyverse", "rgl", "cowplot", "tidyr")

# Vérification des packages installés
for (i in packages) {
  if (!requireNamespace(i, quietly = TRUE)) {
    # Installation du package s'il n'est pas installé
    install.packages(i)
  }
}

# Chargement des packages
library(deSolve) # package de la fonction ode pour les equation différentielles ordinaires
library(tidyverse) # pour les manupilation et les graphiques
library(rgl) # pour les graphiques en 3D
library(cowplot)
library(tidyr)

1 Visualisation des fonction \(k\), \(a\) et \(s\) selon le phénotype

Tout d’abord regardons l’allure des fonctions \(k\), la capacité de charge pour chaque phénotype, et la fonction de compétition \(a\) pour différentes valeurs de \(\sigma^2\).

a <- function(x1,x2, sigmaA){
  return(exp(-.5*(x1-x2)**2/sigmaA))
}

K <- function(x){
  return(K0 - lambda*(x - x0)**2)
}

x0=1
K0=1
lambda = 1

x = seq(0.1,1.9,.1)

sigma2 = c(0.1, 0.5, 1, 10)
plot_list = list()
for (i in 1:length(sigma2)){
  sigmaA = sigma2[i]
  data1 = data.frame(X=x, val_K = K(x), 
                     A_x5=a(x[5],x, sigmaA),
                     A_x9=a(x[9],x, sigmaA),
                     A_x10=a(x[10],x, sigmaA), 
                     A_x11=a(x[11],x, sigmaA),
                     A_x15=a(x[15],x, sigmaA)) %>%
    pivot_longer(-X, names_to="func",values_to = "value")

  plot_list[[i]] = ggplot(data=data1, aes(x=X, y=value, col=func))+
    geom_line() +
    xlab("Phénotype")+
    ggtitle(paste0("Allure des fonctions pour sigma²=", sigmaA))+
    theme_bw()+
    theme(line = element_blank(), 
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        axis.ticks =  element_line(colour = "black"),
        axis.text.x = element_text(colour = "black", size=10),
        axis.text.y = element_text(colour = "black", size=10),
        axis.title=element_text(size=10),
        legend.text = element_text(size = 7),  
        legend.title = element_text(size = 7)) +
    scale_color_manual(
      name = "Fonction",
      values = c("A_x10" = "red","A_x9" = "orange","A_x11" = "magenta", "A_x5" = "blue", "A_x15" = "green", "val_K" = "black"),
      labels = c("A(1,x)","A(0.9,x)","A(1.1,x)", "A(0.5,x)", "A(1.5,x)", "K(x)"))
    
  
}

plot_grid(plotlist = plot_list)

Lorsque \(\sigma^2\) est faible la courbure pour la fonction de compétition est assez faible et plusieurs phénotypes semblent pouvoir cohabiter à des valeurs de K non maximales (en quelque sorte rentrer sous K), plus \(\sigma^2\) augmente, plus la fonction de compétition s’applatit et la coexistence devrait devenir plus compliquée. On s’attend donc à une relation négative entre \(\sigma^2\) et le nombre de phénotype pouvant coexister.

On va maintenant faire une représentation en 3D pour la fitness d’invasion, qui s’exprime : \(s(x_{1},x_{2})=r.(1-a(x_{1},x_{2}).\frac{K(x_{2})}{K(x_{1})})\)

a <- function(x1,x2, sigmaA){
  return()
}

S <- function(x1,x2){
  r=1
  a = exp(-.5*(x1-x2)**2/sigmaA)
  return(r*(1-a*K(x2)/K(x1)))
}
x0=1
K0=1
lambda = 1
sigmaA = 1

x1 <- seq(0.1,1.9,.1)
x2 <- seq(0.1,1.9,.1)
#create z-values
z = outer(x1, x2, S)

#create 3D plot
persp(x1, x2, z, xlab='x', ylab='y', zlab='s(x,y)',
      main='Invasion fitness', col='lightblue', shade=.4, theta = 60, phi = 25, ticktype='simple')

2 Simulations de la dynamique adaptative

2.1 Définition des fonctions utilisées dans la suite

• Créations des fonctions de base pour la simulation

# Fonction pour calculer K
calculeK <- function(X,K0,x0, lambda){
  return(max(K0 - lambda*(X - x0)**2, 10**(-9)))
}


# Fonction pour obtenir les valeurs de la matrice A
calculeA <- function(X, nbSpecies, sigmaA_squared){
  D = matrix(data=X,ncol=nbSpecies,nrow=nbSpecies) # matrice
  C = t(D) - D
  return(exp(-C**2/(2*sigmaA_squared)))
}


# Fonction pour calculer la matrice B
calculeB <- function(X, nbSpecies, K0, x0, lambda){
  vectK = sapply(X, FUN= function(x) { calculeK(x, K0, x0, lambda) })
  B = matrix(data = 1/vectK, nrow=nbSpecies, ncol=nbSpecies, byrow=TRUE)
  return(B)
}


# Fonction par calucle la matrice M qui est le produit de A et B
calculeM = function(X, nbSpecies, sigmaA_squared, K0=1, lambda=1){
  x0 = min(X) + (max(X) - min(X)) / 2  
  
  A = calculeA(X, nbSpecies, sigmaA_squared) # matrice A
  B = calculeB(X, nbSpecies, K0, x0, lambda) # matrice B
  
  M=A*B # produit terme a terme de A avec B
  return(M)
}


# Fonctions pour définition l'equation differentielle 
calculeGradientLV = function(t, Ni, params){
  with(as.list(params), {
    
    dNi = r*Ni*(1 - Ni%*%M)

    return(list(dNi))
  })
}


# Fonctions pour resoudre l'equation differentielle avec ode du package deSolve
resoutEqd <- function(parameters, densite_t0, vect_time, vect_phenotypes){
  df_densite = ode(y = densite_t0, times = vect_time, func = calculeGradientLV, parameters) %>%
    as.data.frame() %>%
    pivot_longer(-time, values_to = "density", names_to = "species")

  df_densite$X = vect_phenotypes[as.numeric(df_densite$species)]
    
  return(df_densite)# retourne sortie modifier de ode
}


# Définition de la fonction pour simuler le dynamique de la poppulation
simuleDynamique <- function(nbSpecies, sigmaA_squared, vect_temps, K0=1, r=1){
  n0 = rep(K0/nbSpecies, nbSpecies) 
  vect_phenotypes = seq(0, 1, length.out = nbSpecies)

  M = calculeM(vect_phenotypes, nbSpecies, sigmaA_squared, K0)
  params = list(r = r, M = M)
  df_densite = resoutEqd(params, n0, vect_temps, vect_phenotypes)
  
  return(df_densite) # la densité de la population obtenu avec ode
}

• Les fonctions pour les représentations graphiques

# définition des thèmes
main_theme = theme_bw()+
  theme(line = element_blank(), 
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        axis.ticks =  element_line(colour = "black"),
        axis.text.x = element_text(colour = "black", size=10),
        axis.text.y = element_text(colour = "black", size=10),
        legend.position="none",
        axis.title=element_text(size=10))

# fonction pour illustrer le dynamique de la population selon le temps
plotDynamique <- function(df_density, plot_theme=main_theme, title= ""){
  plot = ggplot(df_density)+
    geom_line(aes(time, density, col = species)) +
    ggtitle(label = title) +
    xlab("Temps") +
    ylab("Phénotypes") +
    plot_theme
  return(plot)
}

# fonction pour illustrer l'évolution de la dynamique selon les phénotypes
plotEvolutiontrait <- function(df_density, plot_theme=main_theme){
  plot = ggplot(data = df_density) +
    geom_tile(aes(x = X, y = time, fill = density))+
    scale_fill_gradient2(low = "white" , high = "red")+
    xlab("Phénotypes") +
    ylab("Temps") +
    plot_theme
  
  return(plot)
}

2.2 Simulation sans mutation

vect_temps = 0:2000 # temps de simulation
K0 = 1             # capacite de charge du milieu tout phenotype confondu
sigmaA_squared = 1 # varie entre 0.1 et 10
r=1                # taux de croissance, le même pour toutes les espèces


# Simulation de la dynamique avec 10 phenotypes
nbSpecies = 10 # nombre d'especes
df_densite10 = simuleDynamique(nbSpecies, sigmaA_squared, vect_temps)

plot1 = plotDynamique(df_densite10, title="Evolution des phénotypes (N = 10)" )
plot2 = plotEvolutiontrait(df_densite10)


# Simulation de la dynamique avec 100 phenotypes
nbSpecies = 100 # nombre d'especes
df_densite100 = simuleDynamique(nbSpecies, sigmaA_squared, vect_temps)

plot3 = plotDynamique(df_densite100, title="Evolution des phénotypes (N = 100)")
plot4 = plotEvolutiontrait(df_densite100)

plot_grid(labels = "auto",plotlist=list(plot1,plot2,plot3,plot4))

Avec \(\sigma^2 = 1\), une valeur intermédiaire, on observe que toutes les espèces s’éloignant trop de l’optimum sont éliminées. La tendance semble donc plutôt être l’élimination par compétition que la coexistence de phénotype.

2.3 Coexistence des phénotypes et valeurs de sigma

On peut se poser la question suivante : Combien d’espèces persisteront dans le temps selon la valeur de sigma ?

vect_temps = 0:2000 # temps de simulation
sigmaA_squared = c(0.1, 0.5, 1, 10) # Les differents valeurs de sigma

plot_list = list()

nbSpecies = 100
for (i in 1:length(sigmaA_squared)){
  df_densite = simuleDynamique(nbSpecies, sigmaA_squared[i], vect_temps)
  plot_list[[i]] = plotDynamique(df_densite, title=paste("Cas N= 100, sigma² = ", sigmaA_squared[i]))
  #plot_list[[i]] = plotEvolutiontrait(df_densite))
}

plot_grid(labels = "auto",plotlist = plot_list)

Avec \(N = 100\) On constate qu’à mesure que \(\sigma^2\) augmente, le nombre de phénotypes qui persistent à long terme augmente entre 0.1 et 0.5. Lorsque \(\sigma^2\) est supérieur à \(1\), certains phénotypes subsistent initialement mais finissent rapidement par s’éteindre, et au bout de \(2000\) pas de temps de simulation, un seul phénotype subsiste.

Regardons maintenant le nombre de phénotypes qui persistent selon plusieurs valeurs de \(sigma²\).

density_min = 0.005 # en dessous de cette densité le phénotype est considéré disparu

vect_temps = 0:2000 # temps de simulation
sigmaA_squared = c(seq(0.1,1, by=0.1), seq(2,10)) # les valeurs de sigma² a tester

df_persitance = data.frame()
nbSpecies = 100 # Nombre d'espèces
for (i in 1:length(sigmaA_squared)){
  df_densite = simuleDynamique(nbSpecies, sigmaA_squared[i], vect_temps)
  nb_persistant = df_densite[df_densite$time == max(df_densite$time),] %>%
    filter(density>density_min) %>%
    nrow()
  
  df_persitance = rbind(df_persitance,c(sigmaA_squared[i],nb_persistant))
}

colnames(df_persitance) = c("sigma_square", "nb_persistant")

ggplot(data = df_persitance, aes(x = sigma_square, y = nb_persistant)) +
  geom_point() +
  geom_line() +
  labs(title = "Nombre de phénotypes persistants selon sigma²",
       subtitle = "Cas N = 100",
    x = "Sigma²",
    y = "Nombre de phénotypes"
  ) +
  main_theme

Comme on s’y attendait le nombre d’espèces persistantes diminue lorsque \(\sigma^2\) augmente.

Y a-t-il une limite au nombre d’espèces qui peuvent coexister dans ce modèle ? A priori il n’y pas de limite car on prend une fonction gaussienne, mais une limite pourrait apparaitre si on prenais d’autres distributions.

2.4 Simulation avec ajout de mutations dans la dynamique

Définition des fonction pour simuler la mutation.

Pour simuler des mutations, le principe va être le suivant :

1. On tire les temps auxquels les mutations vont avoir lieu de manière aléatoire

2.a On simule la dynamique entre t0 et le premier temps de mutation

2.b Au temps de mutation, on récupère les densités finales puis on applique la mutation, qui consiste à “décaler” les phénotypes vers la droite ou la gauche en déplaçant une certaines proportion de la densité dans le phénotype juste à droite ou juste à gauche.

2.c On retourne à l’étape 2.a et on simule la dynamique entre le temps de mutation et le suivant en prenant comme condition initiale le vecteur des densités mutées.

Définition des fonctions utiles :

# fonction pour obtenir les temps auxquels il y aura une mutation 
getMutationTime <- function(vect_temps, taux_mutation){
  SimulationDuration = max(vect_temps)
  # Determination des instants auxquels il y aura une mutation
  nb_mutations_moy = SimulationDuration * taux_mutation 
  nb_mutations_sd = nb_mutations_moy / 4
  nb_mutations = max(1,floor(rnorm(1, mean=nb_mutations_moy, sd=nb_mutations_sd)))
  
  temps_mutation = c(0,sort(sample(x = vect_temps, size = nb_mutations)),SimulationDuration)
  return(temps_mutation)
}


# fonction pour le décalage (droite ou gauche) de la mutation
mutationParDecalage <- function(df_densite, proportionMutant) {
  
  n0 = df_densite$density[df_densite$time==max(df_densite$time)] 
  nbSpecies = length(n0)
  # mutation : decalage d'une partie des effectifs
  p = runif(1) # tirage aléatoire entre 0 et 1

  if (p < .7){
  # decalage à droite
  effectifMutants = n0[1:(nbSpecies-1)]*proportionMutant
    
  n0[1:(nbSpecies-1)] = n0[1:(nbSpecies-1)] - effectifMutants
  n0[2:nbSpecies] = n0[2:nbSpecies] + effectifMutants
  }
  else{
  # decalage a gauche
  effectifMutants = n0[2:nbSpecies]*proportionMutant
    
  n0[1:(nbSpecies-1)] = n0[1:(nbSpecies-1)] + effectifMutants
  n0[2:nbSpecies] = n0[2:nbSpecies] - effectifMutants
  }
  return(n0)
}


# fonction pour la simulation de la mutation
simuleDynamiqueMutation <- function(nbSpecies, sigmaA_squared, temps_mutation, K0=1, r=1){
  
  # Calcule des phenotypes et de M
  vect_phenotypes = seq(0, 1, length.out = nbSpecies)
  M = calculeM(vect_phenotypes, nbSpecies, sigmaA_squared, K0)
  params = list(r = r, M = M)
  
  # Initialisation de n0 et de df_densite
  n0 = c(K0, rep(0, nbSpecies-1))
  df_densite = data.frame()

  # Simulation de la dynamique pour tout les temps de mutation
  for (i in 2:length(temps_mutation)){
    if (temps_mutation[i]-temps_mutation[i-1]>1){
      # Entre 2 mutation on simule en resolvant l'equation differentielle
      sample_vect_temps = seq(temps_mutation[i-1],temps_mutation[i]-1)
      df_densite = rbind(df_densite,resoutEqd(params, n0, sample_vect_temps, vect_phenotypes))
      # A chaque temps contenu dans temps_mutation une mutation a lieu 
      n0 = mutationParDecalage(df_densite, proportionMutant) # no devient le no avec mutation
    }
    
  }
  return(df_densite) # pop tot avec mutation
}

Visualisation de la dynamique avec 10 et 100 espèces, et un taux de mutation de 1/100, soit en moyenne une mutations tous les 100 pas de temps.

vect_temps = 0:2000 # Temps de simulation
sigmaA_squared = 0.4 # simga²
proportionMutant = 0.1 # Proportion du nombre de phénotype mutant
taux_mutation = 1/100 # taux de mutation


# Simulation de la dynamique avec 10 phenotypes
nbSpecies = 10
temps_mutation = getMutationTime(vect_temps, taux_mutation)
df_densite = simuleDynamiqueMutation(nbSpecies, sigmaA_squared, temps_mutation)

plot1 = plotDynamique(df_densite, title="Cas N = 10") 
plot1 = plot1 + geom_vline(xintercept=temps_mutation, linetype='18', color='red')

plot2 = plotEvolutiontrait(df_densite)



# Simulation de la dynamique avec 100 phenotypes
nbSpecies = 100
temps_mutation = getMutationTime(vect_temps, taux_mutation)
df_densite100 = simuleDynamiqueMutation(nbSpecies, sigmaA_squared, temps_mutation)
  
plot3 = plotDynamique(df_densite100, title="Cas N = 100") 
plot3 = plot3 + geom_vline(xintercept=temps_mutation, linetype='18', color='red')

plot4 = plotEvolutiontrait(df_densite100)

plot_grid(plotlist=list(plot1,plot2,plot3,plot4))

Lorsque l’on fait une simulation avec des mutation, on s’attend à ce que les mutations ramènent progressivement les valeurs de phénotypes vers l’optimum. Dans certains cas, en dynamique adaptative on peut observer des points de branchement. Pour cela il faut qu’un point singulier de la dynamique soit attracteur mais instable, et qu’il puisse être envahi par des phénotype de valeur inférieure et supérieure à sa valeur.

On a l’impression d’oberver un point de branchement avec nos 10 phénotypes. En revanche, on observe que avec 100 phénotypes, on n’a pas le temps de rejoindre l’optimum. On va refaire une simulation avec 100 phénotypes en prenant 2 taux de mutations différents.

vect_temps = 0:10000
sigmaA_squared = 0.4
proportionMutant = 0.1


# Simulation de la dynamique avec 50 phenotypes et taux_mutation = 1/200
nbSpecies = 50
taux_mutation = 1/200
temps_mutation = getMutationTime(vect_temps, taux_mutation)
df_densite = simuleDynamiqueMutation(nbSpecies, sigmaA_squared, temps_mutation)

plot1 = plotDynamique(df_densite, title="Taux de mutation = 0.005") 
plot1 = plot1

plot2 = plotEvolutiontrait(df_densite)



# Simulation de la dynamique avec 100 phenotypes et taux_mutation = 1/50
nbSpecies = 50
taux_mutation = 1/50
temps_mutation = getMutationTime(vect_temps, taux_mutation)
df_densite100 = simuleDynamiqueMutation(nbSpecies, sigmaA_squared, temps_mutation)
  
plot3 = plotDynamique(df_densite100, title="Taux de mutation = 0.02") 
plot3 = plot3

plot4 = plotEvolutiontrait(df_densite100)

plot_grid(plotlist=list(plot1,plot2,plot3,plot4))

Avec un taux de mutation faible, la densité est toujours très resserrée autour de peu de phénotype dominant (le trait est bien rouge et pas étalé) et on se déplace très lentement vers l’optimum. Si le taux de mutation est plus élevé en revanche, on observe un point de branchement et plus de phénotypes qui coexistent en plus faible densité.