L’oscillateur de Van der Pol

Au début des années 1920, le physicien néerlandais Balthazar van der Pol expérimentait avec des tubes à vide d’un type particulier, appelés triodes. L’intérêt de ce genre de tubes à vide est qu’ils peuvent servir d’amplificateurs dans un circuit électronique, par exemple dans un poste récepteur de radio, ou encore de portes logiques. C’est pourquoi, avant l’invention en 1947 des transistors à semi-conducteurs, les ordinateurs étaient construits avec des tubes à vide.

Le circuit électrique suivant est un modèle (légèrement simplifié) d’une triode, appelé oscillateur de van der Pol. Ses éléments de base sont une résistance, de valeur $R$, une bobine d’inductance $L$, un condensateur de capacité $C$ ; l’élément $r(v)$ désigne également une résistance, mais avec la particularité que sa valeur dépend de la différence de potentiel $v$ entre ses extrémités. Enfin, l’élément marqué $u$ est une source de tension.

Le potentiel électrique $v$ à l’endroit indiqué sur la figure évolue au cours du temps $t$ selon l’équation différentielle du second ordre
\[ \frac{\text{d}^2 v}{\text{d} t^2} -\frac{1}{RC} \frac{\text{d}v}{\text{d}t} + \frac{1}{C} \frac{\text{d}}{\text{d}t} \left( \frac{v}{r(v)} \right) + \frac{1}{LC} v = 0 \]

Si $r(v)$ était constante, on aurait affaire à une version de circuit RLC. L’équation différentielle serait alors linéaire, et on pourrait écrire sa solution en termes de fonctions sinus, cosinus et exponentielle. Les circuits RLC servent d’oscillateurs, et sont utilisés par exemple comme horloge d’un ordinateur ou pour la génération de sons dans un synthétiseur.

Dans le cas de la triode, en revanche, la résistance $r(v)$ décroît lorsque $v$ augmente. On suppose d’habitude qu’elle peut s’écrire sous la forme $r(v)=RV_0^2/v^2$ pour une constante $V_0$ (mesurée en Volts). Dans ce cas, $v/r(v)=v^3/(RV_0^2)$, dont la dérivée vaut $(3v^2/(RV_0^2)) \text{d}v/\text{d}t$. La variable sans dimension $x=\sqrt{3}v/V_0$ obéit alors à l’équation différentielle
\[ \frac{\text{d}^2 x}{\text{d} t^2} -\frac{1}{RC} (1-x^2) \frac{\text{d}x}{\text{d}t} + \frac{1}{LC} x = 0 \]
qui est appelée équation de van der Pol [vdP1]. En raison du terme en $x^2$, ce n’est plus une équation linéaire, et il n’existe pas de formule explicite pour ses solutions.

Pour étudier sa dynamique, on peut cependant recourir à une méthode purement géométrique, dans la tradition d’Henri Poincaré. Pour faire cela, on commence par réécrire l’équation du second ordre sous la forme d’un système de deux équations du premier ordre. Ce procédé n’est pas unique, mais une possibilité parmi d’autres est de considérer le système
\[ \begin{cases} \displaystyle \frac{\text{d}x}{\text{d}t} &=& x - \displaystyle \frac13 x^3 - y \\ \displaystyle \frac{\text{d}y}{\text{d}t} &=& \displaystyle \frac1T x \end{cases} \]
où le paramètre $T$ est égal à $L/(R^2C)$.

On peut déjà obtenir beaucoup d’informations utiles, simplement en regardant les signes des membres de droite du système. En effet, si $y > x - \frac13 x^3$, alors $\text{d}x/\text{d}t$ est négatif et donc $x$ diminue, alors que si $y < x - \frac13 x^3$, $x$ augmente. De même, puisque $T>0$, la variable $y$ augmente si $x>0$ et diminue si $x<0$. Le plan $(x,y)$ peut donc être subdivisé en $4$ parties, comme dans la figure ci-dessous. Les flèches indiquent le sens de variation des coordonnées dans chaque domaine. Par exemple dans la zone bleue, $x$ et $y$ diminuent. Dans la zone verte, $x$ augmente et $y$ diminue, et ainsi de suite. Ces domaines sont délimités par la courbe d’équation $y = x - \frac13 x^3$ et la droite d’équation $x=0$, qu’on appelle les nullclines du système [1].

Comme le suggèrent les flèches, les solutions de l’équation vont avoir tendance à tourner dans le sens contraire des aiguilles d’une montre. La simulation interactive ci-dessous vous permettra d’explorer la dynamique, en dessinant la solution partant d’un point donné. Cliquez d’abord sur « Lancer », puis au point choisi comme valeur initiale. La solution partant de ce point va se dessiner dans le plan $(x,y)$, et la valeur de $x$ en fonction du temps sera affichée en dessous. Vous pouvez modifier le paramètre $T$ au cours de la simulation ($T$ est donné par $10^S$ où $S$ est la position du curseur ; ainsi $T=0.1$ lorsque le curseur est à gauche, $T=1$ lorsqu’il est au centre, et $T=10$ lorsque le curseur est à droite).

On s’aperçoit que les solutions tendent toutes à s’approcher d’une courbe fermée, qu’on appelle un cycle limite, ou encore une orbite périodique attractive. Pour de petites valeurs de $T$, les courbes $x(t)$ ressemblent aux oscillations sinusoïdales d’un circuit RLC linéaire, comme la courbe bleue ci-dessous. En revanche, lorsque $T$ augmente, le cycle limite a tendance à s’aplatir, avec des parties presque horizontales parcourues rapidement, et des parties obliques parcourues plus lentement. La courbe $x(t)$ est alors constituée de parties presque verticales, et de parties à la pente plus douce (courbe rouge ci-dessous). Ce comportement, très différent des oscillations sinusoïdales d’un circuit RLC linéaire, a été appelé oscillations de relaxation par van der Pol [vdP2].

Il est intéressant de noter que Van der Pol, avec son collègue Van der Mark, ont proposé que ces oscillations pouvaient servir de modèle aux battements du cœur [vdPvdM].

Neurones

Notre système nerveux est constitué, comme chacun sait, de cellules appelées neurones. La figure ci-dessous donne une représentation schématique d’un neurone.

Tout neurone est connecté à un nombre souvent important d’autres neurones, par l’intermédiaire des extrémités appelées dendrites, ainsi que par un axone. L’axone peut être très long [2], et transmet des impulsions électrochimiques, appelées potentiels d’action. Voici un exemple de mesure de potentiels d’action ; plus précisément, il s’agit de la différence de potentiel électrique à travers la membrane, mesurée à plusieurs reprises en un point de l’axone pendant une dizaines de secondes [BELW].

Comme on le voit, le potentiel montre des pics espacés irrégulièrement, qu’on appelle des spikes. Ces spikes ont toujours à peu près la même forme, schématisée ci-dessous. L’information transmise entre neurones est donc contenue dans les intervalles de temps entre spikes : des spikes plus rapprochés dans le temps traduisent un signal plus intense.

En 1952, Alan Lloyd Hodgkin et Andrew Huxley mirent au point un modèle de circuit électrique, reproduisant très fidèlement le potentiel d’action de nombreux types de neurones [HH]. Ce travail leur valut le Prix Nobel de médecine en 1963, partagé avec Sir John Carew Eccles.

Le modèle de Hodgkin et Huxley représente la membrane de l’axone par un circuit électrique de capacité $C$, ainsi que par trois types de courant associés à des canaux ioniques. Ces canaux sont de petits orifices dans la membrane, qui laissent passer plus ou moins de ions sodium (Na) ou potassium (K) selon la valeur du potentiel électrique $v$ [3]. On a coutume de spécifier les conductances $g_{\text{Na}}$, $g_{\text{K}}$ et $g_{\text{L}}$ plutôt que les résistances, qui sont simplement leurs inverses. On trouvera par exemple ici quelques animations illustrant le fonctionnement de ces canaux.

Le modèle de Hodgkin—Huxley donne lieu à un système de $4$ équations différentielles. Ses variables sont le potentiel de membrane $v$, ainsi que les proportions de canaux ioniques ouverts de trois types différents [4]. Nous n’écrirons pas ici ces équations, assez compliquées, notamment parce que la dépendance des conductances $g_{\text{Na}}$ et $g_{\text{K}}$ du potentiel $v$ fait intervenir des fonctions non polynomiales.

Le modèle de FitzHugh—Nagumo

Étudier mathématiquement un système de $4$ équations différentielles est difficile, même avec des outils géométriques. Un tel système peut par exemple présenter des solutions chaotiques [5]. Pour cette raison, plusieurs modèles simplifiés ont été imaginés pour décrire les potentiels d’action, le plus connu étant probablement le modèle de FitzHugh—Nagumo.

Ce modèle fut proposé en 1961 par Richard FitzHugh [F] [6], et Jin-ichi Nagumo l’interpréta comme un circuit électrique en 1962 [NAY]. Le système est une généralisation de l’oscillateur de van der Pol, et s’écrit
\[ \begin{cases} \displaystyle \frac{\text{d}v}{\text{d}t} &=& v - \displaystyle \frac13 v^3 - w + I\\ \displaystyle \frac{\text{d}w}{\text{d}t} &=& \displaystyle \frac1T (v+a-bw) \end{cases} \]
La variable $v$ représente le potentiel de membrane, alors que $w$, appelée variable de récupération, décrit l’effet combiné des canaux ioniques. Les grandeurs $I$ (qui représente un courant de stimulation, provenant par exemple des autres neurones), $a$ et $b$ sont des paramètres du modèle. Remarquez que si $I=a=b=0$, on retrouve le système de Van der Pol.

La simulation ci-dessous vous permet d’expérimenter l’effet des différents paramètres sur le modèle. Avant de lire la suite, prenez le temps d’explorer les dynamiques possibles des équations.

La caractéristique la plus importante du modèle est que selon la position relative des nullclines, il présente deux types qualitativement différents de dynamiques :

• Soit toutes les trajectoires convergent vers un point d’intersection des nullclines. On dit que ce point est un point d’équilibre stable, et qu’il forme un attracteur global du système.
• Soit les trajectoires s’approchent d’un cycle limite, et le système présente des oscillations comme l’oscillateur de van der Pol ; pour des valeurs convenables des paramètres $I$, $a$, $b$ et $T$, ces oscillations reproduisent la forme des spikes observés dans des neurones.

Plusieurs trajectoires typiques sont représentées ci-dessous, dans un cas où le point d’intersection des nullclines est stable. Pour des valeurs initiales en dessous d’une courbe seuil (notée « quasi-threshold » dans la figure qui suit), on obtient des spikes, dont les différentes phases ont des interprétations physiologiques précises (qu’il n’est toutefois pas nécessaire de comprendre pour lire la suite).

Une seconde caractéristique du modèle est l’excitabilité. Elle signifie que lorsque le système est dans son état stable, mais proche du seuil d’apparition d’oscillations, de petites perturbations suffisent à produire des spikes. C’est une caractéristique commune à tous les modèles de neurones.

Bien qu’il soit une simplification du modèle de Hodgkin—Huxley, le modèle de FitzHugh—Nagumo permet de reproduire un certain nombre de comportements observés. Par exemple, si l’on fait varier la variable $I$, le neurone peut passer de l’état au repos à l’état actif, avec émission de spikes, puis retourner à l’état au repos. On trouvera quelques animations de ces transitions sur cette page.

Les canards

Une question importante du point de vue des mathématiciens est : Comment passe-t-on d’un point d’équilibre stable à un cycle limite, lorsqu’on varie un ou plusieurs paramètres ? La théorie générale des équations différentielles montre que cette transition doit être continue : on ne peut pas sauter d’un point d’équilibre pour une valeur des paramètres à un grand cycle limite pour une valeur infiniment proche des paramètres.

La transition entre équilibre stable et cycle limite s’effectue par ce qu’on appelle une bifurcation de Poincaré—Andronov—Hopf. En modifiant petit à petit un paramètre, le point d’équilibre devient instable, et il se crée un petit cycle limite, qui grandit progressivement.

L’animation ci-dessous montre ce qu’il se passe une fois le point d’équilibre devenu instable. Pour $I=b=0$ et $T=10$, on fait varier petit à petit la valeur de $a$. Pour chaque valeur de $a$, l’animation montre une trajectoire du système, toujours issue du même point (situé en bas à droite de l’image).

Ce qui est frappant, c’est que pour des valeurs de $a$ voisines de $0.9875$, le cycle limite grandit très vite. Sa forme est extrêmement sensible à la valeur du paramètre $a$. Plus $T$ est grand, plus ce phénomène est marqué. Le logo de l’article donne une autre illustration du phénomène : on y a superposé les cycles limites obtenus pour différentes valeurs de $a$, espacées de $0.00001$.

Lorsque $a$ est proche de sa valeur critique, le cycle limite épouse la partie centrale de la nullcline $w = v - v^3/3$. C’est la région appelée « no man’s land » dans la figure précedente. Dans cette région, la dynamique est extrêmement sensible aux conditions initiales : une petite perturbation peut faire dévier la solution vers la gauche ou vers la droite.

Ce phénomène surprenant a été découvert au début des années 1980 par les mathématiciens français Eric Benoît, Jean-Louis Callot, Francine Diener et Marc Diener [BCDD]. Ils ont appelé les solutions particulières qui restent proches de la partie centrale de la nullcline des canards. Pourquoi canards ? La figure suivante en fournit une explication : lorsque $T$ est très grand, le cycle limite tend vers une courbe formée de segments horizontaux et de bouts de nullcline, dont la forme évoque celle d’un canard.

L’étude mathématique des canards présente de nombreuses difficultés, et a donné lieu à un grand nombre de travaux depuis leur découverte. On aimerait par exemple, lorsque $T$ est très grand, obtenir une expression approchée de la solution sous forme d’un développement limité en $1/T$. Comme la forme limite, lorsque $T$ tend vers l’infini, n’est pas lisse du tout [7], il faut être astucieux pour obtenir un tel développement.

Quel est le rôle des canards dans la dynamique des neurones ? Nous avons vu que la plupart du temps, soit le neurone est au repos, soit il est actif et émet des spikes. Pour les valeurs exceptionnelles des paramètres où il existe un canard ($a=0.9875$ dans l’animation), la dynamique du neurone est intermédiaire : le potentiel oscille périodiquement, mais sans parvenir à produire de vrais spikes. Lorsque $T$ est grand, les canards sont très difficiles à attraper, car ils n’existent que dans un tout petit intervalle de paramètres [8]. Nous le voyons dans l’animation : une modification de $a$ de $\pm 0.0001$ suffit à faire basculer le système soit dans l’état au repos, soit dans l’état actif.

Pourquoi alors s’intéresser aux canards ? C’est précisément parce qu’ils servent de frontière entre deux comportements qualitativement différents qu’il est important de savoir les localiser. Si l’on sait trouver les canards, alors on sait prédire pour quelles valeurs des paramètres le neurone sera soit au repos, soit actif. [9]

Dès que l’on considère des systèmes à $3$ variables au lieu de $2$, les canards deviennent nettement moins élusifs : Dans de nombreux cas, les canards existent dans un large domaine de l’espace des paramètres [B]. Au lieu de présenter seulement deux comportements possibles (actif et au repos), les systèmes de dimension $3$ ou plus peuvent faire coexister plusieurs régimes actifs, avec des modes d’oscillation différents, en plus d’un état au repos. Selon les valeurs initiales des variables, on observe alors l’un ou l’autre comportement, et les canards servent à nouveau de frontière entre ces dynamiques. Ce phénomène, qui est d’origine géométrique, fera peut-être l’objet d’un prochain article sur Images des Maths.

Bibliographie

[B]
E.Benoît. Canards et enlacements, Publications de l’Institut des Hautes Études Scientifiques, 72:63—91, 1990.

[BCDD]
E. Benoît, J.-L. Callot, F. Diener, et M. Diener, Chasse au canard, Collect. Math., 31:37—119, 1981.

[BELW]
B.G. Burton, M.N. Economo, J. Lee, et J.A. White, Development of theta rhythmicity in entorhinal stellate cells of the juvenile rat, Journal of Neurophysiology 100 : 3144—3157, 2008.

[F]
R. FitzHugh, Impulses and physiological states in models of nerve membrane, Biophysical Journal, 1, 445—466, 1961.

[HH]
A. L. Hodgkin et A. F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, The Journal of physiology 117 (4) : 500—544, 1952.

[NAY]
J. Nagumo, S. Arimoto, et S. Yoshizawa, An active pulse transmission line simulating nerve axon, Proceedings of the Institute of Radio Engineers, 50, 2061-2070, 1962.

[vdP1]
B. van der Pol, A theory of the amplitude of free and forced triode vibrations, Radio Review, 1, 701-710, 754—762, 1920.

[vdP2]
B. van der Pol, On relaxation-oscillations, The London, Edinburgh and Dublin Phil. Mag. & J. of Sci., 2(7), 978—992 (1926).

[vdPvdM]
B. van der Pol et J. van der Mark, The heartbeat considered as a relaxation oscillation, and an electrical model of the heart. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science Ser.7, 6, 763-775, 1928.