Pyhton et modele d'Hodgkin Huxley - Modéliser un neurone

[invite0ad9d3b8]

Bonjour à tous,

Le modele d'Hodgkin et Huxley modélise le comportement électrique d'un neurone à partir de raisonnement mathématique et physique sur l'electricité, et de données expérimentales. Je ne vais pas en refaire la démonstration, ça serait un peu long ^^

J'essaie d’implémenter ce modele avec python, mais je n'y arrive pas.

Voilà mon code :

 Cliquez pour afficher

import matplotlib.pyplot as plt
import numpy as np

# TOUS LES ALPHAS

def alh(U):
    return 0.07 * np.exp(-U/20.)

def alm(U):
    return 0.1*(25.-U)/((np.exp((25.-U)/10.)) -1.)

def aln(U) :
    return 0.01*(10.-U) / ((np.exp((10.-U)/10.))-1.)

# TOUS LES BETAS

def beh(U):
    return 1./(np.exp((30.-U)/10.)+1.)

def bem(U):
    return 4.* np.exp(-U/18.)

def ben(U) :
    return 0.125 * np.exp(-U/80.)

# FONCTION PROBA OUVERTURE

def open(al, pp, be, dt=0.010):
    return dt * ((al * ( 1 - pp )) - (be * pp)) 
 

# FONCTION ITERATIVE EQUATION HH

def hh(Tmax, I=2., dt=0.010, C=1.):
    # constantes du courant de fuite
    EL = 10.6
    gL = 0.3
    
    # constantes du courant K
    gK = 36.
    EK = 12.
    
    # constantes du courant NA
    gNa = 120.
    ENa = -115.
    
    # vecteurs de stockage
    tps = [0.]
    V=[0.]
    m=[0.]
    n = [0.]
    h = [0.]

    
    # boucle sur le temps 
    
    for t in range(1,int(Tmax/dt)):
        dpn = open(aln(V[t-1]), n[t-1], ben(V[t-1]))
        dpm = open(alm(V[t-1]), m[t-1], bem(V[t-1]))
        dph = open(alh(V[t-1]), h[t-1], beh(V[t-1]))
        dV = dt * (I - ((gL * (V[t-1] - EL)) - (gK * ((n[t-1] + dpn)**4.) * (V[t-1] - EK)) - (gNa * (m[t-1]+dpm)**3. * (h[t-1]+dph) * (V[t-1]-ENa))))/C 
        n.append(n[t-1] + dpn)
        m.append(m[t-1] + dpm)
        h.append(h[t-1] + dph)
        V.append(V[t-1] + dV)
        tps.append(tps[t-1] + dt)
    print h
    return tps,V 


# PLOT
tt, Vm = hh(10.)
print tt
plt.plot(tt, Vm)
plt.show()

Dans ces conditions, le potentiel, stocké dans V, doit etre égal à 0 (conditions de repos). Or le plot qui sort est égal à zéro pendant un petit moment, puis il plong tout d'un coup. Quand je print V, la valeur augmente petit à petit, de maniere infinitesimale, puis d'emballe, devient inf, puis nan nan nan .... Quand je print le temps, tps, alors qu'il devrait etre égal à 0 ; 0,01 : 0,02 ; 0,03..., il n'est pas régulier et ressemble plutot à 0 ; 0.01 ; 0.0199999999999, 0.0300000002. Je me pose donc en premier lieu la question suivante : pourquoi le temps n'est il pas régulier ?
Il y a ensuite beaucoup de questions à poser sur le fonctionnement de mon code, mais il faudrait rentrer dans les détails pour en parler, ce que je ne compte pas faire si personne n'est intéressé pour m'aider ^^
-----

[invite4492c379]

(...)

Dans ces conditions, le potentiel, stocké dans V, doit etre égal à 0 (conditions de repos). Or le plot qui sort est égal à zéro pendant un petit moment, puis il plong tout d'un coup. Quand je print V, la valeur augmente petit à petit, de maniere infinitesimale, puis d'emballe, devient inf, puis nan nan nan .... Quand je print le temps, tps, alors qu'il devrait etre égal à 0 ; 0,01 : 0,02 ; 0,03..., il n'est pas régulier et ressemble plutot à 0 ; 0.01 ; 0.0199999999999, 0.0300000002. Je me pose donc en premier lieu la question suivante : pourquoi le temps n'est il pas régulier ?
Il y a ensuite beaucoup de questions à poser sur le fonctionnement de mon code, mais il faudrait rentrer dans les détails pour en parler, ce que je ne compte pas faire si personne n'est intéressé pour m'aider ^^

Hello,

donc sans entrer dans les détails, tu es confronté à un problème classique : la représentation des réels dans les ordinateurs.
Généralement, un réel est représenté en base 2 en virgule flottante. Pour avoir des détails tu peux te référer au standard IEEE 754. Dans ton cas 0.01 ne peut-être représenté dans ce système car 0.01 n'a pas de représentation binaire finie (d'où les fluctuations). Ce standard permet de définir plusieurs type de nombres, les réels représentables, les infinis et les NaN(=Not a Number, généralement le résultat de division par 0 par exemple).
Tes calculs dépassent les capacités de représentation du standard utilisé par python.

[polo974]

Hello,
...
Tes calculs dépassent les capacités de représentation du standard utilisé par python.

Seulement la visualisation est un peu dans les choux (10 chiffres significatifs (mal) affichés mais il y en a presque 16 d'utiles).
C'est un pb (connu) des routines d'affichage de python (ce pb semble résolu ces dernières versions...)
Mais un arrondi qui part de loin loin après la virgule peut remonter jusqu'au premier chiffre dès lors que la retenue se propage...

petit conseil annexe: change le nom de ta fonction open, c'est celui d'une fonction standard d'ouverture de fichier, ça pourrait induire en erreur...

[invite0ad9d3b8]

Salut,
Merci pour vos réponses. Cette petite imprécision est donc normale. Elle ne change rien de toute façon. Mes problèmes étaient surtout mathématique, et conceptuel au final. En y réfléchissant, j'ai mieux intégré l'équation d'Hodgkin Huxley, et mon programme fonctionne finalemlent ! Le PA se génère tout seul, à V=0, même si c'est normalement au repos, d'après ce que j'ai compris cette activation automatique est normale. J'ai donc obtenu le potentiel d'action d'un neurone de calmar géant, avec émotion ^^
Pour ceux que ça interesse :

import matplotlib.pyplot as plt
import numpy as np

# TOUS LES ALPHAS

def alh(U):
    return 0.07 * np.exp(-U/20.)

def alm(U):
    return 0.1*(25.-U)/((np.exp((25.-U)/10.)) -1.)

def aln(U) :
    return 0.01*(10.-U) / ((np.exp((10.-U)/10.))-1.)

# TOUS LES BETAS

def beh(U):
    return 1./(np.exp((30.-U)/10.)+1.)

def bem(U):
    return 4.* np.exp(-U/18.)

def ben(U) :
    return 0.125 * np.exp(-U/80.)

# FONCTION PROBA OUVERTURE

def ouverture(al, pp, be):
    return (al * ( 1 - pp )) - (be * pp)
 

# FONCTION ITERATIVE EQUATION HH

def hh(Tmax, I=2., dt=0.010, C=1.):
    # constantes du courant de fuite
    EL = 10.6
    gL = 0.3
    
    # constantes du courant K
    gK = 36.
    EK = -12.
    
    # constantes du courant NA
    gNa = 120.
    ENa = 115.
    
    # vecteurs de stockage
    tps = [0.]
    V=[0.]
    m=[0.]
    n = [0.]
    h = [1.]
    K = [0.]
    Na = [0.]
    
    # boucle sur le temps 
    
    for t in range(1,int(Tmax/dt)):
        dpn = dt * ouverture(aln(V[t-1]), n[t-1], ben(V[t-1]))
        dpm = dt * ouverture(alm(V[t-1]), m[t-1], bem(V[t-1]))
        dph = dt * ouverture(alh(V[t-1]), h[t-1], beh(V[t-1]))
        cK = (gK * ((n[t-1])**4.) * (V[t-1] - EK))
        cNa = (gNa * (m[t-1])**3. * (h[t-1]) * (V[t-1]-ENa))
        K.append(cK)
        Na.append(cNa)
        dV = dt * (I - ((gL * (V[t-1] - EL)) + K[t-1] + Na[t-1]))/C
        n.append(n[t-1] + dpn)
        m.append(m[t-1] + dpm)
        h.append(h[t-1] + dph)
        V.append(V[t-1] + dV)
        tps.append(tps[t-1] + dt)
    return tps,V,K,Na,n,m,h 


# PLOT
tt, Vm, lK, lNa, n2, m2, h2 = hh(30.)
plt.subplot(1,3,1)
plt.plot(tt, Vm)
plt.subplot(1,3,2)
plt.plot(tt, lNa, label = 'Courant Na')
plt.plot(tt, lK, label = 'Courant K')
plt.legend()
plt.subplot(1,3,3)
plt.plot(tt,n2, label = 'n')
plt.plot(tt,m2, label = 'm')
plt.plot(tt,h2, label = 'h')
plt.legend()
plt.show()

Je joins le résultat :



A gauche le joli PA, au milieu les courants de sodium et de potassium, et à droite les probabilités d'ouverture des portes des canaux sodique et potassique.

[A voir en vidéo sur Futura]

[JPL]

Bravo : le PA est aussi beau qu'un vrai ! Voire plus