# Un exemple classique de modèle mathématique en neurosciences

Christophe Pouzat

► To cite this version:

Christophe Pouzat. Un exemple classique de modèle mathématique en neurosciences. 2025. hal-05054356v2

HAL Id: hal-05054356

https://hal.science/hal-05054356v2

Preprint submitted on 11 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 - Public Domain Dedication 4.0 International License
# Un exemple classique de modèle mathématique en neurosciences

Christophe Pouzat

IRMA CNRS UMR 7501

Université de Strasbourg

christophe.pouzat@math.unistra.fr

11 juin 2025

Les neurosciences présentent une vaste collection de questions et problèmes d'échelles — spatiale et temporelle — et de nature — déterministe et stochastique — variées. Une classe de modèles, que je qualifierai de réductioniste, rend compte de données expérimentales en décomposant le phénomène observé en entités partiellement ou complètement indépendantes; cette décomposition reflète ou suggère de nouvelles expériences et permet de faire progresser notre compréhension en cherchant un accord quantitatif entre données expérimentales et prédictions théoriques. Une autre classe de modèles, que je qualifierai de minimaliste-généraliste, rend compte de phénomènes observés dans divers contextes expérimentaux, — comme le présence d'un régime métastable, c'est l'aspect « général » — en sélectionnant parmi les propriétés neuronales connues, le plus petit nombre possible — c'est l'aspect « minimaliste » — permettant au phénomène d'intérêt de se manifester / d'être reproduit. Ces deux classes doivent être considérées comme les deux extrêmes entre lesquels les modèles mathématiques « utiles » en neurosciences peuvent être situés. On l'aura compris, la littérature est vaste et ce n'est dans les quelques pages qui suivent que je vais pouvoir lui faire justice. Je vais ainsi me concentrer sur un sujet très « classique » : la génération du potentiel d'action et sa propagation. Nous allons donc discuter d'un phénomène à l'échelle du
neurone, c'est-à-dire avec une des trois dimensions spatiales allant de la *centaine de microns à quelques mètres*, avec un temps caractéristique de l'ordre de la *milliseconde*. Je vais essentiellement présenter l'aspect *déterministe* du problème, en explorant des modèles réductionistes de complexités croissantes, rarement présentés dans la littérature.

# 1 Introduction

Parallèlement à la typologie sommaire des modèles présentée dans le résumé, nous voyons mis en oeuvre un large éventail d'« outils » mathématiques : des équations différentielles ordinaires (EDO) comme des équations aux dérivées partielles (EDP) [26, 20, 22] — c'est ce dont nous allons principalement discuter dans cet article —; des modèles markoviens, de diffusion ou d'intensité (stochastique) [3, 15]; des modèles de graphes [28].

Mais avant de nous plonger dans les aspects mathématiques, il est important que nous ayons une description de la réalité dont nous essayons de rendre compte. Nous allons discuter de *neurones* qui sont des *cellules particulières*. En effet, si la plupart des cellules de notre organisme ont une taille de 10 à 20 µm, les neurones (et les cellules musculaires) peuvent être très grands<fnref n="1" /> [5, figure 4.1E]. Les neurones sont classés en différents types suivant leur morphologie, leur physiologie et la région du système nerveux dans laquelle ils se trouvent. La figure 1, empruntée à Santiago Ramón y Cajal<fnref n="2" /> [4, figure 24, p 116], va nous permettre de présenter certaines propriétés des neurones, mais avant de discuter de celles-ci, il est important — pour le mathématicien intéressé à dialoguer avec des biologistes — de comprendre que cette figure est un exemple de *modèle* comme les biologistes en manipulent tous les jours. La figure 1 est un modèle, comme le sont un plan de métro ou une carte routière. On y distingue des « éléments remarquables », les neurones, qui sont de différents types : dans la muqueuse olfactive (A), on voit des neurones sensoriels; les cellules mitrales (C) et les grains (D) sont les deux principaux types neuronaux du *bulbe olfactif* (le premier relais du système olfactif); les cellules pyramidales (F) sont le principal type neuronal du *cortex*<fnref n="3" />. Ces neurones forment des connexions, des *synapses*, bien spécifiques : les neurones

<footnote n="1">1. Nous parlons de Biologie, il y a donc des exceptions comme les neurones intrinsèques de la rétine.</footnote>

<footnote n="2">2. Contemporain de Poincaré et Hilbert et leur « équivalent » en neurobiologie.</footnote>

<footnote n="3">3. Le cortex, mot latin qui signifie écorce, est la fine couche externe qui enveloppe le cerveau; c'est là que se trouve la plupart des neurones.</footnote>
récepteurs entrent en contact avec les cellules mitrales (C) dans les glomé-
rules (B); les cellules mitrales entrent en contact avec les grains (D) et les
cellules pyramidales (F); les flèches indiquent la direction du « flux de l'in-
formation » (pour employer un terme contemporain). On voit ainsi qu'une
abstraction sous forme de graphe (collection de points et d'arêtes) devrait im-
pliquer différents types de sommets et devrait être dirigée. Au-delà de l'abs-
traction en graphe, la figure illustre une propriété « banale » des neurones qui
sont des cellules asymétriques, chez qui on distingue typiquement trois par-
ties (voir la partie encadrée en bas à droite de la figure) : le corps cellulaire,
structure dont le diamètre est le plus grand (c'est là que se trouve le noyau);
les dendrites <fnref n="4" />, structures arborescentes (orientées vers le bas depuis le corps
cellulaire sur la figure) et l'axone, structure la plus fine (diamètre typique-
ment inférieur au micromètre chez les vertébrés), pouvant être très étendue
(plusieurs mètres pour les moto-neurones de girafes ou de diplodocus) <fnref n="5" />. Ca-
jal utilise cette figure pour illustrer le point suivant : si nous comparons le
flux de l'information à l'orientation des cellules, nous constatons que les den-
drites sont du côté où arrive l'information et que l'axone est du côté où elle va.
Dans le chapitre « Inductions physiologiques tirées de la morphologie et des
connexions des neurones », dont fait partie cette figure, il étudie de façon sys-
tématique les systèmes sensoriels et moteurs et constate que « du point de vue
du neurone », l'information arrive (ou est collectée) par les dendrites et part
(ou est distribuée) par l'axone, ce qu'il résume par la notion de polarisation
dynamique [4, chap. V, pp 110-137].

Dans les quelques pages suivantes, nous allons d'abord présenter de façon qualitative le phénomène que nous souhaitons modéliser : une impulsion de forme fixe, se propageant le long de l'axone à vitesse constante, dès lors qu'une stimulation suffisamment forte est appliquée. Nous introduirons ensuite un modèle quantitatif de l'axone qui nous mènera à une version de l'équation de la chaleur. Ce cadre nous permettra de situer et de décrire précisément les contributions de Hodgkin et Huxley avant d'introduire le modèle d'EDP qui porte leurs noms. Comme ce modèle est trop compliqué pour être résolu analytiquement, nous commencerons par suivre les pas de Zeldovich et Frank-Kamenetsky qui se concentrent, au prix d'un simplification, sur la propagation du « front » de l'impulsion. L'analyse du portrait de phase de ce modèle simple, nous permettra d'établir l'existence d'une unique solution que nous déterminerons d'abord numériquement par la méthode du tir, puis analyti-

<footnote n="4">4. Terme dérivé du mot grec pour « arbre ».</footnote>

<footnote n="5">5. Tous les neurones n'ont pas d'axone; les grains (D) sur la figure n'en ont pas, tout comme la plupart des types cellulaires de la rétine.</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Schéma du circuit olfactif par Santiago Ramón y Cajal (1909), modifié par l'ajout de la partie encadrée en bas à droite.</figcaption></figure>
quement, en restreignant quelque peu le cadre de notre question. Nous enta-
merons alors une « remontée » vers des modèles plus réalistes avec Markine et
Tchizmadjev. Leur modèle, bien qu’analytiquement très simple, admet deux
solutions ayant une forme d’impulsion, une lente et instable, une rapide et
stable. Nous discuterons alors de la simplification la plus connue du modèle
de Hodgkin et Huxley, due à Fitzhugh et Nagumo, et nous mentionnerons les
versions stochastiques récentes avant de conclure.

2 Le Potentiel d'action

2.1 Propriétés expérimentales

Nous venons de voir que l'axone — un cylindre de petit diamètre (< µm chez les vertébrés) et potentiellement très long — constitue la « voie de sortie » de la plupart des neurones. Si nous nous mettons dans la peau d'un neurophysiologiste, nous devons alors essayer de répondre à la question suivante : comment un signal peut-il se propager du corps cellulaire jusqu'à l'extrémité de l'axone de façon fiable et rapide ? Une signalisation basée sur la diffusion de molécules est omniprésente en biologie, mais beaucoup trop lente pour constituer une réponse sérieuse à notre question [6]. Mais les neurones, comme toute cellule, présentent une différence de potentiel électrique entre leur intérieur et leur extérieur ; le premier est typiquement inférieur de 60 à 70 millivolts (mV) au second. Cette différence de potentiel est appelée *potentiel de membrane*. Cela a motivé très tôt les physiologistes à mesurer les changements du potentiel de membrane au cours du temps en différents points de l'axone [27]. Ils ont alors observé le phénomène illustré sur la figure 2 — pour éviter des ennuis de *copyright*, j'utilise des simulations numériques du modèle de Hodgkin et Huxley [16] à la place de véritables données expérimentales, mais 70 ans d'expériences nous montrent que ce modèle constitue une excellente première approximation et reproduit toutes les caractéristiques qualitatives d'un véritable axone. Cette figure montre les « réponses » — changements du potentiel de membrane, la valeur de « repos » étant conventionnellement fixée à 0 —, en différentes positions le long de l'axone, à une « stimulation »<fnref n="6" /> appliquée en une seule position, sachant qu'il est parfaitement possible de stimuler et d'enregistrer une réponse en une même position. La partie supérieure de la figure montre deux réponses (paires de trajectoires noires d'une part et grise d'autre part) à deux stimulations qu'un physiolo-

<footnote n="6">6. Injection d'un courant qui tend à diminuer la différence de potentiel entre l'intérieur et l'extérieur du neurone.</footnote>
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Les caractéristiques marquantes du potentiel d’action. Solutions numériques du modèle proposé par Hodgkin et Huxley [16] pour « l’axone géant » de calmar. Adaptation de [7, figure 2].</figcaption></figure>

giste qualifierait de « sous le seuil ». La réponse à chacune des deux stimulations est montrée en deux positions distinctes : au site d'injection (les deux trajectoires, noires et grises qui croissent dès que le temps est positif et dont les valeurs maximales sont au-dessus de 15 mV) et 1 centimètre plus loin (les deux autres trajectoires). L'intensité du courant injecté croit des trajectoires noires aux grises. Le phénomène remarquable illustré ici est la très forte aténuation de l'amplitude du changement de potentiel de membrane avec la distance : un changement de 15 mV au site d'injection ne donne qu'une déviation de l'ordre de 1 mV à 1 cm. Si j'avais rajouté les trajectoires du potentiel de membrane à 2 cm du site d'injection, on ne verrait rien à l'échelle de cette figure. Sachant qu'un calmar fait plusieurs dizaines de centimètres (voir plusieurs mètres) et que l'axone géant est employé pour déclencher une réponse de fuite très importante pour la survie de l'animal, on voit que ce type de signalisation ne peut « faire l'affaire ». La partie inférieure de la figure illustre ce
qui se passe lorsque l'intensité de la stimulation est augmentée (remarquez le changement d'échelle des ordonnées entre les deux parties de la figure). Ici les trajectoires du potentiel de membrane sont illustrées en 5 positions distinctes, de gauche à droite : au site de stimulation, puis à 1, 2, 3 et 4 centimètres de celui-ci. La trajectoire grise est la même que celle de la partie supérieure de la figure au site de stimulation. Les couleurs bleue et orange correspondent à deux intensités de stimulations différentes. On voit qu'un phénomène tout à fait nouveau apparaît : une impulsion est générée au site de stimulation et se propage à vitesse constante, sans déformation le long de l'axone; c'est ce que les physiologistes nomment potentiel d'action. On voit également, en compara-rant les trajectoires bleue et orange, que si l'intensité de stimulation est assez grande (« au-dessus du seuil », comme dirait un physiologiste), c'est le temps auquel l'impulsion apparaît qui change avec l'intensité, mais ni sa forme, ni sa vitesse ne sont modifiées. Le potentiel d'action semble ainsi être un bon candidat pour une transmission fiable et rapide sur de « longues » distances.

### 2.1.1 Pourquoi le calmar ?

Le calmar, du fait de son axone géant (son diamètre peut aller jusqu'à 1 mm) et de sa synapse géante d'où est issu l'axone précité, a été un animal très employé par les neurophysiologistes des années 30 aux années 80 du siècle dernier — la grande taille des neurones permettant des enregistrements de haute précision avec des technologies relativement rudimentaires<fnref n="7" />. Toute notre compréhension de base de la biophysique du potentiel d'action et de la transmission synaptique vient de cet animal; pensez y la prochaine fois que vous aurez du calmar dans votre assiette !

## 2.2 Le modèle de Hodgkin et Huxley

Nous décrivons maintenant un des plus spectaculaire succès de la physio- logie du XXe siècle, le modèle proposé par Hodgkin et Huxley en 1952 [16]. Ce modèle est construit sur un circuit (électrique) équivalent représentant un « petit morceau » d'axone, couplé à une description quantitative de la dépen- dance vis-à-vis du potentiel de membrane et vis-à-vis du temps des conduc- tances spécifiques de deux ions : les ions potassium (K⁺), dont la concentration est élevée à l'intérieur (397 milli-mole par litre, mM) et faible à l'extérieur de

<footnote n="7">7. Une recherche sur YouTube avec les mots clés *squid* et *Hodgkin* donne accès à de courtes vidéos qui permettent de se faire une idée très concrète de ce type d'expériences; voir en particulier la vidéo « *The squid and its giant nerve fiber* ».</footnote>
l'axone (20 mM); les ions sodium (Na⁺), de concentration faible à l'intérieur (50 mM) et élevée à l'extérieur de l'axone (437 mM).

### 2.2.1 « Circuit équivalent » de l'axone

Bien que l'axone soit physiquement une structure cylindrique, une approximation unidimensionnelle s'avère appropriée — on peut parfaitement négliger les différences de potentiel électrique en deux points intérieurs de l'axone situés sur un même rayon —. La membrane qui sépare l'intérieur de l'extérieur de l'axone est une bi-couche lipidique (comme toutes les membranes cellulaires) qui, lorsqu'elle est « au repos », c'est-à-dire en l'absence de stimulation, est faiblement perméable aux ions potassium (K⁺) et chlorure (Cl⁻) du fait de la présence de canaux ioniques et de « pompes » membranaires <fnref n="8" /> (sans ces pompes et canaux, la membrane serait essentiellement imperméable) [30, Chap. 11]. Des expériences effectuées dans les années 30 et 40 [17, 27] on montré qu'un « petit » morceau d'axone de longueur Δx et de rayon a, peut être modélisé, du point de vue électrique, par une capacité cΔx (la membrane est un isolant) en parallèle avec une conductance gΔx (du fait des canaux et des pompes). Deux petits morceaux voisins sont connectés par une résistance longitudinales rΔx (le milieu intérieur se comporte comme un milieu purement résistif dans le domaine de fréquence qui nous intéresse ici) et la résistance du milieu extra-cellulaire peut être négligée <fnref n="9" />. Ceci nous amène au modèle présenté en figure 3.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Le circuit équivalent d'un axone</figcaption></figure>

Sur cette figure, quatre « petits morceaux » d'axone de longueur Δx [cm] <fnref n="10" />

<footnote n="8">8. Ces connaissances sur l'organisation moléculaires de la membrane ont été acquises après les années 60, ce qui signifie qu'Hodgkin et Huxley ne les avaient pas; ils parlent ainsi de « conductances » ou de « perméabilité » et non de canaux ioniques.</footnote>

<footnote n="9">9. Résistance et conductance sont les deux faces d'une même pièce : l'une est l'inverse de l'autre; la littérature de neurophysiologie a tendance à employer « conductance » pour la membrane et « résistance » pour le milieu intracellulaire. Cela permet accessoirement d'éviter l'emploi d'indices pour distinguer les deux résistances ou les deux conductances, ce qui donne des équations plus compactes.</footnote>

<footnote n="10">10. Lorsqu'une quantité, numérique ou symbolique, est suivie d'une indication entre « [] »,</footnote>
sont représentés, le deuxième apparaît en rouge. Les quantités g, c et r dépendent du rayon a de l'axone et sont précisées dans la liste suivante :

* g [S cm⁻¹] (S est le symbole du siemens, l'unité de conductance du système international) est la conductance de la membrane d'un morceaux d'axone de longueur unité et de rayon a; elle est dérivée de la conductance spécifique d'une surface de membrane unité, G [S cm⁻²], en multipliant cette dernière par la surface d'un morceau de longueur unité : $g = G2\pi a$;

* c [F cm⁻¹] (F est le symbole du farad, l'unité de capacité du système international) est la capacité de la membrane d'un morceaux d'axone de longueur unité et de rayon a; elle est dérivée de la capacité spécifique d'une surface de membrane unité, C [F cm⁻²] en multipliant cette dernière par la surface d'un morceau de longueur unité : $c = C2\pi a$.

* r [Ω cm⁻¹] est la résistance longitudinale d'un morceau de longueur unité et de rayon a; elle est dérivée de la résistance spécifique du milieu intérieur, ρ [Ω cm] en divisant cette dernière par la surface du disque générateur du cylindre / axone : $r = \rho/\pi a^2$.

Si nous utilisons la variable x pour la position le long de l'axone et que nous considérons un morceau de longueur Δx, la loi d'Ohm nous donne au premier ordre en Δx le relation suivante entre potentiel de membrane V(x, t) [mV] et courant longitudinal i(x, t) [mA] :

$$V(x, t) - V(x + \Delta x, t) = i(x + \Delta x/2, t) r \Delta x,$$

souvenez-vous de vos cours de Physique : le courant électrique, comme l'eau d'une rivière, va des valeurs élevées du potentiel vers les valeurs faibles. (Comme le potentiel de membrane est exprimé en mV et que rΔx est exprimé en Ω, i est exprimé en mA.) La conservation de la charge électrique implique (au premier ordre en Δx) :

$$i(x, t) - i(x + \Delta x, t) = \left( c \frac{\partial V(x + \Delta x/2, t)}{\partial t} + gV(x + \Delta x/2, t) \right) \Delta x.$$

Cette équation signifie que si le courant change entre l'entrée et la sortie d'un « petit morceau » d'axone, la différence doit nécessairement correspondre au courant à travers la membrane (qui est compté positivement lorsqu'il est dirigé de l'intérieur vers l'extérieur). Ce dernier à deux composantes : le courant

<footnote>cette dernière correspond à l'unité de mesure de la première.</footnote>
capacitif (premier terme dans la parenthèse du membre de droite) et le courant résistif (second terme). En combinant ces deux dernières équations pour éliminer *i* et en passant à la limite lorsque $\Delta x \to 0$, nous obtenons une version de l'équation de la chaleur :

$$ \frac{1}{rc} \frac{\partial^2 V}{\partial x^2}(x, t) - \frac{\partial V}{\partial t}(x, t) = \frac{g}{c} V(x, t). \qquad (1) $$

Lorsque que de petites déviations de V par-rapport à la valeur de repos ($V = 0$) sont considérées — comme les deux trajectoires noires sur la partie supérieur de la figure 2 —, g peut être considéré comme un scalaire. Pour rendre compte de déviations plus importantes, l'approximation scalaire doit être abandonnée et g doit être vu comme une fonction de la position et du temps.

### 2.2.2 Ce qu'Hodgkin et Huxley ont fait (et ce qu'ils n'ont pas fait)

La contribution expérimentale essentielle d'Hodgkin et Huxley est leur décomposition / description de la conductance membranaire lors du potentiel d'action. De façon un peu contre-intuitive, ils emploient pour cela une approche qui supprime la propagation du potentiel d'action. Au moyen d'un fil conducteur inséré sur toute la longueur de l'axone, ils imposent à celui-ci un potentiel spatialement uniforme (il n'y a donc plus de propagation du potentiel d'action). Comme le laisse entendre cette dernière phrase, ils contrôlent le potentiel de membrane et mesurent le courant membranaire qui en résulte ; les physiologistes disent qu'ils travaillent en « voltage imposé ». En substituant la choline, qui ne passe jamais à travers la membrane, aux ions sodium et en exploitant des différences de cinétiques (fig. 4, partie droite), ils montrent que des conductances spécifiques aux ions sodium ($Na^+$), d'une part, et aux ions potassium ($K^+$), d'autre part, sont modulées dans le domaine de potentiel ($[-10, 100]$ mV) du potentiel d'action. Le membre de droite de l'équation (1) prend alors la forme :

$$ \frac{g}{c} V(x, t) = \frac{1}{c} [g_{Na}(x, t)(V(x, t) - V_{Na}) + g_{K}(x, t)(V(x, t) - V_{K}) + g_{F}(V(x, t) - V_{F})], $$

où $g_{Na}(x, t)$ est la conductance spécifique aux ions sodium et $V_{Na} = 115$ mV est le potentiel d'inversion<fnref n="11" /> de ceux-ci ; $g_{K}(x, t)$ est la conductance spécifique

<footnote n="11">11. C'est le potentiel auquel le courant de diffusion, dû à la différence de concentrations des ions entre l'intérieur et l'extérieur de l'axone, et le courant de conduction, dû à l'effet d'une différence de potentiel (donc d'un champ électrique) sur un ion s'annulent [30, 26].</footnote>
aux ions potassium et $V_K = -12$ mV leur potentiel d'inversion; $g_F$ est une conductance « de fuite », elle ne dépend ni de la position, ni du temps et est utilisée pour avoir un courant membranaire nul lorsque $V = 0$, son potentiel d'inversion est $V_F = +10,6$ mV. Ils proposent la description suivante pour les conductances potassique et sodique :

$$g_K(x, t) = g_{K,max} n^4(x, t)$$

et

$$g_{Na}(x, t) = g_{Na,max} m^3(x, t) h(x, t),$$

où les fonctions « d'activation » $n$, $m$ et $h$ prennent leurs valeurs dans $[0, 1]$ et sont solutions d'équations différentielles du premier ordre dont la forme générale est :

$$\frac{\partial p(x, t)}{\partial t} = - \frac{p(x, t) - p_0(V(x, t))}{\tau_p(V(x, t))}. \qquad (2)$$

Comme leurs indices l'indiquent, $g_{K,max}$ et $g_{Na,max}$ sont les conductances maximales pour les deux ions. La figure 4 montre les graphes des asymptotes ($n_0(V)$, $m_0(V)$, $h_0(V)$) et des constantes de temps ($\tau_n(V)$, $\tau_m(V)$, $\tau_h(V)$) en fonction du potentiel de membrane. Qualitativement, on voit que si le potentiel de membrane augmente, la conductance sodique ($g_{Na} \propto m^3 h$) va aussi augmenter — les physiologistes parlent d'activation —, car $m$ va passer d'une valeur proche de 0 à une valeur proche de 1 avec une cinétique rapide, alors que $h$ va diminuer de 0,9 à 0 avec une cinétique lente. Si le potentiel de membrane reste élevé « longtemps », alors $g_{Na}$ va revenir à 0 — les physiologistes parlent d'inactivation —, car $h$ va atteindre cette valeur. Par contre $g_K \propto n^4$ va croître « lentement » et rester proche de son maximum aussi longtemps que la membrane aura un potentiel élevé.

**Remarque** Hodgkin et Huxley ont effectué des compromis en construisant leur modèle, ils s'en expliquent ainsi [16, p 506] : « Our object here is to find equations which describe conductances with reasonable accuracy and are sufficiently simple for theoretical calculation of the action potential... ». Par « calcul théorique », ils entendent une solution numérique, pas une solution analytique; de plus, nous sommes en 1952, ils ne disposent que de calculatrices mécaniques. En particulier, le choix des puissances ($n^4, m^3$), plutôt que d'équations d'ordre plus élevé est motivé par des considérations pratiques. Il n'est pas surprenant que des modèles de conductances plus précis — basés sur des résultats expérimentaux obtenus avec techniques bien plus sophistiquées — aient été établis depuis [1, 18].
<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Les asymptotes et les constantes de temps des composantes de $g_K$ et $g_{Na}$. Valeurs estimées expérimentalement des 6 fonctions $n_0(V)$, $m_0(V)$, $h_0(V)$ et $\tau_n(V)$, $\tau_m(V)$, $\tau_h(V)$. Les constantes de temps sont très sensibles à la température, celles qui sont représentées correspondent à une température de 18,5 °C.</figcaption></figure>

Le modèle de Hodgkin et Huxley peut maintenant être spécifié, il implique 4 fonctions dépendant de la position le long de l'axone, x, et du temps, t : $V(x, t)$, $n(x, t)$, $m(x, t)$ et $h(x, t)$, qui sont solutions du système :

$$ 
\begin{aligned}
\frac{1}{rc} \frac{\partial^2 V}{\partial x^2}(x, t) - \frac{\partial V}{\partial t}(x, t) &= \frac{1}{c} [g_{Na}(x, t)(V(x, t) - V_{Na}) + \\
& \qquad g_K(x, t)(V(x, t) - V_K) + \\
& \qquad g_F(V(x, t) - V_F)], \\
\frac{\partial n(x, t)}{\partial t} &= - \frac{n(x, t) - n_0(V(x, t))}{\tau_n(V(x, t))}, \\
\frac{\partial m(x, t)}{\partial t} &= - \frac{m(x, t) - m_0(V(x, t))}{\tau_m(V(x, t))}, \\
\frac{\partial h(x, t)}{\partial t} &= - \frac{h(x, t) - h_0(V(x, t))}{\tau_h(V(x, t))}.
\end{aligned}
$$

Comme les outils de calcul (des calculatrices mécaniques) ne leur permet-
taient pas de résoudre numériquement les équations aux dérivées partielles (EDP) ci-dessus, ils ont cherché une solution sous forme d'impulsion se propageant à vitesse constante v. C'est-à-dire qu'ils ont cherché une solution de la forme : $V(x, t) = \tilde{V}(x - vt)$. Si nous introduisons la variable $\xi \equiv x - vt$ [cm], nous pouvons transformer l'équation aux dérivées partielles satisfaite par V en une équation différentielle ordinaire. En effet nous avons :

$$ \frac{\partial V}{\partial x} = \frac{d\tilde{V}}{d\xi} \frac{\partial \xi}{\partial x} = \frac{d\tilde{V}}{d\xi} $$

et

$$ \frac{\partial V}{\partial t} = \frac{d\tilde{V}}{d\xi} \frac{\partial \xi}{\partial t} = -v \frac{d\tilde{V}}{d\xi}. $$

Si nous introduisons simultanément $\tilde{p}(\xi) = p(x, t)$ pour $p = m, n$ et $h, \tilde{g}_{Na} = g_{Na,max} \tilde{m}^3 \tilde{h}$ et $\tilde{g}_K = g_{K,max} \tilde{n}^4$, alors notre système d'EDP devient un système d'équations différentielles ordinaires (EDO) :

$$ \begin{aligned} & \frac{1}{rc} \frac{d^2 \tilde{V}(\xi)}{d\xi^2} + v \frac{d \tilde{V}(\xi)}{d\xi} = \frac{1}{c} [\tilde{g}_{Na}(\xi)(\tilde{V}(\xi) - V_{Na}) + \\ & \qquad \tilde{g}_K(\xi)(\tilde{V}(\xi) - V_K) + \\ & \qquad g_F(\tilde{V}(\xi) - V_F)], \\ & \frac{d \tilde{n}(\xi)}{d\xi} = v \frac{\tilde{n}(\xi) - n_0(\tilde{V}(\xi))}{\tau_n(\tilde{V}(\xi))}, \\ & \frac{d \tilde{m}(\xi)}{d\xi} = v \frac{\tilde{m}(\xi) - m_0(\tilde{V}(\xi))}{\tau_m(\tilde{V}(\xi))}, \\ & \frac{d \tilde{h}(\xi)}{d\xi} = v \frac{\tilde{h}(\xi) - h_0(\tilde{V}(\xi))}{\tau_h(\tilde{V}(\xi))}. \end{aligned} $$

La vitesse v est devenue un paramètre qu'Hodgkin et Huxley ont déterminé (patiemment!) par la méthode du tir (que nous illustrons dans la section suivante avec un cas plus simple). Ils ont ainsi obtenu une vitesse calculée de 18,8 m/s alors que la vitesse expérimentale était 21,2 m/s [16, p 528], pour « axone standard » dont les paramètres sont : $\rho = 35,4$ [ohm cm]; $a = 238$ [$\mu$m]; $r = 20 \times 10^4$ [ohm/cm]; $c = 1,5 \times 10^{-7}$ [farad/cm]. En introduisant $\tilde{W} = d\tilde{V}/d\xi$, on voit que le système ci-dessus devient un système de 5 équations différentielles autonomes du premier ordre. Système qui n'est facile ni à visualiser
ni à analyser. Dans la section suivante nous allons introduire quelques hypothèses simplificatrices qui nous permettront d'étudier analytiquement le front d'onde.

## 2.3 Le front d'onde

L'examen de la partie droite de la figure 4 montre que la constante de temps $τ_m$, déterminant l'évolution de la fonction $\tilde{m}(\xi)$, est typiquement 10 fois plus petite les deux autres constantes de temps, $τ_h$ et $τ_n$. Ceci suggère l'approximation suivante :

$$τ_m \approx 0, \quad τ_h \approx \infty, \quad τ_n \approx \infty.$$

Alors,

$$\tilde{m}(\xi) \approx m_0(\tilde{V}), \quad \tilde{h}(\xi) \approx h_0(0), \quad \tilde{n}(\xi) \approx n_0(0).$$

C'est-à-dire que $\tilde{m}$ atteint instantanément sa valeur d'équilibre, alors que $\tilde{h}$ et $\tilde{n}$ restent à leurs valeurs de départ : $\tilde{h}(0) \approx 0,6$; $\tilde{n}(0)^4 = 0,01$. Sous ces hypothèses, la dynamique du « front d'onde » (fo) est décrite par :

$$\frac{1}{rc} \frac{d^2 \tilde{V}(\xi)}{d\xi^2} + v \frac{d \tilde{V}(\xi)}{d\xi} = \frac{j_{fo}(\tilde{V}(\xi))}{c}, \quad (3)$$

où

$$\begin{aligned} j_{fo}(\tilde{V}) = &g_{Na,max} \tilde{m}^3(\tilde{V}) \tilde{h}(0) (\tilde{V} - V_{Na}) + \\ &g_{K,max} \tilde{n}^4(0) (\tilde{V} - V_K) + \\ &g_F(\tilde{V} - V_F), \end{aligned}$$

dont le graphe est tracé sur la figure 5.

Cette approximation revient, physiologiquement, à bloquer à la fois l'activation des conductances sodiques et l'activation des conductances potasiques<fnref n="12" />. On obtient ainsi un système où le courant sodium,

$$g_{Na,max} \tilde{m}^3(\tilde{V}) \tilde{h}(0) (\tilde{V} - V_{Na}),$$

va très vite dominer les deux autres $g_{K,max} \tilde{n}^4(0) (\tilde{V} - V_K) + g_F(\tilde{V} - V_F)$; ce courant va être négatif tant que $\tilde{V} < V_{Na}$, ce qui va avoir pour effet d'augmenter $\tilde{V}$ au cours du temps, si nous supposons que nous nous plaçons en une

<footnote n="12">12. Il serait aujourd'hui possible de faire cela expérimentalement, même si cela n'a jamais été effectué à ma connaissance.</footnote>
<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Densité de courant membranaire $J_{fo}(\tilde{V})$ en fonction du potentiel pour l'approximation du front d'onde (3). Nous avons $j_{fo} = 2\pi aJ_{fo}$, les paramètres de « l'axone standard » ont été employés : $g_{Na} = 2\pi aG_{Na}$ et $G_{Na} = 120$ [mS/cm$^2$]; $g_K = 2\pi aG_K$ et $G_K = 36$ [mS/cm$^2$]; $g_F = 2\pi aG_F$ et $G_F = 0,3$ [mS/cm$^2$]. Le petit morceau de courbe bleue est un « zoom » sur $J_{fo}(\tilde{V})$ au voisinage de zéro afin de montrer que cette fonction possède 3 zéros dans le domaine couvert par le potentiel d'action. Ces trois zéros sont localisés en $\tilde{V}_0 = 0, \tilde{V}_1 = 2,62$ et $\tilde{V}_2 = 113,92$ [mV].</figcaption></figure>

position déterminée le long de l'axone (et un temps qui croît en un point donné, revient à un $\xi$ qui décroît dans l'équation (3)). Ce processus va persister jusqu'à ce qu'un potentiel $\tilde{V}_2$ soit atteint tel que le courant sodium entrant, $g_{Na,max} \tilde{m}^3(\tilde{V}_2) \tilde{h}(0) (\tilde{V}_2 - V_{Na})$ soit exactement compensé par les deux autres courants sortants, $g_{K,max} \tilde{n}^4(0) (\tilde{V}_2 - V_K) + g_F(\tilde{V}_2 - V_F)$. Une solution de l'équation (3), si elle existe, devrait ainsi correspondre à une vague de potentiel, partant de $\tilde{V} = 0$ pour arriver à $\tilde{V} = \tilde{V}_2$ et se propageant de la gauche vers la droite pour $v > 0$ (ou de la droite vers la gauche si $v < 0$), d'où la dénomination « front d'onde ».

Si nous introduisons $g_{max}$ la pente la plus importante de $J_{fo}(\tilde{V})$ parmi les trois zéros, clairement $g_{max} = \frac{dJ_{fo}(\tilde{V}_2)}{d\tilde{V}}$, nous pouvons normaliser l'équation (3)
de la façon suivante :

$$ \frac{1}{rg_{max}} \frac{d^2 \tilde{V}(\xi)}{d\xi^2} + v \frac{c}{g_{max}} \frac{d \tilde{V}(\xi)}{d\xi} = \frac{j_{fo}(\tilde{V}(\xi))}{g_{max}}. \quad (4) $$

Ici, le membre de droite est le rapport d'un courant par unité de longueur exprimé en mA/cm sur une conductance par unité de longueur exprimée en S/cm, c'est donc une différence de potentiel exprimée en mV. On peut vérifier que $rg_{max}$ correspond à l'inverse d'une distance au carré ($\xi$ est exprimé en cm) et que $\frac{c}{g_{max}}$ est un temps exprimé en seconde ($v$ est une vitesse exprimée en cm/s). Le changement de variable $y = \sqrt{rg_{max}}\xi$ ($y$ est sans dimension), associé aux définitions $U(y) = \tilde{V}(\xi)$, $\tilde{v} = v\sqrt{rc^2/g_{max}}$ ($\tilde{v}$ est sans dimension) et $f(U) = j_{fo}(U)/g_{max}$, nous donnent une version de l'équation (4) dont la variable est sans dimension :

$$ \frac{d^2 U(y)}{dy^2} + \tilde{v} \frac{dU(y)}{dy} = f(U(y)). \quad (5) $$

L'introduction de $W(y) = dU(y)/dy$ nous amène au système d'EDO autonomes du premier ordre :

$$ \frac{dU(y)}{dy} = W(y), \quad (6) $$

$$ \frac{dW(y)}{dy} = f(U(y)) - \tilde{v}W(y). \quad (7) $$

Nous allons voir à présent que même si ce système n'a pas de solution analytique, il est facilement représentable et certaines propriétés de sa solution peuvent être établies.

### 2.3.1 Points singuliers

Le système d'EDO défini par les équations (6) et (7) admet trois points singuliers ($U, W$) : (0,0); ($\tilde{V}_1, 0$) et ($\tilde{V}_2, 0$). Ce système linéarisé aux voisinages des points singuliers ($U_x + \delta U, \delta W$), où $U_x$ est égal à 0, $\tilde{V}_1, \tilde{V}_2$ suivant le point considéré, prend la forme suivante :

$$ \frac{d}{dy} \begin{bmatrix} \delta U \\ \delta W \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ f'(U_x) & -\tilde{v} \end{bmatrix} \begin{bmatrix} \delta U \\ \delta W \end{bmatrix}. $$

La matrice $2 \times 2$ du membre de droite admet deux valeurs propres : $\lambda_{x\pm} = -\tilde{v}/2 \pm \sqrt{f'(U_x) + \tilde{v}^2/4}$. On voit que les deux points extrêmes, (0,0) et ($\tilde{V}_2, 0$),
pour lesquels $f'(U_x) > 0$, sont des points cols (une valeur propre est strictement positive, l'autre est strictement négative), alors que le point intermédiaire, $(\tilde{V}_1, 0)$ est un nœud stable si $\tilde{v} > 0$ et si $\tilde{v}^2/4 > |f'(\tilde{V}_1)|$ (les deux valeurs propres sont strictement négatives) et un foyer stable si la seconde condition est fausse (les deux valeurs propres sont complexes avec une partie réelle négative). Dans le cas où $\tilde{v} > 0$ le nœud est instable de même que le foyer ; si $\tilde{v} = 0$, c'est un centre.

### 2.3.2 Que se passe-t-il si $\tilde{v} = 0$ ?

Nous nous intéressons ici aux trajectoires qui quittent le point (0,0). En posant $\tilde{v} = 0$ dans l'équation (7) puis en « multipliant les deux membres de cette dernière » par ceux de l'équation (6), nous obtenons :

$$W dW = f(U)dU.$$

Ce qui nous donne l'équation suivante satisfaite par $W(U)$ :

$$W^2(U)/2 = \int_0^U f(x)dx. \qquad (8)$$

En nous référant à la figure 5 nous voyons qu'en aucun cas, le $U$ de l'équation (8) ne peut atteindre $\tilde{V}_2$, puisque $\int_0^{\tilde{V}_2} f(x)dx < 0$. Au-delà de $U = \tilde{V}_1$, $W^2(U)/2$ est en effet une fonction strictement décroissante. Il va donc y avoir un $\tilde{V}_1 < U_{max} < \tilde{V}_2$, tel que $W^2(U_{max})/2 = 0$. Ce $U_{max}$ sera la plus grande valeur atteinte par $U$ sur cette trajectoire qui va d'abord « monter » vers $U_{max}$ suivant :

$$W_{\uparrow}(U) = \sqrt{2 \int_{0}^{U} f(x)dx},$$

puis « redescendre » vers 0 suivant :

$$W_{\downarrow}(U) = -\sqrt{2 \int_{0}^{U} f(x)dx}.$$

La partie gauche de la figure 6 montre en orange cette trajectoire obtenue numériquement (par application de la méthode de Simpson sur des intervalles de $U$ de longueur $10^{-4}$). Les lignes orientés bleues représentent le flux associé au champs de vecteurs défini par les équations (6) et (7) (avec $\tilde{v} = 0$). L'épaisseur des courbes est proportionnelle à la vitesse (norme des vecteurs du champ). Le nœuds $(\tilde{V}_1, 0)$ apparaît clairement comme centre des orbites
circulaires. La solution $U(y)$ illustrée sur la partie droite de la figure est obtenue à partir de l'équation (6) qui nous permet de définir $y$ comme fonction de $U$ de la façon suivante :

$$y = \int_{0}^{U} dx/W(x).$$

L'estimation de $W(U)$ obtenue numériquement lors de la construction de la partie gauche de la figure est ici réutilisée et l'intégrale est à nouveau évaluée numériquement.

<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Trajectoire de et vers l'origine à vitesse nulle. À gauche, en orange, les deux parties successives, $W_↑(U)$ et $W_↓(U)$, de la trajectoire solution de (8) avec le champs de vecteurs (en bleu) défini par le système d'EDO (6) et (7); à droite, le profil de la solution en fonction de $y = \sqrt{rg_{max}}\xi$.</figcaption></figure>

### 2.3.3 Propriété « cubique » de $f$ et conséquence : lien entre $\tilde{v}$ et intégrale de $f$

La fonction $f$ est « cubique », au sens où elle a trois zéros dans le domaine qui nous intéresse. Si nous poursuivons les considérations de la section qui nous a amenée au système d'EDO défini par les équations (6) et (7), nous postulons l'existence d'une solution de ce système passant de 0 à $\tilde{V}_2$ lorsque
y passe de $-\infty$ à $+\infty$ dans l'équation (5), nous devons conclure qu'une telle solution serait telle que $U$ tendrait vers une constante pour de grandes valeurs de $|y|$. Nous aurions ainsi $W(\pm\infty) = 0$. Mais les équations (6) et (7) peuvent aussi s'écrire :

$$W dW = f(U)dU - \tilde{v} W^2 dy.$$

Si nous intégrons chacun des deux membres lorsque $y$ va de $]-\infty, +\infty[$, nous obtenons :

$$\int_{0}^{\tilde{V}_2} f(U)dU = \tilde{v} \int_{-\infty}^{\infty} W^2 dy.$$

Comme le membre de gauche est négatif et que l'intégrale du membre de droite est positive, nous devons avoir la vitesse $\tilde{v}$ d'une éventuelle solution strictement négative. Remarquez que le même raisonnement nous donne la contrainte $\tilde{v} > 0$ pour une éventuelle solution allant « dans le sens inverse », c'est-à-dire de $\tilde{V}_2$ à 0.

### 2.3.4 Méthode du tir et portrait de phase

Comme nous ne pouvons pas résoudre notre système analytiquement (du fait de la forme de $f$), nous avons recours à une méthode numérique et nous employons la méthode du tir pour déterminer $\tilde{v}^{13}$. En quelques mots, nous cherchons une trajectoire allant de $(0,0)$ à $(\tilde{V}_2, 0)$. Pour ce faire, nous commençons l'intégration numérique au voisinage immédiat du premier point singulier $(\delta V, \delta V(\sqrt{\tilde{v}^2/4 + f'(0) - \tilde{v}/2}))$, avec $\delta V = 10^{-4}$, puis nous amorçons une mise en œuvre d'une méthode d'Adams-Bashforth d'ordre 5 par 4 pas successifs d'une mise en œuvre d'une méthode de Runge-Kutta d'ordre 4; la taille du pas d'intégration est de $10^{-2}$ (après vérification que des résultats identiques sont obtenus avec un pas 100 fois plus petit). Pour la méthode de tir, nous employons une recherche exploratoire qui nous permet de déterminer deux vitesses telles que $2 \times 10^3$ pas d'intégration donnent des trajectoires dont les valeurs ultimes de $W$ sont négative pour l'une et positive pour l'autre, puis nous itérons une procédure de dichotomie jusqu'à ce que les deux vitesses soient différentes de moins de $10^{-8}$. La partie gauche de la figure 7 montre trois trajectoires ainsi obtenues dont celle correspondant à la meilleure approximation, 1,2112, de $\tilde{v}$ en orange. La partie droite de la figure reproduit cette dernière trajectoire et l'accompagne des lignes de flux définies par le système d'EDO. La conversion de la vitesse non dimensionnée $\tilde{v}$ en vitesse

<footnote n="13">13. Tous les détails de la mise en œuvre des méthodes numériques, ainsi que les programmes correspondant — écrits en Python — sont disponibles sur le dépôt https://gitlab.math.unistra.fr/christophe.pouzat/modelisation-en-neurosciences.</footnote>
<figure><img src="image_8.png" /><figcaption>FIGURE 7 – À gauche, résolution numérique du système d’EDO (6) et (7) par la méthode du tir (en bleu, la vitesse est trop grande, en vert elle est trop petite, en orange elle est « correcte »); à droite, la trajectoire obtenue (orange) et le champs de vecteurs défini par le système d’EDO.</figcaption></figure>

dimensionnée (avec l'équation $\tilde{v} = v\sqrt{rc^2/g_{max}}$) nous donne 59,3 m/s. Cette vitesse est plus grande que les 18,8 m/s obtenus par l'application de la méthode du tir au système de Hodgkin et Huxley, ce qui n'est pas surprenant puisque l'approximation que nous avons considérée ici donne une activation instantanée des courants sodium et exclut le retour du potentiel au niveau de départ (ce retour au niveau de départ a pour effet de ralentir la vitesse de propagation).

Remarquez que $\tilde{v}$ est unique, les trajectoires $W(U; v)$ correspondant à des valeurs de la vitesse $v$ strictement plus grandes que $\tilde{v}$ sont toujours « en-dessous » de $W(U; \tilde{v})$ dans le demi-plan $W > 0$ (par référence à la figure 7). De même, si $v < \tilde{v}$, alors $W(U; v)$ est « au-dessus ». Pour montrer cela, nous utilisons d'abord le fait que la pente des trajectoires qui quittent le point (0,0) pour aller vers $(\tilde{V}_2, 0)$ est $\sqrt{v^2/4 + f'(0)} - v/2$<fnref n="14" />, c'est-à-dire une fonction dé-

<footnote n="14">14. C'est le rapport des coordonnées du vecteur propre correspondant à la valeur propre positive en ce point col.</footnote>
croissante de $v$. Ainsi si $v > \tilde{v}$, alors $W(U; v)$ est au-dessous de $W(U; \tilde{v})$ au voisinage immédiat de $(0,0)$. Pour que $W(U; v)$ rejoigne $W(U; \tilde{v})$, ou « passe en-dessus » de celle-ci, en un $U$ donné, il faudrait qu’en ce point nous ayons $dW(U; v)/dU \ge dW(U; \tilde{v})/dU$ (inégalité stricte en cas de croisement). Mais la combinaison des équations (6) et (7) nous donne :

$$ \frac{dW}{dU} = \frac{f(U)}{W} - v $$

qui est une fonction strictement décroissante de $v$ en tout point de l'espace des phases. Au point de réunion ou de croisement nous devrions ainsi avoir si-multanément $dW(U; v)/dU \ge dW(U; \tilde{v})/dU$ et $dW(U; v)/dU < dW(U; \tilde{v})/dU$, ce qui est impossible.

### 2.3.5 Approximation cubique de $f$ et solution analytique (Zeldovich et Frank-Kamenetsky, 1938)

Une approximation encore plus « radicale » consiste à remplacer la fonction $f$ de l'équation (5) par un polynôme du troisième ordre :

$$ f_{ZFK}(\tilde{V}) = \frac{\tilde{V}(\tilde{V} - \tilde{V}_1)(\tilde{V} - \tilde{V}_2)}{\tilde{V}_2(\tilde{V}_2 - \tilde{V}_1)}, $$

où l'indice ZFK rappelle que cette approximation a d'abord été proposée en 1938 par Zeldovich et Frank-Kamenetsky qui étudiaient la propagation d'un front de flamme [26, p 102]. Avec cette approximation, l'unique $\tilde{v}$ est facilement obtenue en utilisant la plus simple forme fonctionnelle $W(U)$ ayant précisément deux zéros, un en 0 et l'autre en $\tilde{V}_2$, c'est-à-dire $W = KU(U - \tilde{V}_2)$ (ou $K$ est une constante strictement négative à déterminer), pour aboutir à $K = -1/\sqrt{2\tilde{V}_2(\tilde{V}_2 - \tilde{V}_1)}$ :

$$ \tilde{v}_{ZFK} = - \frac{\tilde{V}_2 - 2\tilde{V}_1}{\sqrt{2\tilde{V}_2(\tilde{V}_2 - \tilde{V}_1)}}. $$

La vitesse dimensionnée est alors 33 m/s. On remarque que la vitesse dépend essentiellement du rapport $\tilde{V}_1/\tilde{V}_2$. La solution du système d'EDO est alors :

$$ U_{ZFK}(y) = \frac{\tilde{V}_2}{1 + \exp(K\tilde{V}_2y + \text{cste})}. $$

La figure 8 permet de comparer la solution (numérique) $U(y)$ de la première approximation à la solution (analytique) $U_{ZFK}(y)$ de la seconde.
<figure><img src="image_9.png" /><figcaption>FIGURE 8 – Comparaison des profils de fronts d'onde : $U(y)$ et $U_{ZFK}(y)$.</figcaption></figure>

## 2.4 Pour aller plus loin

### 2.4.1 Stabilité du front d'onde

La stabilité de la solution de l'équation (3) se propageant à vitesse constante a été prouvée par Zeldovich et Barenblatt [31]. Il ne s'agit pas d'une stabilité asymptotique — où une perturbation devient « aussi petite que l'on veut » à mesure que le temps passe —, car une petite perturbation se traduit par un décalage dans le temps de la solution de départ; ce que nous pouvons noter de façon informelle : $\tilde{V}(\xi) + \text{perturbation} \to \tilde{V}(\xi + \xi_0)$. Dans ce cas, une perturbation initiale « persiste » (sous forme d'un décalage), elle ne tend pas vers 0 (cas de la stabilité asymptotique), pas plus qu'elle n'explose (instabilité).

### 2.4.2 Un premier modèle avec repolarisation

Du point de vue physiologique, les fronts d'ondes décrivent un phénomène qui n'a lieu qu'une seule fois — il n'y a pas de retour à l'état de départ —, mais pour que le potentiel d'action remplisse sa fonction (physiologique), les neurones doivent pouvoir en générer un grand nombre. Nous sommes donc très intéressés par des modèles « avec repolarisation ». Nous avons aussi mentionné dans les sections précédentes qu'une des raisons pour lesquelles les vitesse des fronts d'ondes sont plus grandes que la vitesse du potentiel d'action du modèle de Hodgkin et Huxley est l'absence de repolarisation. Le modèle
introduit par Markine et Tchizmadjev en 1967 [23, chap. 4] va nous permettre d'étayer notre propos. Nous reprenons le cheminement qui nous a amené du modèle de Hodgkin et Huxley (Sec. 2.2.2) à l'équation du front-d'onde (3) :

$$ \frac{1}{rc} \frac{\partial^2 V(x, t)}{\partial x^2} - \frac{\partial V(x, t)}{\partial t} = \frac{j_{mc}(x, t)}{c}, $$

mais $j_{mc}(x, t)$ n'est plus une fonction « cubique » du potentiel de membrane. C'est une fonction nulle tant que le potentiel de membrane est inférieur à un potentiel seuil, $V_\theta > 0$. Lorsque $V(x, t) \ge V_\theta$ en un temps $t$, que nous pren-drons comme origine, $j_{mc}(x, t)$ est définie comme une fonction constante par morceaux :

$$ j_{mc}(x, t) = \begin{cases} 0 & , \text{si } t < 0, \\ -j_1 & , \text{si } 0 \le t < \tau_1, \\ j_2 & , \text{si } \tau_1 \le t < \tau_1 + \tau_2, \\ 0 & , \text{si } \tau_1 + \tau_2 \le t. \end{cases} $$

Ici, $j_1, j_2, \tau_1, \tau_2 > 0$ et $j_1\tau_1 = j_2\tau_2$ (la quantité de charge qui traverse la mem-brane pendant un potentiel d'action est alors nulle). Ce modèle propose une formalisation élémentaire de la séquence d'événements décrite par les phy-siologistes : comme l'illustre la Fig. 2, il y a, au moins de façon informelle, un seuil qui sépare deux comportements qualitativement différents du poten-tiel de membrane; la partie rapidement croissante du potentiel d'action est due à un courant entrant (le courant d'ions sodium domine alors largement le courant d'ions potassium); la partie décroissante est due à un courant sortant [16, courbe $I_i$, Fig. 18C, p. 531]$^{15}$. Comme précédemment, nous cherchons une solution sous la forme d'une impulsion se déplaçant à vitesse constante $v$, $V(x, t) = \tilde{V}(x - vt) = \tilde{V}(\xi)$. Nous obtenons alors :

$$ \frac{1}{rc} \frac{d^2 \tilde{V}(\xi)}{d\xi^2} + v \frac{d \tilde{V}(\xi)}{d\xi} = \frac{\tilde{j}_{mc}(\xi)}{c}, $$

où $\tilde{j}_{mc}(\xi) = 0$ si $\tilde{V} < V_\theta$ et, si $\tilde{V} \ge V_\theta$ :

$$ \tilde{j}_{mc}(\xi) = \begin{cases} 0 & , \text{si } \xi > 0, \\ -j_1 & , \text{si } -v\tau_1 < \xi \le 0, \\ j_2 & , \text{si } -v(\tau_1 + \tau_2) < t \le -v\tau_1, \\ 0 & , \text{si } t \le -v(\tau_1 + \tau_2). \end{cases} $$

<footnote n="15">15. La convention courant entrant négatif / courant sortant positif est postérieure à l'article [16], où la convention opposée est employée.</footnote>
Nous cherchons une ou des solutions qui partent de et qui reviennent à 0 ($\tilde{V}(\xi) = 0$ pour $|\xi|$ suffisamment grand). Nous considérons quatre intervalles : $(I) \equiv ]0, \infty[$, $(II) \equiv ]-v\tau_1, 0[$, $(III) \equiv ]-v(\tau_1 + \tau_2), -v\tau_1[$ et $(IV) \equiv ]-\infty, -v(\tau_1 + \tau_2)]$. Sur chacun de ces intervalles, nous devons résoudre une équation du seconde ordre, dont la forme générale sera : $A_i + B_i\xi + C_i \exp(-vrc\xi)$ et notre solution globale doit être dérivable à dérivée continue sur $\mathbb{R}$, ce qui nous génère des contraintes sur les $A_i, B_i, C_i$. La condition $\tilde{V}(\xi) = 0$ pour $|\xi|$ impose : $A(I) = B(I) = 0$ et $A(IV) = B(IV) = C(IV) = 0$. À partir de là, avec un peu de patience, nous obtenons [23, p. 111-112]$^{16}$ :

$$
\begin{align*}
V_{(I)}(\xi) &= \frac{e^{-vrc\xi}}{v^2rc^2} (j_1 + j_2 e^{-v^2rc(\tau_1+\tau_2)} - (j_1 + j_2)e^{-v^2rc\tau_1}), \\
V_{(II)}(\xi) &= \frac{e^{-vrc\xi}}{v^2rc^2} (j_2 e^{-v^2rc(\tau_1+\tau_2)} - (j_1 + j_2)e^{-v^2rc\tau_1}) \\
&\qquad - \frac{j_1\xi}{vc} + \frac{j_1}{v^2rc^2}, \\
V_{(III)}(\xi) &= \frac{j_2 e^{-v^2rc(\tau_1+\tau_2)} e^{-vrc\xi}}{v^2rc^2} + \frac{j_2\xi}{vc} \\
&\qquad + \frac{j_2(\tau_1+\tau_2)}{c} - \frac{j_2}{v^2rc^2}, \\
V_{(IV)}(\xi) &= 0.
\end{align*}
$$

Nous déduisons de $V_{(I)}(0) = V_{(II)}(0) = V_{\theta}$ que :

$$
V_{\theta} = \frac{1}{v^2 r c^2} (j_1 + j_2 e^{-v^2 r c (\tau_1 + \tau_2)} - (j_1 + j_2) e^{-v^2 r c \tau_1}). \quad (9)
$$

v étant un paramètre, nous pouvons obtenir sa, ou ses, valeurs possibles en
construisant le graphe du membre de droite et en cherchant sa ou ses inter-
sections avec une droite horizontale dont l'ordonnée est $V_{\theta}$. En employant
les paramètres de « l'axone standard » de calamar : $r = 2 \times 10^6 \ \Omega/m$; $c =$
$1,5 \times 10^{-5} \ \mathrm{F/m}$; $\tau_1 = 2 \times 10^{-4} \ \mathrm{s}$; $\tau_1 = 2 \times 10^{-3} \ \mathrm{s}$; $j_1 = 3,18 \times 10^{-3} \ \mathrm{A/m}$;
$j_2 = 3,18 \times 10^{-4} \ \mathrm{A/m}$; nous obtenons le graphe de la partie gauche de la fi-
gure 9. Ici la ligne à tirets grise figure $V_{\theta}$, la ligne pointillée grise correspond
au maximum (28,3 mV) et les deux lignes verticales matérialisent les deux
vituesses possibles : 4,89 m/s (bleu) et 16,7 m/s (orange). Il y a donc deux so-
lutions dont les graphes sont présentés sur la partie droite de la figure (en
orange, la solution rapide et en bleu, la solution lente). Nous pouvons voir de
façon informelle que la solution rapide est stable alors que la lente ne l'est pas

<footnote n="16">16. Il y a une erreur dans [23, eq. (4.5), p. 111].</footnote>
[23, p. 112-113]. Considérons en effet la solution rapide (orange sur la figure) et imaginons qu'une petite perturbation augmente légèrement la vitesse, alors le potentiel à l'avant de l'impulsion décroît légèrement (fig. 9, partie gauche, à droite la ligne verticale orange) et est inférieur au seuil de déclenchement du courant entrant ($-j_1$), ce qui va ralentir l'impulsion (il faudra « attendre un peu » pour atteindre le seuil). Si la perturbation initiale diminue un peu la vitesse, alors le seuil sera atteint un peu plus tôt, ce qui aura pour conséquence d'augmenter la vitesse. Si nous considérons maintenant la solution lente, nous voyons que l'inverse se produit, une perturbation qui augmente la vitesse, fera que le seuil est atteint plus vite ce qui va encore augmenter la vitesse, alors qu'une perturbation opposée induira une diminution supplémentaire de la vitesse. L'équation (9) nous permet aussi de voir que la présence de la repolarisation ($j_2 > 0$) ralentit l'impulsion. Si nous réarrangeons le terme entre parenthèse avec $j_1$ d'un côté et $j_2$ de l'autre :

$$V_{\theta} = \frac{1}{v^2 r c^2} (j_1(1 - e^{-v^2 r c \tau_1}) + j_2(e^{-v^2 r c \tau_2} - 1)e^{-v^2 r c \tau_1}),$$

nous voyons que poser $j_2 = 0$ augmente le numérateur; augmentation qu'il faudra compenser par une augmentation de $v$ pour garder le même seuil. Enfin, les vitesses de propagations que nous avons obtenues (4,89 et 16,7 m/s) sont remarquablement bonnes puisqu'une des études numériques du modèle de Hodgkin et Huxley ont produit une impulsion instable et lente, dont la vitesse est 5,66 m/s [19] et une impulsion stable et rapide de vitesse 18,8 m/s [16, p. 528].

### 2.4.3 Le modèle de FitzHugh-Nagumo : le modèle avec repolarisation le plus étudié

Le modèle de Markine-Tchizmadjev paie sa simplicité analytique par un côté ad hoc : nous commençons avec une solution numérique du modèle de Hodgkin et Huxley pour obtenir la courant membranaire associé à un potentiel d'action en fonction du temps; nous construisons une approximation constante par morceaux de ce courant membranaire ($j_{mc}(x, t)$); nous « injectons » cette approximation dans notre équation de diffusion du potentiel membranaire. Partant de la différence d'ordre de grandeur des constantes de temps du modèle de Hodgkin et Huxley (partie droite de la Fig. 4) : l'activation des conductances sodiques ($\tau_m$ inférieur à 0,5 millisecunde) d'une part; l'inactivation de ces conductances ($\tau_h$) et l'activation des conductances potassiques ($\tau_n$) supérieure à 5 millisec secondes d'autre part; FitzHugh [14] et Nagumo et
<figure><img src="image_10.png" /><figcaption>FIGURE 9 – Vitesses de propagation et solutions du modèle de Markine et Tchizmadjev. À gauche, le graphe du membre de droite de (9), avec le seuil (tirets), le potentiel maximal (pointillés) et les deux vitesses possibles (lignes verticales bleu et orange); à droite, les deux solutions, la rapide et stable en orange et la lente et instable en bleu.</figcaption></figure>

ses collaborateurs [24] ont proposé le modèle :

$$ \frac{\partial^2 V(x, t)}{\partial x^2} - \frac{\partial V(x, t)}{\partial t} = f(V(x, t)) + R(x, t), \\ \frac{\partial R(x, t)}{\partial t} = \epsilon (V(x, t) + c - bR(x, t)). $$

Ici, nous employons la version non dimensionnée des variables x et t que nous avons introduites juste avant (5). $\epsilon$, $c$ et $b$ sont des paramètres dont l'ajustement permet d'approcher les résultats du modèle de Hodgkin et Huxley et $f$ est (typiquement) une fonction cubique ou constante par morceaux. L'idée est que $f$ rend compte de l'activation des conductances sodiques, alors que $R$ — la lettre a été choisie car elle décrit la « récupération » (recovery) du potentiel de membrane; c'est-à-dire la repolarisation — décrit la combinaison de l'inactivation des conductances sodiques et de l'activation des conductances potassiques. Nous remarquons que le choix $\epsilon = 0$ nous ramène au cas du front d'onde. Des études numériques résumées dans [26, Sec. 6.3, p. 124-127] ont montré que pour $0 < \epsilon \ll 1$ avec une fonction $f$ cubique, le système
ci-dessus admet deux solutions sous forme d'une impulsion se déplaçant à vi-
tesse constante : une rapide et stable ; une lente et instable (comme le modèle
de Markine-Tchizmadjev). Avec une fonction f constante par morceaux, des
solutions sous formes d'ondes (périodiques) se déplaçant à vitesse constante
sont observées, avec encore une fois une solution rapide et stable et une solu-
tion lente et instable. D'importants articles de John Evans au début des années
70 [9, 10, 11, 12] ont établi des critères de stabilité pour les solutions sous
forme d'impulsion se déplaçant à vitesse constante, aussi bien pour le modèle
de Hodgkin et Huxley que pour le modèle de FitzHugh-Nagumo; Jones [21]
a montré que ces critères étaient remplis pour ce dernier modèle. La littéra-
ture sur le modèle de FitzHugh-Nagumo est aujourd'hui très abondante et
l'ouvrage d'Alwyn Scott [26], souvent cité dans le présent article, comme la
conférence d'Olivier Faugeras [13], en sont d'excellents points d'entrée.

2.4.4 Aspects stochastiques

Comme je l'ai expliqué, Hodgkin et Huxley ignoraient la nature physique
des conductances qu'ils ont décrites, mais nous savons désormais que ces
conductances sont des *canaux ioniques* : des protéines insérées dans la mem-
brane cellulaire qui forment des pores sélectifs à certains ions et dont la confor-
mation (état ouvert ou fermé du pore) dépend du potentiel de membrane.
Nous observons de plus expérimentalement des transitions rapides<sup>17</sup> et ré-
versibles entre les états ouverts et fermés d'un canal donné. Les conductances
déterministes du modèle de Hodgkin et Huxley sont donc, en réalité, des
sommes de toutes petites conductances qui fluctuent fortement [29, 8, 18].
Savoir si le modèle de Hodgkin et Huxley peut être vu comme la limite dé-
terministe d'une somme de processus devient alors pertinente. Cette question
génère aujourd'hui beaucoup d'activité suite à la publication des articles [2,
25].

3 Conclusions

J'espère que les lecteurs suffisamment patients pour en arriver à ces lignes
ne m'en voudront pas trop d'avoir consacré autant d'espace à la modélisa-
tion de ce qui n'est qu'une partie, même si elle est essentielle, d'un neurone!
L'exemple du potentiel d'action est, je crois, *représentatif* du type de mathé-
matiques qui intervient en neurosciences : une succession de modèles allant

<footnote n="17">17. Les temps caractéristiques sont de l'ordre de la dizaine de microseconde à la milliseconde.</footnote>
d'images « grossières » de la réalité à des constructions plus affinées; ces modèles posent un cadre dans lequel nous pouvons faire des mathématiques et où approches analytiques et numériques se complètent et se renforcent. Les ouvrages cités dans mon résumé devraient permettre aux lecteurs souhaitant approfondir sur tel ou tel aspect de la foisonnante recherche en « neuromathématiques » d'aller plus loin.

J'ai commencé par une question posée par les neurophysiologistes : « comment obtient-on la propagation fiable et rapide d'un signal d'un bout à l'autre d'un neurone? » avant de plonger dans les aspects mathématiques et, un peu, numériques de modèles proposés pour y répondre. Mes lecteurs pourraient légitimement se demander, si nous avons effectivement produit des résultats qui intéressent les physiologistes. La réponse est affirmative au moins de quatre points de vue différents. Tout d'abord, les quelques modèles que nous venons de discuter admettent tous une unique solution stable sous forme de front d'onde (Zeldovich et Frank-Kamenetsky) ou d'impulsion (Markine-Tchizmadjev et Fitzhugh-Nagumo). C'est une propriété essentielle pour qui veut comprendre l'organisation de bon nombre de circuits neuronaux où les délais de propagation jouent un rôle prépondérant (c'est le cas de tout système d'anticipation motrice, comme de la localisation d'une source sonore à haute fréquence chez les vertébrés à « grosses têtes » dont les chouettes et nous). Ensuite, les neurophysiologistes sont très friands de modèles d'EDP, réductionnistes, détaillés qui ne peuvent être abordés que numériquement. Mais la meilleure façon de se convaincre que les algorithmes et programmes qu'ils emploient sont raisonnablement fiables est de suivre la voie esquissée ici : obtenir une solution explicite dans un cas simple ; utiliser cette solution pour valider notre mise en œuvre de la méthode du tir ; utiliser la vitesse de propagation obtenue par cette dernière méthode pour vérifier la solution générée par le programme de résolution du système d'EDP. Dans le même ordre d'idées, l'impulsion qui se propage à vitesse constante le long de l'axone dispensent les numériciens de modéliser ce dernier explicitement lorsqu'ils simulent de grands réseaux, ce qui représente un gain de temps (de calcul) considérable. Enfin, l'approche de Hodgkin et Huxley : séparation du courant de membrane en composantes distinctes ; description quantitative de ces dernières ; a constitué un programme de recherche extrêmement fructueux en neurophysiologie. De nouvelles conductances ont été ainsi identifiées, quantifiés et incorporées dans des généralisations du modèle de H & H qui s'y prête particulièrement bien. Les expérimentateurs produisent des quantités toujours plus importantes de données et les modèles quantitatifs et « assimilateurs » sont une vraie rareté en Biologie !
**Je remercie**: Philomène Picart de la bibliothèque inter-universitaire de santé, Université Paris-Cité, pour m'avoir fourni une copie de haute résolution de la figure 24 de S. Ramón y Cajal [4]; Dasha Loukianova pour ses commentaires et son aide sur l'orthographe des noms russes; mes deux « rapportrices », Michèle Romanos et Eva Löcherbach pour leurs critiques et remarques aussi constructives que bienveillantes; mon éditrice, Magali Ribot, pour son immeuble patience à mon égard. Ce travail est soutenu par l'ANR : projet ANR-22-CE45-0027.

Références

[1] C M ARMSTRONG. « Sodium channels and gating currents. » In : *Physiological Reviews* 61.3 (1981), p. 644-683. URL : https://doi.org/10.1152/physrev.1981.61.3.644.

[2] Tim D. AUSTIN. « The emergence of the deterministic Hodgkin–Huxley equations as a limit from the underlying stochastic ion-channel mechanism ». In : *The Annals of Applied Probability* 18.4 (août 2008). ISSN : 1050-5164. URL : http://dx.doi.org/10.1214/07-AAP494.

[3] Mostafa BACHAR, Jerry BATZEL et Susanne DITLEVSEN, éd. *Stochastic Bio- mathematical Models : with Applications to Neuronal Modeling*. Springer Berlin Heidelberg, 2012.

[4] S. RAMÓN CAJAL. Histologie du système nerveux de l'homme et des vertébrés. Édition française, revue et mise à jour par l'auteur, traduite de l'espagnol par le Dr. L. Azoulay. Tome 1. A. Maloine, Paris, 1909. URL : https://www.biusante.parisdescartes.fr/histmed/medica/cote?ryc004x01.

[5] Neil CAMPBELL et al. *Biology : A Global Approach*. 12e éd. Pearson, 2020.

[6] Ann M. CASTELFRANCO et Daniel K. HARTLINE. « Evolution of rapid nerve conduction ». In : *Brain Research* 1641 (2016), p. 11-33.

[7] J.W. COOLEY et F.A. DODGE. « Digital Computer Solutions for Excitation and Propagation of the Nerve Impulse ». In : *Biophys J* 6.5 (sept. 1966), p. 583-599.

[8] Louis J. DEFELICE. *Introduction to Membrane Noise*. Plenum Press, 1981.

[9] JOHN W. EVANS. « Nerve Axon Equations : 1 Linear Approximations ». In : *Indiana University Mathematics Journal* 21.9 (1972), p. 877-885. ISSN : 00222518, 19435258. URL : http://www.jstor.org/stable/24890354 (visité le 25/11/2024).
[10] JOHN W. EVANS. « Nerve Axon Equations : II Stability at Rest ». In : *Indiana University Mathematics Journal* 22.1 (1972), p. 75-90. ISSN : 00222518, 19435258. URL : http://www.jstor.org/stable/24890397 (visité le 25/11/2024).

[11] JOHN W. EVANS. « Nerve Axon Equations : III Stability of the Nerve Impulse ». In : *Indiana University Mathematics Journal* 22.6 (1972), p. 577-593. ISSN : 00222518, 19435258. URL : http://www.jstor.org/stable/24890484 (visité le 25/11/2024).

[12] JOHN W. EVANS. « Nerve Axon Equations : IV The Stable and the Unstable Impulse ». In : *Indiana University Mathematics Journal* 24.12 (1975), p. 1169-1190. ISSN : 00222518, 19435258. URL : http://www.jstor.org/stable/24890975 (visité le 25/11/2024).

[13] Olivier FAUGERAS. *Mathematical neuroscience*. IHP. Audiovisual resource. 6 mars 2023. URL : https://dx.doi.org/10.57987/IHP.2023.T1B.WS1.001.

[14] Richard FITZHUGH. « Impulses and Physiological States in Theoretical Models of Nerve Membrane ». In : *Biophysical Journal* 1.6 (1961), p. 445-466. ISSN : 0006-3495. URL : https://www.sciencedirect.com/science/article/pii/S0006349561869026.

[15] Antonio GALVES, Eva LÖCHERBACH et Christophe POUZAT. *Probabilistic Spiking Neuronal Nets*. Springer Cham, 2024.

[16] A.L. HODGKIN et A.F. HUXLEY. « A quantitative description of membrane current and its application to conduction and excitation in nerve ». In : *The Journal of Physiology* 117.4 (1952), p. 500-544. URL : https://physoc.onlinelibrary.wiley.com/doi/abs/10.1113/jphysiol.1952.sp004764.

[17] Alan Lloyd HODGKIN et William Albert Hugh RUSHTON. « The electrical constants of a crustacean nerve fibre ». In : *Proceedings of the Royal Society of London. Series B - Biological Sciences* 133.873 (1946), p. 444-479. URL : https://royalsocietypublishing.org/doi/abs/10.1098/rspb.1946.0024.

[18] R HORN et C A VANDENBERG. « Statistical properties of single sodium channels. » In : *Journal of General Physiology* 84.4 (oct. 1984), p. 505-534. ISSN : 0022-1295. URL : https://doi.org/10.1085/jgp.84.4.505.

[19] A. F. HUXLEY. « Can a nerve propagate a subthreshold disturbance? » In : *The Journal of Physiology* 148.suppl (1959), p. 80-81.
[20] Eugene M. IZHIKEVICH. Dynamical Systems in Neuroscience : The Geometry of Excitability and Bursting. The MIT Press, juill. 2006. ISBN : 9780262276078.
URL : http://dx.doi.org/10.7551/mitpress/2526.001.0001.

[21] Christopher K. R. T. JONES. « Stability of the travelling wave solution of the FitzHugh-Nagumo system ». In : Trans. Amer. Math. Soc. 286.2 (1984), p. 431-469. ISSN : 0002-9947. URL : https://doi.org/10.2307/1999806.

[22] James KEENER et James SNEYD. Mathematical Physiology. I : Cellular Physiology. Interdisciplinary Applied Mathematics. Springer New York, 2009.
ISBN : 9780387758473. URL : http://dx.doi.org/10.1007/978-0-387-75847-3.

[23] V. S. MARKIN, V. F. PASTUSHENKO et Y. A. CHIZMADZHEV. *Theory of Excitable Media*. Wiley, 1987.

[24] J. NAGUMO, S. ARIMOTO et S. YOSHIZAWA. « An Active Pulse Transmission Line Simulating Nerve Axon ». In : Proceedings of the IRE 50.10 (1962), p. 2061-2070.

[25] K. PAKDAMAN, M. THIEULLEN et G. WAINRIB. « Fluid limit theorems for stochastic hybrid systems with application to neuron models ». In : *Advances in Applied Probability* 42.3 (2010), p. 761-794.

[26] Alwyn SCOTT. Neuroscience. A Mathematical Primer. Springer New York, NY, 2002. ISBN : 978-0-387-95402-8. URL : https://doi.org/10.1007/b98897.

[27] Alwyn C. SCOTT. « The electrophysics of a nerve fiber ». In : Rev. Mod. Phys. 47 (2 avr. 1975), p. 487-533. URL : https://link.aps.org/doi/10.1103/RevModPhys.47.487.

[28] Olaf SPORNS. Networks of the Brain. The MIT Press, oct. 2010. ISBN : 9780262289788. URL : https://doi.org/10.7551/mitpress/8476.001.0001.

[29] A.A. VERVEEN et H.E. DERKSEN. « Fluctuation phenomena in nerve membrane ». In : Proceedings of the IEEE 56.6 (1968), p. 906-916. ISSN : 0018-9219. URL : http://dx.doi.org/10.1109/proc.1968.6443.

[30] M. VOLKENSTEIN. Biophysique. Éditions Mir, 1985. URL : https://archive.org/details/m-volkenstein-biophysique-mir-1985.

[31] Y.B. ZELDOVICH et G.I. BARENBLATT. « Theorv of flame propagation ». In : Combustion and Flame 3 (1959), p. 61-74. ISSN : 0010-2180. URL : https://www.sciencedirect.com/science/article/pii/0010218059900070.