Démystifier l'Ordinateur en Biologie

1) reproduction: modèle médical: chimpanzé modèle de l'homme (modèle réduit)

2) modèle mathématique: les modèles les plus simples sont des équations exprimant une variable en fonction d'autres variables

Evidence de la notion de modèle - vision du monde:

- modèle Newtonien de la gravité, ce modèle, sa géométrie et sa mathématique font partie de la vision du monde d'un étudiant de la fin du XXème siècle.

- modèle quantique de la nature: basé sur une algèbre moins immédiate, ce modèle n'est appréhendable que par équations non "visualisables" (sauf pour des spécialistes...).

Un modèle ne peut s'imposer ou être évident à priori, sauf à avoir une vision du monde dans lequel il se traduit immédiatement (ex: mécanique Newtonienne, forces, vecteurs).

Un modèle n'est que la représentation dans notre esprit (ou sur une page d'équation) d'une partie de la réalité.

Attention: Différence entre construction du modèle et ajustement du modèle.

On peut en général ajuster "assez bien" un modèle (même si il est manifestement faux). Mais un modèle qui décrit parfaitement le phénomène est un bon modèle. L'ordinateur joue un rôle central dans la modélisation dans le sens où il permet la simulation de modèle, et l'ajustement des modèles.

Problème de l'ajustement : erreurs de mesures, mais surtout petites variations aléatoires intrinsèques, ce qui fait qu'on n'aura jamais un modèle s'ajustant parfaitement à tous les points expérimentaux.

Exemple: surface d'une feuille par rapport à ses longueurs et largeurs (x = L*l). On attend (si la forme est la même) une relation du type Surface = b * x, le modèle est donc déterministe.

En fait, beaucoup de facteurs influencent la croissance et la forme finale de la feuille : on aura finalement une relation du type Surface = b * x + e.

Le e représente les déviations par rapport à la droite, et on obtient donc un modèle stochastique.

A partir des données expérimentales, on ne pourra obtenir qu'une estimation du paramètre b.

Le modèle de feuille vu ci dessus ne représente qu'un aspect de la feuille.

On pourrait (Lindenmayer,1975) décrire la croissance de la feuille par un modèle de croissance simple qui peut donner lieu à une description complexe : le point de départ est qu'une feuille est un assemblage de cellules dans différents états (a,b,c,d,k). La croissance a lieu par pas successifs, par les cellules du bord de la feuille, et à chaque pas, des règles de remplacement simples vont être appliquées à ces cellules du bord de la feuille :

a Þ cbc, b Þ dad, c Þ k, d Þ a, k Þ k.

1 {a}

2 {c, b, c}

3 {k, d, a, d, k}

4 {k, a, c, b, c, a, k}

5 {k, c, b, c, k, d, a, d, k, c, b, c, k}

6 {k, k, d, a, d, k, k, a, c, b, c, a, k, k, d, a, d, k, k}

7 {k, k, a, c, b, c, a, k, k, c, b, c, k, d, a, d, k, c, b, c, k, k, a, c, b, c, a, k, k}

8 {k, k, c, b, c, k, d, a, d, k, c, b, c, k, k, k, d, a, d, k, k, a, c, b, c, a, k, k, d, a, d, k, k, k, c, b, c, k, d, a, d, k, c, b, c, k, k}

Si on interprète a et b comme des pointes, c et d comme les bords de la pointe, et k comme les bords creux de la feuille, on peut voir qu'à la génération 6, on retrouve au milieu le motif de la génération 4 {k, a, c, b, c, a, k}. De même, à la génération 8, on retrouve au milieu le motif de la génération 6 {k, k, d, a, d, k, k, a, c, b, c, a, k, k, d, a, d, k, k}.

Distinction entre modèle dynamique et modèle non dynamique :

Modèle dynamique: fait intervenir le temps (exemple des feuilles)

Modèle non dynamique: équilibre (le temps n'intervient pas).

En fait, il y a un continuum entre modèles dynamiques et modèles non dynamiques. Les échelles de temps de relaxation des phénomènes varient selon ceux ci: électro-mécaniques (décharges neuronales, ms), métabolisme (enzymes, sec à minutes), épigénétiques (régulation à court terme des concentrations d'enzymes, minutes à heures), phénomènes de développement (différentiation, heures à années) et d'évolution (mouvement de matériel génétique à travers une population, mois à années).

Le calcul informatique est nécessaire en modélisation biologique pour deux choses:

- la simulation d'un modèle (mettre en jeu des équations dynamiques, voir d'une manière qualitative si le modèle répond comme "la réalité").

- le calcul des paramètres d'un modèle (une fois le modèle construit, on ajuste les paramètres sur les données expérimentales).

Biologie des populations:

Modèles de croissance:

Croissance exponentielle

Colonie au temps t+delta_t = colonie au temps t + incrément de croissance pendant (t,t+delta_t)

ou symboliquement:

N(t+delta_t) = N(t) + delta_N

L'accroissement delta_N est proportionnel à la population, c'est à dire que delta_N = r N delta_t

r est le taux de croissance.

En considérant un intervalle de temps infinitésimal, (i.e. en faisant tendre delta_t vers 0, soit dt), on a donc l'équation différentielle suivante : dN/dt = r N ou dN/N = r dt et donc si on intègre des deux côtés:

Ce modèle simple décrit bien la croissance de E. coli dans la phase de démarrage d'une culture.

Ce modèle est continu, c'est à dire que l'on ne différencie pas les individus (qui devraient être des nombres entiers en ce cas), ni le temps de générations.

Si on prend un insecte se reproduisant annuellement, on aurait :

Ni+1 = R Ni

En démarrant avec N0, on aurait au bout de i générations, Ni= R^i N0, le paramètre R est appelé le taux de croissance géométrique (Note: l'operateur ^ signifie "puissance"). Si on l'écrit sous la forme R = elog(R) , on obtient : Ni = e^i log(R) N0, et la population s'incrémente de façon exponentielle avec i, le nombre de générations.

Ce modèle de croissance est évidemment idéal dans le sens où une vraie population est limitée par la niche écologique qu'elle habite (ou qu'elle peut envahir). Manque de nourriture, manque de place ou accumulation de toxiques limitent donc la population.

On arrive à un modèle de croissance limitée.

Croissance limitée

Modèle chimique monomoléculaire:

R -> C

Si R(t) est le réactif au temps t, (k constante de la réaction en concentration/unité de temps). La solution est , ou R(0) est le réactif au temps t=0.

Si C(t) est le produit de la réaction, on a .

Graphe de C(t) en fonction du temps:

Modèle de population à croissance limitée:

Le modèle précédent (croissance exponentielle) nous donnait . Maintenant le taux de croissance décroît linéairement avec la taille de la population, ce qui donne :

Solution analytique:

-> log (N) - log (K-N) = r t + c (c: constante d'intégration)

Si N(t) = N(0) à t=0, on a :

log (N) - log (K-N) = r t + log (N(0)) - log(K-N(0)) = log(exp(rt)) + log (N(0)) - log(K-N(0))

et donc:

C'est l'équation de Verlhust. Elle est parfois écrite sous la forme :

avec K/N(0) - 1 = exp(hr)

Graphe de N(t) en fonction de t (avec différentes populations initiales N(0)):

On remarque que lorsque la population initiale excède la population maximale théorique, elle diminue et tend à se rapprocher de celle ci.

On peut encore raffiner le modèle logistique. En effet, ce modèle impliquait que les facteurs limitant le taux d'accroissement agissaient immédiatement. Or, de manière plus réaliste, on peut imaginer qu'à un moment donné, le taux de variation de la population dépend non plus de la population à ce même moment, mais de la population à un moment antérieur. Ceci peut conduire à un modèle plus réaliste, exemple : population d'herbivores qui s'accroîtrait en mangeant toutes les réserves, les naissances sont aussi fortes alors que les facteurs limitant (nourriture) n'est plus là). On parle du modèle logistique avec delai :

En ce cas, on ne pourra plus utiliser de méthode analytique pour intégrer, on utilisera une méthode numérique (voir plus loin).

Cela donne plusieurs comportements, allant de la courbe "standard" de Verlhust (rT=0.30), à une variation atténuée au fil du temps (rT=1.2), ou à des cycles limites stables (rT=1.8)

.

Modèles discrets de croissance :

(Vus en TP tableur au DOB).

Ces modèles sont des modèles discrets de croissance.

Modèle discret exponentiel :

Ni+1 = R Ni

Modèle discret à croissance limitée:

Le paramètre K mesure l'effet des facteurs limitant (note: Ni ne peut être plus grand que K, sinon il y aurait croissance négative!).

Si K est grand, la croissance est exponentielle (au moins au début, quand Ni est négligeable devant K).

Souvent on écrit ce modèle comme:

Xi+1 = R Xi (1-Xi)

et donc la variation s'exprime comme:

Delta_Xi = Xi+1 - Xi = RXi - RXi^2 - Xi = R Xi (1-Xi-1/R) = R Xi ((R-1)/R - Xi)

et donc Delta_Xi = R Xi (XE -Xi) avec XE = (R-1)/R

Le paramètre XE est X à l'équilibre, car si Xi = XE , la population ne varie pas.

On a vu en TD tableur les différents cas selon le paramètre R.

Si 1<R<3, stable

Si 3<R<3.45, cycle limite.

Si 3.57 < R < 4, régime chaotique.

Donc, la discrétisation amène un comportement différent du modèle continu.

Ajustement de points expérimentaux à un modèle théorique:

Notion d'ajustement de modèle:

Problème:

On construit un modèle, c'est- à-dire qu'on a une fonction théorique f d'une variable mesurée ou dérivée expérimentalement x. Par exemple, la feuille du début du cours, x = L*l (longueur*largeur)) et sa surface théorique S(x) = b*x

Par la suite, on appellera les mesures expérimentales les yi, (dans le cas de la feuille, les surfaces mesurées par pesée par exemple). A chaque yi correspond un xi. Si le modèle est bon, et que les spécimens suivent exactement le modèle (pas de bruit dans la mesure, ni de variation à la règle), les yi devraient être égaux aux f(xi). Ce n'est évidemment jamais le cas (même en physique !).

Le modèle qu'on va effectivement avoir c'est:

f(xi) = yi + e avec e représentant une variation aléatoire, petite devant yi. Le modèle théorique f(x) possède des paramètres a,b,c,...Dans l'exemple de la feuille f(x) = b*x, on a un seul paramètre, le paramètre b.

Pour que le modèle théorique représente la "réalité" qu'on veut modéliser, il va falloir ajuster les paramètres du modèle pour "coller" à cette réalité expérimentale. C'est-à-dire trouver des valeurs de paramètres qui superposent la courbe théorique f(x) avec les yi mesurés.

Le problème qu'on va donc se poser, c'est, étant donné un modèle f(xi), faire varier a,b,c,... pour que les f(xi) soient le plus "proches" possibles des yi, et ce tous ensembles.

Choix du critère:

On pourrait prendre la somme des différences yi - f(xi) et chercher à minimiser cette somme: cela n'est pas correct et nous n'atteindrons pas notre but, car une différence négative va compenser une différence positive, ce qui donnera un résultat manifestement faux.

Exemple: Sur la figure la différence algébrique entre la droite rouge et les points est du même ordre que celle de la droite verte avec les points, alors que l'ajustement de la droite verte sur les points est visiblement meilleure que l'alignement de la droite rouge.

Pour corriger cet effet algébrique, on pourrait prendre la somme des différences en valeur absolue : cela pose des problèmes numériques.

On va en fait prendre la somme des carrés des écarts yi-f(xi). Ces écarts aux carrés seront toujours positifs et d'autant plus grands que les points expérimentaux seront éloignés du modèle (=points théoriques).

Voilà une analogie physique de la somme des carrés des écarts: relâchement et équilibrage des ressorts

Utilisation du carré des écarts lors d'un ajustement de courbe

Une explication intuitive de la minimisation de la somme des carrés des écarts est donnée dans la figure ci-dessus

Tout se passe comme si l'on attachait des ressorts (de longueur initiale nulle) entre les points (expérimentaux) et la fonction théorique. On recherche la courbe théorique (ici une droite) de paramètres a et b (pente, intersection à l'origine) qui minimisent l'ensemble des tensions des ressorts (l'énergie d'un ressort est bien en k * x^2). C'est en fait la droite qui s'ajuste le mieux sur les points (expérimentaux).

Or un ressort étiré a une force de rappel proportionnelle à sa longueur, qui mesure ici y-f(x), y étant l'ordonnée du point (croix), f(x) la fonction cherchée (ici fonction linéaire). Et l'énergie de tension pour chaque ressort est proportionnelle à (y-f(x))^2: ce qui fait que la minimisation de la somme des carrés des écarts à la droite revient à équilibrer et à minimiser l'énergie de tension des ressorts.

Justification théorique du carré des écarts:

Si on considère que les points sont disposés suivant une distribution gaussienne autour de leur position théorique (courbe théorique), on peut évaluer la probabilité de chaque point d'être entre les valeurs y et y+dy, c'est une gaussienne centrée sur la "vraie" valeur (valeur théorique f(x)) : elle est . La variance dépend du bruit ou du nombre de mesures faites. On considère que les erreurs sont indépendantes et identiquement distribuées (iid), et donc la probabilité de l'ensemble des données est le produit des probabilités d'avoir chaque point :

Le but est de trouver les parmètres de f(x) telle que la probabilité d'avoir tous ces points expérimentaux soit maximale.

Si on prend le log, on a :

Le log étant une fonction monotone croissante, si on veut maximiser la probabilité P, il faut maximiser log(P) et étant donné que est constante, on doit donc minimiser la somme (puisque le facteur est constant, et que la somme intervient de manière négative dans l'expression).

Maximiser P (c.-a-d. minimiser la somme des carrés des écarts) c'est maximiser la probabilité qu'on ait le modèle f(x) étant données les données yi.

Note: autres critères. Minimisation des Valeurs absolues. Minimisation du Maximum de différences. Comparaison avec Minimisation du Carré des écarts.

La droite (noire, 3ème en partant du bas) trouvée avec le minimum de la différence du carré des écarts est "attirée" par les points aberrants (la différence est magnifiée par le carré de y-f[x]). Celle (rouge, 2ème en partant du bas) qui minimise cette différence divisée par le carré des valeurs est moins attirée par les points aberrants.

La droite (verte, 1ère en partant du bas) donnée par la valeur absolue de la différence des écarts passe bien par la majorité des points (la différence avec les points aberrants n'est pas suffisamment attractive).

La droite (bleue, 4ème en partant du bas) donnant la différence minimum de la puissance infinie (ici 500) des écarts passe au mieux des points extrêmes.

Minimisation :

On cherche à minimiser la fonction S, somme du carré des écarts entre les yi et les f(xi).

Elle s'exprime sous la forme S = ·(yi - f(xi))^2

En anglais, on parle de "least square fitting": ajustement par les moindres carrés.

Et c'est bien les paramètres de f(x) qu'on va faire varier - pas n'importe comment au hasard, mais avec une méthode adaptée - pour ajuster la courbe sur les points (et pas l'inverse !).

Donc si la fonction f(x) a plusieurs paramètres a,b,c,… on va chercher le minimum de S en cherchant le point où sa dérivée(par rapport aux paramètres) s'annule en variant le modèle f (c'est à dire ses paramètres a,b,c,…).

Pour un seul paramètre, a, il faut faire:

dS/da = 0, puis calculer le paramètre a à partir de cette équation.

Si on a deux paramètres, a et b, il faudra dériver S sur les 2 paramètres:

dS/da = 0 et dS/db = 0, d'où on tire 2 équations à 2 inconnues (a et b), d'où on tire a et b.

On montre que le minimum d'une fonction à plusieurs paramètres est aussi un minimum sur chacun des paramètres. Note : la dérivée d'une somme est la somme des dérivées.

régression linéaire:

Dans le cas d'une fonction linéaire (f(x) = ax+b), on va montrer qu'on peut calculer a et b analytiquement (dans le cas général, fonction non linéaire, il faudra utiliser une méthode numérique).

. (on note et les estimations de a et b)On trouvera le minimum de cette fonction en cherchant le point où les 2 dérivées partielles (Attention: dérivées par rapport aux paramètres) s'annulent en même temps: on a un système de 2 équations à 2 inconnues:

(1)

(2)

(2) donne (on note la moyenne des yi et la moyenne des xi)

(1) donne en remplaçant par ,

d'où:

(3)

On a d'autre part, par définition, (la somme algébrique des différences autour de la moyenne est nulle), donc si on soustrait (qui est nul) au numérateur et (qui est aussi nul) dénominateur de (3), on obtient :

(on remarque en outre que la droite ajustée passe par le centroïde des points, puisqu’on a )

Note : la moyenne est l'estimateur qui minimise le carré des écarts :

Démonstration :

Ce qui est la définition de l’estimateur de la moyenne (cqfd).

Modèles linéaires et non linéaires

Modèles linéaires :

Les modèles linéaires sont ceux qui sont fonction linéaire des paramètres du modèle.

Exemple:

f(x) = a *x + b

f(x) = a *x^2 + b*x + c

Les autres sont des modèles non linéaires:

exemples:

f(x) = exp(-a x)

Note :

Résolution numérique versus résolution analytique :

On connaît la résolution analytique d'une équation du 2ème degré (et même du 3ème degré ou 4ème degré), c'est à dire une équation de la forme:

y = f(x) = a x^2 + bx + c

Mais on ne sait pas résoudre analytiquement :

y = f(x) = x exp(-x) ou même y = f(x) = x cos(x)

Dans ce cas, on utilise une méthode de résolution numérique.

Par exemple, dans le cas général d'une équation différentielle :

dN/dt= f(N(t),t) avec N(0) à t=0, la résolution numérique utilise l'équation aux différences:

Ñ(t+dt) = Ñ(t) + dt (approximation du premier ordre, développement de Taylor)

Où Ñ(t) est une valeur calculée qui est une approximation de la vraie valeur N(t).

On effectue une itération en partant de la valeur connue de N à t=0, c'est à dire N(0), pour donner N(dt), N(2dt), N(3dt), etc... L'approximation est d'autant meilleure que le pas dt est petit. Cette méthode revient, connaissant N(t), à calculer N(t+dt), ceci est une application pratique de la dérivée: et on voit graphiquement pourquoi l'approximation est meilleure quand le pas est petit.