Polisano Kévin Etude mathématique et simulation numérique de la dynamique neuronale

De Ensiwiki
Aller à : navigation, rechercher


Etude mathématique et simulation numérique de la dynamique neuronale

Labo INRIA
Equipe BIPOP
Encadrants arnaud.tonnelier@inria.fr


Etudiant

Polisano Kévin, filière MMIS - IRVM


Introduction

Neurone.jpg

La modélisation est devenue un outil essentiel en neurosciences notamment pour explorer la nature du "code neuronal", question qui reste très largement ouverte encore aujourd'hui. Différents projets internationaux ont pour objectifs la modélisation et la simulation de réseaux de neurones de grande taille et cherchent à reproduire, à analyser et à comprendre le fonctionnement de structures neuronales. Il convient alors d'utiliser des modèles de neurones réalistes (biologiquement plausibles) et dont la simulation soit peu coûteuse.

Les modèles neuronaux de type intègre-et-tire (integrate-and-fire neural models) constituent une classe de modèles très largement utilisés en neuroscience computationnelle. Ces modèles qui s'écrivent comme des systèmes differentiels non-réguliers possèdent une dynamique très riche (points fixes, excitabilité, oscillations, chaos). Un élément essentiel des modèles de type intègre-et-tire est l'équation contrôlant le seuil d'un potentiel d'action.

L'objectif de cette IRL était de proposer une généralisation de cette équation basée sur la réduction dimensionnelle de modèles détaillés.


Eléments de prérequis

Avant de travailler sur cet axe de recherche il est nécessaire de comprendre le fonctionnement des modèles d'origine utilisés pour décrire le neurone. Les éléments prérequis sont donc les suivants :

  • Présentation du modèle de Hodgkin et Huxley
  • Réduction dimensionnelle du modèle et étude mathématique
  • Les modèles intègre-et-tire

Le modèle de Hodgkin et Huxley

Au début des années 1950 deux scientifiques, Alan Hodgkin et Andrew Huxley ont effectué des expériences sur un axone géant de calamar. Cet axone est remarquable de par son diamètre extraordinairement large (de l'ordre du millimètre), ce qui leur permis de réaliser des manipulations techniquement infaisables sur des neurones traditionnels trop petits. La contraction d'une tentacule est due à un courant électrique généré par les neurones du calamar, lorsqu'il repère un prédateur par exemple. Hodgkin et Huxley ont interprété ces potentiels d'actions en terme de migration d'ions sodium Na^+ et potassium K^+ à travers la membrane du neurone. Basé sur des mesures d'électrodes ils ont proposé un modèle mathématique connu sous le nom de modèle de Hodgkin-Huxley (HH).


Echanges ioniques à travers la membrane
Allure d'un potentiel d'action
Le neurone est une cellule nerveuse comportant donc une membrane, plus ou moins perméable aux ions. Celle-ci contient en fait des canaux régulant le flux d'ions passant à travers la membrane. La variation de concentration de ces ions de par et d'autre de la membrane induit une différence de potentiel appelé potentiel membranaire. Typiquement au repos le potentiel membranaire d'un neurone est de -70 mV. Lorsqu'un flux d'ions positifs pénètre dans le cellule, le potentiel augmente et peut donner lieu à un potentiel d'action, on dit que le neurone décharge (il transmet l'information sous forme de courant électrique).
Circuit électrique modélisant un neurone
Hodgkin et Huxley ont modélisé le neurone par un circuit électrique, dans lequel on peut établir une équation différentielle régissant le potentiel membranaire V.

Par ailleurs chaque canal est constitué de "portes" qui doivent toutes être ouvertes pour laisser passer les ions (et donc le courant). On observe 3 types de portes et on dénote par m, n et h respectivement la fraction de portes ouvertes de ce type. Ces probabilités varient également au cours du temps suivant une équation différentielle. Le comportement neuronal est ainsi décrit par un système dynamique à 4 variables :


\left\{\begin{array}{l}
c_M\frac{dV}{dt}=-\bar{g}_{Na}m^3h(V-E_{Na})-\bar{g}_Kn^4(V-E_K)-\bar{g}_L(V-E_L)
\\\\\frac{dn}{dt}=\phi \frac{n-n_{\infty}(V)}{\tau_{n}(V)}
\\\\\frac{dm}{dt}=\phi \frac{m-m_{\infty}(V)}{\tau_{m}(V)}
\\\\\frac{dh}{dt}=\phi \frac{h-h_{\infty}(V)}{\tau_{h}(V)}
\end{array}\right.

Réduction dimensionnelle et étude mathématique

En remarquant que n+h\simeq \text{constant} et que m(t)\simeq m_{\infty}(V(t)) on élimine 2 variables dans le système, qui est maintenant de la forme :


 \left\{\begin{array}{l}
\frac{dV}{dt}=f(V,n)
\\\\\frac{dn}{dt}=g(V,n)
\end{array}\right.


L'intérêt de ce modèle réduit est que l'on peut représenter les trajectoires dans le plan de phase et procéder à une étude mathématique plus approfondie. Ci-dessous on simule le modèle pour 2 conditions initiales (de potentiel) différentes (courbes rouge et bleue) et 2 échelons de courants différents :



Remarquez qu'on a fait également figurer dans le plan de phase des courbes verte et rose qui correspondent respectivement aux courbes f(V,n)=0 et g(V,n)=0 appelés V-nullcline et n-nullcline. Le point fixe du système se situe à l'intersection de celles-ci (là où les 2 dérivées s'annulent).

On constate qu'en fonction de la valeur du courant appliqué la réponse du neurone n'a pas le même comportement : à 60 mA la trajectoire converge vers le point fixe (stable, retour à l'état de repos) tandis qu'à 100 mA elle décrit une courbe fermée appelée cycle limite (instable, réponse périodique du neurone). L'étude d'une jacobienne nous renseigne sur la stabilité du point fixe, caractérisée par sa position sur la V-nullcline (dépendant du courant appliqué). Cette sensibilité du système dynamique au paramètre I_{app} se nomme bifurcation de Hopf. La stabilité change en fonction de I_{app} comme en témoigne le graphe ci-dessous :


Bifurcation de Hopf

Modèles intègre-et-tire

Le principal problème du modèle de Hodgkin et Huxley est qu'il se prête mal à la simulation d'un réseau de neurones en raison de la complexité de ses équations. En particulier le modèle décrit précisément l'allure d'un spike, et durant ce laps de temps le potentiel croît et décroît brusquement, ce qui est très couteux dans les résolutions numériques à pas adaptatif. En fin de compte on connait l'allure d'un spike qui est sensiblement toujours de la même forme. Ce qui nous intéresse surtout c'est de savoir quand ils sont émis. C'est pourquoi d'autres modèles mathématiques plus simples ont été proposés de façon à décrire au mieux le comportement de certains types de neurones tout en gardant une complexité raisonnable. C'est le cas d'une famille de modèles appelés intègre-et-tire qui ont l'avantage d'être mathématiquement exploitables et capables de reproduire une grande variété de comportements du neurone.

La forme générale de ces modèles est du type:

Interprétation graphique des paramètres du modèle IF
\left\{\begin{array}{l}\dot v=F(v)-w+I_{app}\\\dot w=a(bv-w)\\si~v\geqslant v_{seuil}\\ alors~ v\leftarrow v_{reset}~et~w\leftarrow w+d\end{array}\right.


v représente le potentiel membranaire du neurone et w représente une variable d'adaptation (comme n précédemment). Avec ce modèle on considère qu'un potentiel d'action est émis dès lors que le potentiel dépasse un certain seuil v_{seuil} (égal à 30 mA par exemple, comme ci-contre), puis on réinitialise les variables. De cette façon on simule que ce qui est nécessaire (la forme du spike ne nous intéressant pas), le potentiel est "tronqué" mais on connait l'instant de spike, c'est l'essentiel.

Problématique : Tout le problème réside dans le choix du seuil, c'est-à-dire fixer une limite au delà de laquelle on peut prétendre qu'un potentiel d'action sera émis. Tel que le seuil (fixe) est défini dans le modèle intègre-et-tire, cette frontière, dans le plan de phase, est une droite verticale v=v_{seuil}. Ainsi dès qu'une trajectoire dans le plan de phase franchit cette droite, on enregistre l'instant de spike, on réinitialise les variables et on continue la simulation. Or il n'y a aucune raison biologique qui justifierait que ce seuil soit fixe. C'est à partir de ce constat que nous avons élaboré un nouveau type de seuil, qui constitue la valeur ajoutée de cette IRL.


Travail réalisé

Recherche d'un seuil plus naturel

Dans le but de visualiser l'allure de la véritable frontière séparant la zone de spike et la zone de non-spike, j'ai discrétisé le plan de phase en une grille 500x500, et pour chacune de ces 25 000 conditions initiales (V(0),n(0)), j'ai simulé le modèle de Hodgkin-Huxley, et j'ai déterminé \max_{t}V_i(t). En fonction de cette dernière valeur j'ai coloré dans le plan le point (V(0),n(0)) suivant une échelle de couleurs (du bleu si le max est faible, autrement dit s'il n'y a pas de spike, vers le rouge sinon). On obtient la coloration suivante (ainsi que les nullclines en vert et jaune) :


Visualisation de l'allure du seuil


On voit alors clairement qu'un seuil du type v=\text{constante} ne correspond pas à la réalité, il semblerait que la fonction seuil possède à peu près le même comportement que la V-nullcline (sauf dans les potentiels faibles). J'ai donc approximer la délimitation par une droite affine, tangente au point d'inflexion de la courbe. Ce nouveau seuil n'est donc plus fixe mais régit par l'équation n=0.0188V+1.4322.

Introduction de ce nouveau seuil au modèle de Hodgkin-Huxley

Simulation du modèle de Hodgkin-Huxley avec et sans seuil

Si on simule de nouveau le modèle de Hodgkin-Huxley réduit muni de ce seuil et soumis à un courant aléatoire :


\left\{\begin{array}{l}
      \frac{dV}{dt}=f(V,n)+I_{alea}(V,t)
\\\frac{dn}{dt}= g(V,n)\\
si~n<0.0188V+1.4322\\ alors~ v\leftarrow v_{reset}\\et~n\leftarrow n_{reset}
\end{array}\right.


on s'assure que l'introduction du seuil ne modifie pas les temps de spike, et effectivement les potentiels tronqués (en rouge) sont bien syncronisés avec ceux d'origine (en bleu) comme illustré ci-contre. On a donc un gain en temps de simulation puisqu'on ne simule pas la forme du spike coûteuse (autrement dit dans le plan de phase on ne simule jamais sous la droite n=0.0188V+1.4322).

Extraction d'un modèle intègre-et-tire à partir de celui de Hodgkin-Huxley

  • Pour gagner davantage en temps de simulation j'ai ensuite cherché à approcher les équations de Hodgkin et Huxley par des équations plus simples. J'ai donc essayé dans un premier temps d'approcher au sens des moindre carrés les nullclines par des polynômes, dans la zone qui nous intéresse (i.e au dessus de la droite seuil) soit la partie convexe pour la V-nullcline sur la figure. Problème : en simulant le système avec ces nouvelles équations f(V,n)=1,27.10^{-3}V^2+0,16V+5,21-n+I_{app} et g(V,n)=0,0159V+1,343-n, je n'obtenais pas du tout le même comportement, et pour cause : la seule approximation des nullclines ne peut rendre compte de la dynamique du système puisqu'elles caractérisent seulement les points où la dérivée s'annule.


  • J'ai donc ensuite déterminer des surfaces quadriques z=\alpha V^2+\beta V+\gamma+\mu n approchant au mieux les surfaces f(V,n) et g(V,n) pour la norme de L^2. Problème : les surfaces z=f(V,n) et z=g(V,n) sont maintenant bien approchées, mais les nullclines résultantes (obtenues en coupant par le plan z=0) n'ont plus du tout l'allure d'origine et elles se coupent vers -40 mV, le potentiel de repos ne correspond donc plus du tout, on perd la dynamique du système.


  • Pour pallier à ce problème il faudrait effectuer une minimisation sous contrainte, c'est-à-dire approximer au mieux les surfaces tout en conservant des nullclines pertinentes.



Discussion : application à des données biologiques

Parallèlement à l'extraction du modèle intègre-et-tire, je préparais les algorithmes visant à confronter celui-ci à des données biologiques, issues d'une expérience sur un vrai neurone soumis à un courant donné en entrée pendant 60 secondes. Ces données ont fait l'objet d'un concours proposé par l'INCF dont le but du jeu était le suivant : les participants disposaient du courant sur toute la durée de l'expérience, mais seulement de la réponse du neurone biologique sur les 40 secondes. Ils pouvaient alors ajuster au mieux leur modèle de façon à le faire correspondre le plus précisément possible à la réponse du neurone ; puis (et c'est là où ça devient intéressant) prédire le comportement du neurone sur les 20 secondes restantes. En confrontant alors les résultats des participants au réel comportement du neurone ils étaient en mesure de désigner le meilleur modèle mathématiques décrivant l'activité de ce neurone. Nous comptions donc soumettre notre modèle intègre-et-tire muni du nouveau seuil à ce challenge et analyser les résultats.


Principe de simulation et de prédiction de spikes

Conclusions

Après avoir découvert deux des modèles les plus populaires décrivant le comportement neuronal, à savoir le modèle de Hodgkin et Huxley et le modèle intègre-et-tire, j'ai travaillé plus spécifiquement sur une des caractéristiques des modèles intègre-et-tire : l'équation du seuil. Les avantages de ce nouveau type de seuil sont multiples :

  • Biologiquement plus plausible
  • Gain de temps en simulation
  • Amélioration du pouvoir prédictif

Nous avons pu vérifier les deux premiers points, en revanche pour quantifier l'amélioration du pouvoir prédictif il faudrait dans des perspectives futures :

  • Terminer l'extraction du modèle intègre-et-tire à partir de celui de Hodgkin-Huxley (en exploitant la minimisation sous contrainte)
  • Soumettre ce modèle aux données biologiques du challenge et analyser statistiquement les résultats
  • Appliquer ce nouveau type de seuil aux modèles déjà existants et l'étendre aux réseaux de neurones

Références

[1] Mathematical Foundations of Neuroscience, G.Bard Ermentrout & David H.Terman, 2010.

[2] Simple Model of Spiking Neurons, Eugene M.Izhikevich, 2003.

[3] Which Model to Use for Cortical Spiking Neurons, Eugene M. Izhikevich, 2004.

[4] Made to order spiking neuron model equipped with a multi-timescale adaptive threshold, Ryota Kobasyashi, Yasuhiro Tsubo & Shigeru Shinomoto, 2009.

[5] A benchmark test for a quantitative assessment of simple neuron models, Renaud Jolivet, Ryota Kobayashi, Alexander Rauch, Richard Naud, Shigeru Shinomoto & Wulfram Gerstner, 2007.

[6] Adaptive Exponential Integrate-and-Fire Model as an Effective Description of Neuronal Activity, Romain Bette & Wulfram Gerstner, 2005.

[7] Importance of the Cutoff Value in the Quadratic Adaptive Integrate-and-Fire Model, Jonathan Touboul, 2009.

[8] Bifurcation Analysis of general class of Nonlinear integrate-and-fire Neurons, Jonathan Touboul, 2007.

[9] Spiking Dynamics of Bidimensional Integrate-and-Fire Neurons, Jonathan Touboul, 2009.

[10] How Good Are Neuron Models?, Wulfram Gerstner & Richard Naud, 2009.

Documents additionnels

Transparents de la soutenance

Rapport