Méthode de Turing et morphogénèse : illustration des équations de réaction-diffusion

L'appliquette sert surtout dans le cadre de la biologie théorique, voire la géographie ou la géologie si l'on s'intéresse à l'étude de la formation des dunes. Il n'entre donc pas vraiment dans le cadre astrophysique, du moins en l'état, à deux dimensions.

Rappels de cours

Bibliographie

Explications

Description du Modèle

Le modèle de Turing repose sur l'interaction entre éléments infinitésimaux, les cellules. Celles-ci peuvent être considérées comme analogues aux cellules vivantes de la peau d'un animal ou d'un embryon d'animal, par exemple. Dans un autre domaine, on pourrait considérer des grains de sable s'organisant dans des dunes.

On considère un réseau, une grille, de cellules qui peuvent être dans deux états possibles: un état activé ou désactivé. Cela s'accorde bien avec le fait que chez les mammifèes par exemple une tache noire serait dûe à l'activation de mélanine qui colore la peau. On part de l'idée de la coexistence de deux types de sécrétions pour les cellules :

Un effet activateur à courte portée
Un effet inhibiteur à longue portée

Ces deux effets, issus des cellules elle-mêmes, font que l'état d'une cellule à un instant donné dépend de celui des cellules autour d'elle. On peut le voir sur la figure ci-dessous: la cellule 2 reçoit l'inhibiteur de la cellule 1 (car elle est dans sa zone gris foncé, sa distance à la cellule 1 est donc inférieure à R2). Cependant la cellule 2 est hors de la zone d'activateur sécrété par la cellule 1 (zone gris clair), car sa distance à la cellule 1 est supérieure à R1. Elle n'est donc pas incitée à s'activer par la cellule 1.

Rayons d'action

Equation d'Evolution

Ces équations traduisent la diffusion de l'activateur et de l'inhibiteur sécrétés par une cellule et affectant ses différentes voisines.
On note :

R1 le rayon d'action de l'activateur
R2 le rayon de l'inhibiteur
J1 le coefficient de diffusion de l'activateur
J2 le coefficient de diffusion de l'inhibiteur
h le palier d'activation

$cellule_i(t+1) = signe\Biggl[h + J1 \sum \limits_{r_{ij}<R1}{cellule_j(t)} + J2 \sum \limits_{R1 \leq r_{ij}<R2}{cellule_j(t)}\Biggr]$

La formule ci-dessus signifie qu'on va examiner l'état des cellules (j) dans un rayon R1 pour prendre en comptre leur effet activateur et l'état des cellules (j) dans un rayon R2 pour prendre en comptre leur effet inhibiteur. Cela sert à calculer l'état de la cellule (i) qui est la cellule en cours d'actualisation.

Modes de mise à jour de la grille

On en a considéré trois, complémentaires. Le mode asynchrone actualise les cellules une par une de façon aléatoire, le paramètre "durée" de l'appliquette déterminant le nombre de cellules actualisées avant chaque nouvel affichage. Le mode synchrone calcule les états de toutes les cellules de l'image (n+1) à partir de celles de l'image n simultanément.

Le mode semi-asynchrone est le plus utile: il actualise l'image en calculant l'évolution cellule par cellule de façon asynchrone, mais considérant aléatoirement les cellules, il n'affiche la prochaine image que lorsque toutes les cellules ont été actualisées de la sorte. Cela permet d'éviter la présence d'une première constante de temps lors de l'apparition des structures qui n'était que le temps nécessaire pour visiter aléatoirement toute la grille.

Ci-dessous, on peut voir la stabilisation d'une image lors de l'apparition des taches; l'ordonnée représente la somme des valeurs de toutes les cellules (la valeur +1 est associée à une cellule activée, -1 à une cellule désactivée). Cette somme permet de suivre la stabilisation de l'automate au niveau global.

Stabilisation

Intérêt des Géométries

Les cellules étudiées sont placées sur des canevas, des maillages les liant entre elles. Cependant on peut voir ci-dessous un exemple d'émergence de structure de type empreinte digitale avec de forts effets horizontaux et verticaux: l'anisotropie est issue de la structure de la grille simple (rectangulaire). Donc la nécessité d'introduire d'autres structures apparaît.

simple

Structure fine de grille simple (Ract=1,Rinh=2, Jact=1, Jinh=-0.1, Offset=0)

Voici ce que donne la grille hexagonale. On voit que les structures horizontales et verticales liées à la géométrie de la grille simple ont disparu. Les structures sont également plus épaisses, ce qui est encore une preuve de l'influence du rayon d'action de l'inhibiteur sur la structure.

Hexagonale

Premier exemple de grille hexagonale (Ract=1,Rinh=5, Jact=1, Jinh=-0.1, Offset=0)

Un troisième maillage, celui de cellules de Voronoï (ou triangulation de Delaunay), vise à s'affranchir des artefacts dûs à la régularité du maillage.

Voronoy

Structure de maillage non régulier (Ract=3,Rinh=4, Jact=6.5, Jinh=-3, Offset=0)

Pour en savoir plus sur ces différents maillages, voir ici : Maillages. Ces grilles ont également été utilisées pour développer une appliquette concernant le modèle d'Ising.

Outils d'Analyse

Dans le but de déterminer la taille caractéristique des taches d'animaux ainsi obtenues, deux méthodes ont été développées. La première est basée sur l'étude de la fonction de corrélation de la grille et la deuxième sur l'analyse de la transformée de Fourier de l'image obtenue.

Fonction de Corrélation

$G(k) = \frac{1}{N(k)} \sum \limits_{l=1}^{N(k)}{\Biggl(\bigg|\sum \limits_{ i \in \leq k \atop j \in \leq k}{cellule_{ij}}\bigg|\frac{1}{k^2}\Biggr)} \quad \mbox{avec}\quad N(k) = \frac{dimx}{k} \times \frac{dimy}{k}$

Le calcul de la fonction de corrélation des taches se base sur une méthode simple consistant à faire grossir un élément structurant et analyser "ce qui passe au travers". Par exemple, on peut déterminer la taille des pores d'une passoire en testant le passage de grains de sable de différents diamètres. La courbe représentative des grains ayant traversé en fonction de leur diamètre est maximale pour des grains de petite taille, puis diminue car des grains se bloquent mutuellement. Elle ne devient nulle que lorsque les grains ont atteint un diamètre tout juste supérieur à celui des pores de la passoire.

Cet effet est celui qui a été exploité ici. Une cellule activée (noire, par exemple) se voit affectée d'une valeur égale à

, tandis qu'une cellule non activée reçoit la valeur

. On fait alors grossir un carré de taille

(avec

variant de 1 pixel à la moitié de l'image), que l'on déplace à travers toute l'image (il y a N(k)

déplacements qui dépendent de dimx

correspondant aux dimensions de la grille étudiée). A chaque déplacement, on additionne les valeurs des cellules incluses dans cet élément structurant (en valeur absolue et en normalisant par le nombre de cellules : k^2

). On reporte la moyenne de toutes les sommes correspondant à tous les déplacements du carré en fonction de sa taille

(voir Figure ci-dessous).

Exemple de fonction de corrélation

La courbe rouge représente la moyenne de la somme des cellules pour tous les déplacements d'un carré de taille

en fonction de

. La courbe bleue est sa dérivée calculée grossièrement comme étant le taux de variation entre les points consécutifs. On part de petits carrés de taille 1 pixel (le premier point de cette courbe est à 0 mais n'a bien entendu aucune signification car un élément structurant de taille nulle n'est pas défini), cette situation est schématisée en (a). Les cellules incluses dans un carré de 1 pixel de coté sont forcément au nombre d'une unité :

, ce qui, en valeur absolue et une fois normalisé, donne le résultat attendu f(k=0)=1

Plus la taile du carré croît, et plus la somme de son contenu diminue (situation (b)) : les valeurs des cellules activées et non activées commencent à se compenser mais une composante domine toujours l'autre. Arrive alors le moment où la taille du carré, peu importe sa position sur l'image, englobe quasiment autant de cellules activées que de désactivées (situation (c)). La courbe rouge rencontre alors une inflexion relativement bien prononcée car dès lors, la croissance du carré ne change pas fortement la valeur de la fonction de corrélation.

On remarque alors la présence de deux régimes distincts : une forte décroissance de la fonction de corrélation suivie d'une période plus plate où elle se rapproche d'une valeur constante (qu'on peut noter

). Ce genre de courbe peut être approximée par une fonction de type $f(x) = (1-a)e^{(-x/b)} - a$ . Cependant, l'ajustement d'une telle fonction sur ces points expérimentaux pour déterminer les coefficients a et b n'est pas une chose aisée. Pour contourner cette difficulté, deux ajustements linéaires ont été préférés. On détermine le point approximatif séparant les deux ajustements en utilisant la dérivée de la fonction de corrélation (courbe bleue) : on définit ce point approximatif x_c

là où elle s'annule.

Ensuite, on réalise un ajustement linéaire par la méthode des moindes carrés sur les points compris entre x=0

. Sur la deuxième série de points, une moyenne de leurs valeurs est calculée et donne directement l'ordonnée à l'origine d'une fonction affine constante. L'intersection entre ces deux droites détermine la taille caractéristique du système issue de la fonction de corrélation.

Transformée de Fourier

Transformée de Fourier (TF) en échelle logarithmique de "Réseaux imbriqués" avec une maille hexagonale

Section moyennée sur tous les angles de la TF précédente

Pareille aux figures d'interférences lumineuses qui se construisent peu ou prou par coïncidence d'ondes électromagnétiques modulées par un masque, la figure ci-dessus, appelée transformée de FOURIER (ou TF), est une représentation dans l'espace de FOURIER du motif qui se développe. Les zones lumineuses et sombres correspondent à des interférences constructives et destructives rendues possibles grâce aux périodicités de la figure transformée. Ainsi, la transformée de FOURIER nous renseigne sur les façons dont se répètent des motifs caractéristiques qui, bien souvent, échappent à la vue.

Hormis les directions privilégiées introduites par le choix du maillage (cf dernier paragraphe de cette section), on n'observera d'anisotropie claire dans le motif que pour certaines configurations bien spécifiques. Cela étant, on a réalisé un moyennage anneau par anneau de la courbe d'intensité de la TF afin de pouvoir la représenter en fonction de la distance au pixel central*1 (cf le spectre de puissance ci-dessus). Si l'interprétation de la TF nécessite une certaine capacité de lecture dans l'espace de FOURIER, ce genre de section intégrée sur les angles permet d'avoir un aperçu plus rapide des périodicités du motif. Ce spectre de puissance (en échelle logarithmique pour mieux faire ressortir les rebonds successifs appelés harmoniques au-delà du 1er pic - le fondamental) est le résultat d'une telle opération pratiquée sur une TF analogue à celle ci-dessus. On voit clairement apparaître le fondamental (premier anneau éclairé sur la TF) qui peut être précisément localisé dans l'espace de FOURIER, en termes de pulsations spatiales. Ces dernières sont reliées aux dimensions dans l'espace direct (en pixels) par la relation suivante :

$\omega=\frac{2\pi}{l}$

où l correspond à la dimension spatiale associée à cette pulsation. Cette relation indique que plus les dimensions évoquées dans l'espace direct sont grandes, plus les pulsations spatiales correspondantes sont petites. De cette façon, les à-coups et autres discontinuités dans le motif contribueront à augmenter les valeurs de la TF à hautes pulsations alors que les grandes structures contribuent à la TF aux basses pulsations. D'après la définition ci-dessus, L (la taille totale de l'échantillon dont on étudie la TF) fixe donc la résolution de la TF (impossible de voir d'éventuels pics de moins d'un pixel de large sur cette figure!). Ceci s'explique naturellement par l'impossibilité d'observer, dans l'espace direct, la répétition d'un motif plus grand que la zone étudiée*2. Ainsi, cette figure nous permet d'affirmer qu'il existe une fréquence spatiale caractéristique, donnée par le maximum de la figure ie le 1er pic, qui rend compte de la répétition d'un motif bien spécifique de dimension ~ 35 pixels (l'abscisse du fondamental). Comme pour une fonction sinus, il arrive qu'une fonction périodique ait une allure similaire sur 2 demi-périodes successives*3 (ou plus généralement, sur n n-ièmes de période) ce qui "leurre" la TF et crée des maximas locaux dans la figure ci-dessus (au nombre de 2 ici). Ces harmoniques nous renseignent sur la régularité du motif, mais celui-ci peut lui-même comporter des périodicités internes! Bien que la résolution de cette applet aura tôt fait d'entraver la route de celui qui mettra sérieusement cette proposition à l'épreuve, les quelques points qui ponctuent les grandes structures du hanneton foulon ne sont sûrement pas étrangers aux soubresauts que connaît la section aux hautes pulsations.

Deux TF du même motif sans périodicité, mais avec un maillage différent

Les 2 figures juste ci-dessus illustrent l'influence de la maille sur la TF (et donc aussi, bien que de façon moins évidente, dans le motif en lui-même!). Aux grandes pulsations, on observe le comportement du motif aux échelles les plus petites fixées ici par le pixel dont la géométrie définie un maillage. Naturellement, la forme du maillage se fait sentir ainsi que les pulsations qui lui sont inhérentes. Ici, les périodicités ne sont ni spatiales ni temporelles mais angulaires ; en effet, un carré reste inchangé par des rotations de $\pi/2$ (=> 2 directions privilégiées) et un hexagone reste inchangé par des rotations de $\pi/3$ (=> 3 directions privilégiées). Ces propriétés se retrouvent dans la TF où on reconnaît les axes de symétrie de la maille (on laissera au lecteur le loisir de deviner quelle TF est issu de quel maillage…). Parce qu'il comporte plus d'axes de symétrie, le maillage hexagonal rend plus compte d'une modélisation isotrope de la morphogénèse. A haute pulsation, il convient alors de lire les résultats avec précaution afin de faire la part des choses entre les périodicités du signal auquel on s'intéresse et celles de la stratégie d'échantillonnage choisie par l'observateur. Pour reprendre l'analogie exoplanétologique, avec un satellite comme Kepler qui intègre la luminosité sur près d'une demi-heure, malheur à celui qui s'intéresserait à des corps sur une orbite de quelques heures sans prendre garde aux problèmes d'intégration… En revanche, dès que le motif comporte des périodicités de taille bien supérieures à la cellule élémentaire, on observe que, dans la TF, les artefacts dues aux symétries de la maille sont noyés dans le signal issu de ces périodicités (cf la première des deux TF de la configuration guépard des savanes). Pour ce qui est du cas intermédiaire, on peut consulter la deuxième TF ci-après où se mêlent TF des périodicités du motif et de la maille (dans le cas général, on parle de "convolution" entre le maillage, le signal et la fenêtre de mesure).

Deux TF de la configuration guépard des savanes, en maille carrée et hexagonale

N.B. : le premier "F" de "FFT" (pour Fast FOURIER Transform) explique pourquoi les grilles sont carrées (si l'option corps est choisie en entrée, autrement la TF n'est pas calculée) et de taille une puissance de 2. En effet, comme esquissé plus haut, la transformation de FOURIER est une opération qui jouit de propriétés de symétrie très riches et qu'il serait dommage de ne pas exploiter. En "pliant successivement en 2 la TF" (un nombre de fois d'autant plus grand que la taille de la grille est proche d'une puissance de 2), on s'évite des calculs numériquement chronophages en devinant des coefficients par symétrie.

*1 pour les puristes, la TF ci-contre est une translation de la moitié de la dimension du motif dans les deux directions et n'a pas été normalisée (mais la moyenne du motif a été soustraite avant de réaliser la transformation de FOURIER).

*2 de même, l'exoplanétologie suggère une analogie temporelle éloquente : si des satellites comme Kepler ou Corot observent en continu des étoiles en espérant déceler une diminution de leur luminosité due au transit d'une exoplanète, c'est pour disposer d'une base de temps (analogue temporelle de L ci-dessus) suffisante pour détecter la répétition d'un transit d'orbite en orbite (dont la durée - 1 an dans le cas de la Terre - est l'analogue temporelle de la période spatiale de 35 pixels ci-dessus).

*3 on notera que sin(x)

et $sin(2\pi-x)$ ne diffèrent que d'un moins…

Quelques résultats

Robustesse

Pour tester la reproductibilité des mesures, on peut relancer plusieurs fois la simulation pour différents paramètres. Comme le montre la figure ci-dessus, qui représente la distribution de longueurs de corrélation pour plusieurs itérations pour les mêmes paramètres, il existe un effet de discrétisation. Cependant celui-ci peut être plus dû aux effets de calcul de la longueur par ajustement sur la courbe qu'à un effet d'états initiaux aléatoires.

Taille des taches

Comme il se doit, l'applet Java permet de modéliser des taches semblables à celles d'un guépard (mais pas celles d'un léopard, plus complexes!). En témoigne le set Guépard des savanes qui donne, après quelques itérations, un pelage tacheté. Les longueurs de corrélation et de périodicité calculées s'accordent autour d'une dimension caractéristique d'environ 7 pixels pour cette configuration mais il est possible de modifier cette longueur. En considérant les rayons de l'activateur et l'inhibiteur, nous avons esquissé une étude de l'espace des paramètres (voir section suivante). Par exemple, le set de paramètres ci-dessous donne des taches de dimension beaucoup plus faible. On est alors à la limite du critère de SHANNON qui interdit toute détection significative de dimension caractéristique de moins de 2 pixels. La fréquence spatiale associée est donnée par un rayon, pour le cercle dans la TF, grand comme la moitié d'un côté de son cadre (presque atteint). En-deçà de cette limite de 2 pixels, on s'expose à provoquer un repliement de spectre dont la présence brouille la lecture de la TF (bien souvent de façon subreptice pour l'utilisateur novice de la TF!). En conformité avec les remarques de la section précédente, c'est sans inquiétude que l'on note l'ampleur du rôle de la maille dans la TF à des échelles si petites dans l'espace direct : non seulement des lobes latéraux apparaissent dans la TF, et ce en dépit de l'échelle linéaire choisie pour la représenter, mais une plongée à l'oeil nu dans le motif lui-même laisse entrevoir une figure symptomatique d'un pavage hexagonal (cf rosace hexagonale ci-dessous).

Configuration dérivée du Guépard des savanes, dans laquelle on réalise les déformations que peut introduire le maillage

Exploration de l'espace des paramètres

Cette applet Java a avant tout un rôle pédagogique. Nous voulons croire qu'elle participe, à son modeste niveau, à illustrer les concepts parfois délicats qui sous-tendent les automates cellulaires et les équations de réaction-diffusion. Pourtant, un autre support numérique (eg un programme C++) serait l'occasion pour un lecteur curieux d'en savoir plus d'envisager une entreprise délicate pour une applet, une exploration de l'espace des paramètres. L'utilisateur se convaincra facilement de la sensibilité du motif aux paramètres ; il est facile de passer à côté de configurations pertinentes sans s'en apercevoir et difficile de prétendre avoir échantillonné de façon exhaustive sinon significative l'ensemble des configurations permises par ce formalisme (et les humbles capacités de calcul des machines).

Ne résistant pas à l'idée d'amorcer l'exploration d'un nouveau monde, nous nous sommes aventurés à esquisser une étude de corrélation entre l'évolution des rayons (activateur et inhibiteur) et la valeur de la longueur de corrélation associée (après stabilisation). Les résultats sont représentés ci-après :

Longueur de corrélation en fonction des rayons

Si les graphiques ci-dessus (réalisés avec une version plus ancienne de l'applet dans laquelle subsistait un formalisme plus archaïque) s'avèrent jamais avoir une vertu, ce sera bien celle d'avoir suscité la vocation exploratrice chez le chaland enthousiasmé. Une véritable exploration de l'espace des paramètres nécessite (du temps,) un talent et une expérience certaines pour anticiper les régions susceptibles d'être le lieu d'une découverte et appréhender les résultats pertinents.

Le marsupilami, un animal pas comme les autres

Un résultat célèbre obtenu par Alan TURING alors qu'il étudiait la morphégénèse, est illustré avec la configuration "Marsupilami corrigé". Avec son modèle, TURING parvint à expliquer pourquoi, alors que les mêmes mécanismes cellulaires locaux étaient à l'oeuvre sur l'ensemble du corps d'un animal comme le guépard, ce dernier présentait des taches circulaires sur le corps et des bandes orthogonales à la direction privilégiée sur la queue (cf figure ci-dessous). Dans la situation Marsupilami corrigé (comme toujours en géométrie "queue"), on a court-circuité les calculs de la fonction de corrélation et de la TF pour des raisons techniques. L'intérêt de cette configuration réside avant tout dans l'allure du motif en géométrie queue par rapport à celui obtenu, avec les mêmes paramètres, en géométrie corps. Bien que le mécanisme local d'activation des cellules reste le même, tout se passe comme si l'évolution des motifs dépendait de la géométrie globale du système! En effet, les mémes paramètres vont donner des bandes en géométrie queue et des taches en géométrie corps.

Remarquer la différence corps/queue

Crédit : Malene Thyssen (guépard)

L'explication réside dans les conditions aux limites imposées au système. Ici, le terme queue est abusif en ce sens que, si la géométrie corps modélise bel et bien une sphère ( ~ corps réel de guépard ou de dalmatien), la géométrie queue correspond en fait à un tore (conditions au limites verticales mais aussi horizontales - le début et la fin de la queue sont reliés!). Néanmoins, la configuration "Marsupilami corrigé" a été réglée pour être sensible aux effets de bords haut/bas mais pas droite/gauche : on modélise bien la queue cylindrique d'un animal et on retrouve ainsi le célèbre résultat d'Alan TURING.

Motif dérivé de la configuration marsupilami (rayon de l'inhibiteur plus petit)

En diminuant le rayon de l'inhibiteur, on finit par perdre cette périodicité (les cellules ne se voient plus d'un bord à l'autre) et des taches réapparaissent sur la queue (cf queue ci-dessus). Ainsi, l'observation de la transition taches/bandes sur le pelage d'un animal ou dans un paysage renseigne sur le rayon de diffusion de l'entité qualifiée d'activatrice.

N.B. : l'amateur de bande dessiné n'aura pas manqué de remarquer que la queue du marsupilami original est tacheté, pas rayée, en dépit de taches similaires à celles du guépard sur le corps. Un bien étrange animal aux rayons fluctuants...

Animaux

Ci-dessous, on peut voir quelques exemples de pelage ou carapaces d'animaux dont la structure est proche de certaines simulations. Les images seront volontairement prises grandes (configuration la plus puissante) pour le lecteur malchanceux dont la machine ne permettrait pas l'affichage de ces grandes configurations.

Guépard (Acinonyx Jubatus)

Les paramètres sont : R1=2, R2=10, J1=1, J2=-0.1, h=-18.

Guépard royal (Acinonyx jubatus rex)

Les paramètres sont : R1=5, R2=10, J1=1, J2=-0.1, h=0.

Hanneton Foulon (Polyphylla fullo)

Les paramètres sont : R1=3, R2=4, J1=6.5, J2=-3, h=0.

Cet animal, encore appelé hanneton des pins, prouve que les mammifères ne sont pas forcément les seuls concernés par le modèle de Turing.

Poisson-chirurgien d'Hawaï

Ci-dessus, on peut observer l'apparence de l'ensemble ichtyen pour une densité faible de cellules de Voronoy. Cela permet d'observer l'agencement désordonné de la grille. Ci-dessous, le même ensemble est généré avec la densité la plus élevée prévue. On voit que l'apparence est proche de celle des autres géométries (et du poisson) avec une meilleure isotropie. Cela se vérifie sur la transformée de Fourier correspondante.

En observant la photographie du poisson, on remarque d'ailleurs une différence entre les motifs de sa tête, ses nageoires et son corps. La simulation est plus proche d'une géométrie tête (avec ses conditions aux limites périodiques), d'où la plus forte ressemblance avec la tête de l'animal.

Dunes

Dunes du désert Namib en Namibie

Les paramètres sont : R1=2, R2=14, J1=0.5, J2=-0.1, h=0.

Certains géophysiciens ont appliqué le modèle de Turing avec succès à la formation de dunes (voir ici).

Méthode de Turing et morphogénèse : illustration des équations de réaction-diffusion

Table des matières

Liste des paramètres d'entrée de l'applet

Liste des paramètres de sortie de l'applet

Introduction

Champ

Objectifs