Source code for brainpy._src.dyn.neurons.hh

from functools import partial
from typing import Any, Sequence
from typing import Union, Callable, Optional

import brainpy.math as bm
from brainpy._src.context import share
from brainpy._src.dyn.base import NeuDyn, IonChaDyn
from brainpy._src.initialize import OneInit
from brainpy._src.initialize import Uniform, variable_, noise as init_noise
from brainpy._src.integrators import JointEq
from brainpy._src.integrators import odeint, sdeint
from brainpy._src.mixin import Container, TreeNode
from brainpy._src.types import ArrayType
from brainpy.check import is_initializer
from brainpy.types import Shape

__all__ = [
  'HHTypedNeuron',
  'CondNeuGroupLTC',
  'CondNeuGroup',
  'HHLTC',
  'HH',
  'MorrisLecarLTC',
  'MorrisLecar',
  'WangBuzsakiHHLTC',
  'WangBuzsakiHH'
]


[docs] class HHTypedNeuron(NeuDyn): pass
[docs] class CondNeuGroupLTC(HHTypedNeuron, Container, TreeNode): r"""Base class to model conductance-based neuron group. The standard formulation for a conductance-based model is given as .. math:: C_m {dV \over dt} = \sum_jg_j(E - V) + I_{ext} where :math:`g_j=\bar{g}_{j} M^x N^y` is the channel conductance, :math:`E` is the reversal potential, :math:`M` is the activation variable, and :math:`N` is the inactivation variable. :math:`M` and :math:`N` have the dynamics of .. math:: {dx \over dt} = \phi_x {x_\infty (V) - x \over \tau_x(V)} where :math:`x \in [M, N]`, :math:`\phi_x` is a temperature-dependent factor, :math:`x_\infty` is the steady state, and :math:`\tau_x` is the time constant. Equivalently, the above equation can be written as: .. math:: \frac{d x}{d t}=\phi_{x}\left(\alpha_{x}(1-x)-\beta_{x} x\right) where :math:`\alpha_{x}` and :math:`\beta_{x}` are rate constants. .. versionadded:: 2.1.9 Modeling the conductance-based neuron model. Parameters ---------- size : int, sequence of int The network size of this neuron group. method: str The numerical integration method. name : optional, str The neuron group name. """ def __init__( self, size: Shape, keep_size: bool = False, C: Union[float, ArrayType, Callable] = 1., A: Union[float, ArrayType, Callable] = 1e-3, V_th: Union[float, ArrayType, Callable] = 0., V_initializer: Union[Callable, ArrayType] = Uniform(-70, -60.), noise: Optional[Union[float, ArrayType, Callable]] = None, method: str = 'exp_auto', name: Optional[str] = None, mode: Optional[bm.Mode] = None, init_var: bool = True, input_var: bool = True, spk_type: Optional[type] = None, **channels ): super().__init__(size, keep_size=keep_size, mode=mode, name=name, ) # attribute for ``Container`` self.children = bm.node_dict(self.format_elements(IonChaDyn, **channels)) # parameters for neurons self.input_var = input_var self.C = C self.A = A self.V_th = V_th self.noise = init_noise(noise, self.varshape, num_vars=1) self._V_initializer = V_initializer self.spk_type = ((bm.float_ if isinstance(self.mode, bm.TrainingMode) else bm.bool) if (spk_type is None) else spk_type) # function if self.noise is None: self.integral = odeint(f=self.derivative, method=method) else: self.integral = sdeint(f=self.derivative, g=self.noise, method=method) if init_var: self.reset_state(self.mode) def derivative(self, V, t, I): # synapses I = self.sum_current_inputs(V, init=I) # channels for ch in self.nodes(level=1, include_self=False).subset(IonChaDyn).unique().values(): I = I + ch.current(V) return I / self.C def reset_state(self, batch_size=None): self.V = variable_(self._V_initializer, self.varshape, batch_size) self.spike = variable_(partial(bm.zeros, dtype=self.spk_type), self.varshape, batch_size) if self.input_var: self.input = variable_(bm.zeros, self.varshape, batch_size) for channel in self.nodes(level=1, include_self=False).subset(IonChaDyn).unique().values(): channel.reset_state(self.V.value, batch_size=batch_size)
[docs] def update(self, x=None): # inputs x = 0. if x is None else x if self.input_var: self.input += x x = self.input.value x = x * (1e-3 / self.A) # integral V = self.integral(self.V.value, share['t'], x, share['dt']) + self.sum_delta_inputs() # check whether the children channels have the correct parents. channels = self.nodes(level=1, include_self=False).subset(IonChaDyn).unique() self.check_hierarchies(self.__class__, **channels) # update channels for node in channels.values(): node(self.V.value) # update variables if self.spike.dtype == bool: self.spike.value = bm.logical_and(V >= self.V_th, self.V < self.V_th) else: self.spike.value = bm.logical_and(V >= self.V_th, self.V < self.V_th).astype(self.spike.dtype) self.V.value = V return self.spike.value
[docs] def clear_input(self): """Useful for monitoring inputs. """ if self.input_var: self.input.value = bm.zeros_like(self.input)
def return_info(self): return self.spike
[docs] class CondNeuGroup(CondNeuGroupLTC): def derivative(self, V, t, I): for ch in self.nodes(level=1, include_self=False).subset(IonChaDyn).unique().values(): I = I + ch.current(V) return I / self.C
[docs] def update(self, x=None): # inputs x = 0. if x is None else x x = self.sum_current_inputs(self.V.value, init=x) return super().update(x)
[docs] class HHLTC(NeuDyn): r"""Hodgkin–Huxley neuron model with liquid time constant. **Model Descriptions** The Hodgkin-Huxley (HH; Hodgkin & Huxley, 1952) model [1]_ for the generation of the nerve action potential is one of the most successful mathematical models of a complex biological process that has ever been formulated. The basic concepts expressed in the model have proved a valid approach to the study of bio-electrical activity from the most primitive single-celled organisms such as *Paramecium*, right through to the neurons within our own brains. Mathematically, the model is given by, .. math:: C \frac {dV} {dt} = -(\bar{g}_{Na} m^3 h (V &-E_{Na}) + \bar{g}_K n^4 (V-E_K) + g_{leak} (V - E_{leak})) + I(t) \frac {dx} {dt} &= \alpha_x (1-x) - \beta_x, \quad x\in {\rm{\{m, h, n\}}} &\alpha_m(V) = \frac {0.1(V+40)}{1-\exp(\frac{-(V + 40)} {10})} &\beta_m(V) = 4.0 \exp(\frac{-(V + 65)} {18}) &\alpha_h(V) = 0.07 \exp(\frac{-(V+65)}{20}) &\beta_h(V) = \frac 1 {1 + \exp(\frac{-(V + 35)} {10})} &\alpha_n(V) = \frac {0.01(V+55)}{1-\exp(-(V+55)/10)} &\beta_n(V) = 0.125 \exp(\frac{-(V + 65)} {80}) The illustrated example of HH neuron model please see `this notebook <../neurons/HH_model.ipynb>`_. The Hodgkin–Huxley model can be thought of as a differential equation system with four state variables, :math:`V_{m}(t),n(t),m(t)`, and :math:`h(t)`, that change with respect to time :math:`t`. The system is difficult to study because it is a nonlinear system and cannot be solved analytically. However, there are many numeric methods available to analyze the system. Certain properties and general behaviors, such as limit cycles, can be proven to exist. References ---------- .. [1] Hodgkin, Alan L., and Andrew F. Huxley. "A quantitative description of membrane current and its application to conduction and excitation in nerve." The Journal of physiology 117.4 (1952): 500. .. [2] https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model .. [3] Ashwin, Peter, Stephen Coombes, and Rachel Nicks. "Mathematical frameworks for oscillatory network dynamics in neuroscience." The Journal of Mathematical Neuroscience 6, no. 1 (2016): 1-92. **Examples** Here is a simple usage example: .. code-block:: python import brainpy as bp neu = bp.dyn.HHLTC(1) # raise input current from 4 mA to 40 mA inputs = bp.inputs.ramp_input(4, 40, 700, 100, 600,) runner = bp.DSRunner(neu, monitors=['V']) runner.run(inputs=inputs) bp.visualize.line_plot(runner.mon['ts'], runner.mon['V'], show=True) Parameters ---------- size: sequence of int, int The size of the neuron group. ENa: float, ArrayType, Initializer, callable The reversal potential of sodium. Default is 50 mV. gNa: float, ArrayType, Initializer, callable The maximum conductance of sodium channel. Default is 120 msiemens. EK: float, ArrayType, Initializer, callable The reversal potential of potassium. Default is -77 mV. gK: float, ArrayType, Initializer, callable The maximum conductance of potassium channel. Default is 36 msiemens. EL: float, ArrayType, Initializer, callable The reversal potential of learky channel. Default is -54.387 mV. gL: float, ArrayType, Initializer, callable The conductance of learky channel. Default is 0.03 msiemens. V_th: float, ArrayType, Initializer, callable The threshold of the membrane spike. Default is 20 mV. C: float, ArrayType, Initializer, callable The membrane capacitance. Default is 1 ufarad. V_initializer: ArrayType, Initializer, callable The initializer of membrane potential. m_initializer: ArrayType, Initializer, callable The initializer of m channel. h_initializer: ArrayType, Initializer, callable The initializer of h channel. n_initializer: ArrayType, Initializer, callable The initializer of n channel. method: str The numerical integration method. name: str The group name. """ def __init__( self, size: Union[int, Sequence[int]], sharding: Any = None, keep_size: bool = False, mode: bm.Mode = None, name: str = None, method: str = 'exp_auto', init_var: bool = True, # neuron parameters ENa: Union[float, ArrayType, Callable] = 50., gNa: Union[float, ArrayType, Callable] = 120., EK: Union[float, ArrayType, Callable] = -77., gK: Union[float, ArrayType, Callable] = 36., EL: Union[float, ArrayType, Callable] = -54.387, gL: Union[float, ArrayType, Callable] = 0.03, V_th: Union[float, ArrayType, Callable] = 20., C: Union[float, ArrayType, Callable] = 1.0, V_initializer: Union[Callable, ArrayType] = Uniform(-70, -60.), m_initializer: Optional[Union[Callable, ArrayType]] = None, h_initializer: Optional[Union[Callable, ArrayType]] = None, n_initializer: Optional[Union[Callable, ArrayType]] = None, ): # initialization super().__init__(size=size, sharding=sharding, keep_size=keep_size, mode=mode, name=name, method=method) # parameters self.ENa = self.init_param(ENa) self.EK = self.init_param(EK) self.EL = self.init_param(EL) self.gNa = self.init_param(gNa) self.gK = self.init_param(gK) self.gL = self.init_param(gL) self.C = self.init_param(C) self.V_th = self.init_param(V_th) # initializers self._m_initializer = is_initializer(m_initializer, allow_none=True) self._h_initializer = is_initializer(h_initializer, allow_none=True) self._n_initializer = is_initializer(n_initializer, allow_none=True) self._V_initializer = is_initializer(V_initializer) # integral self.integral = odeint(method=method, f=self.derivative) # model if init_var: self.reset_state(self.mode) # m channel # m_alpha = lambda self, V: 0.1 * (V + 40) / (1 - bm.exp(-(V + 40) / 10)) m_alpha = lambda self, V: 1. / bm.exprel(-(V + 40) / 10) m_beta = lambda self, V: 4.0 * bm.exp(-(V + 65) / 18) m_inf = lambda self, V: self.m_alpha(V) / (self.m_alpha(V) + self.m_beta(V)) dm = lambda self, m, t, V: self.m_alpha(V) * (1 - m) - self.m_beta(V) * m # h channel h_alpha = lambda self, V: 0.07 * bm.exp(-(V + 65) / 20.) h_beta = lambda self, V: 1 / (1 + bm.exp(-(V + 35) / 10)) h_inf = lambda self, V: self.h_alpha(V) / (self.h_alpha(V) + self.h_beta(V)) dh = lambda self, h, t, V: self.h_alpha(V) * (1 - h) - self.h_beta(V) * h # n channel # n_alpha = lambda self, V: 0.01 * (V + 55) / (1 - bm.exp(-(V + 55) / 10)) n_alpha = lambda self, V: 0.1 / bm.exprel(-(V + 55) / 10) n_beta = lambda self, V: 0.125 * bm.exp(-(V + 65) / 80) n_inf = lambda self, V: self.n_alpha(V) / (self.n_alpha(V) + self.n_beta(V)) dn = lambda self, n, t, V: self.n_alpha(V) * (1 - n) - self.n_beta(V) * n def reset_state(self, batch_size=None, **kwargs): self.V = self.init_variable(self._V_initializer, batch_size) if self._m_initializer is None: self.m = bm.Variable(self.m_inf(self.V.value), batch_axis=self.V.batch_axis) else: self.m = self.init_variable(self._m_initializer, batch_size) if self._h_initializer is None: self.h = bm.Variable(self.h_inf(self.V.value), batch_axis=self.V.batch_axis) else: self.h = self.init_variable(self._h_initializer, batch_size) if self._n_initializer is None: self.n = bm.Variable(self.n_inf(self.V.value), batch_axis=self.V.batch_axis) else: self.n = self.init_variable(self._n_initializer, batch_size) self.spike = self.init_variable(partial(bm.zeros, dtype=bool), batch_size) def dV(self, V, t, m, h, n, I): I = self.sum_current_inputs(V, init=I) I_Na = (self.gNa * m * m * m * h) * (V - self.ENa) n2 = n * n I_K = (self.gK * n2 * n2) * (V - self.EK) I_leak = self.gL * (V - self.EL) dVdt = (- I_Na - I_K - I_leak + I) / self.C return dVdt @property def derivative(self): return JointEq(self.dV, self.dm, self.dh, self.dn)
[docs] def update(self, x=None): t = share.load('t') dt = share.load('dt') x = 0. if x is None else x V, m, h, n = self.integral(self.V.value, self.m.value, self.h.value, self.n.value, t, x, dt) V += self.sum_delta_inputs() self.spike.value = bm.logical_and(self.V < self.V_th, V >= self.V_th) self.V.value = V self.m.value = m self.h.value = h self.n.value = n return self.spike.value
def return_info(self): return self.spike
[docs] class HH(HHLTC): r"""Hodgkin–Huxley neuron model. **Model Descriptions** The Hodgkin-Huxley (HH; Hodgkin & Huxley, 1952) model [1]_ for the generation of the nerve action potential is one of the most successful mathematical models of a complex biological process that has ever been formulated. The basic concepts expressed in the model have proved a valid approach to the study of bio-electrical activity from the most primitive single-celled organisms such as *Paramecium*, right through to the neurons within our own brains. Mathematically, the model is given by, .. math:: C \frac {dV} {dt} = -(\bar{g}_{Na} m^3 h (V &-E_{Na}) + \bar{g}_K n^4 (V-E_K) + g_{leak} (V - E_{leak})) + I(t) \frac {dx} {dt} &= \alpha_x (1-x) - \beta_x, \quad x\in {\rm{\{m, h, n\}}} &\alpha_m(V) = \frac {0.1(V+40)}{1-\exp(\frac{-(V + 40)} {10})} &\beta_m(V) = 4.0 \exp(\frac{-(V + 65)} {18}) &\alpha_h(V) = 0.07 \exp(\frac{-(V+65)}{20}) &\beta_h(V) = \frac 1 {1 + \exp(\frac{-(V + 35)} {10})} &\alpha_n(V) = \frac {0.01(V+55)}{1-\exp(-(V+55)/10)} &\beta_n(V) = 0.125 \exp(\frac{-(V + 65)} {80}) References ---------- .. [1] Hodgkin, Alan L., and Andrew F. Huxley. "A quantitative description of membrane current and its application to conduction and excitation in nerve." The Journal of physiology 117.4 (1952): 500. .. [2] https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model .. [3] Ashwin, Peter, Stephen Coombes, and Rachel Nicks. "Mathematical frameworks for oscillatory network dynamics in neuroscience." The Journal of Mathematical Neuroscience 6, no. 1 (2016): 1-92. **Examples** Here is a simple usage example: .. code-block:: python import brainpy as bp import matplotlib.pyplot as plt neu = bp.dyn.HH(1,) inputs = bp.inputs.ramp_input(4, 40, 700, 100, 600, ) runner = bp.DSRunner(neu, monitors=['V']) runner.run(inputs = inputs) plt.plot(runner.mon['ts'], runner.mon['V']) plt.plot(runner.mon.ts, inputs.value) # show input current plt.legend(['Membrane potential/mA', 'Input current/mA'], loc='upper right') plt.tight_layout() plt.show() The illustrated example of HH neuron model please see `this notebook <../neurons/HH_model.ipynb>`_. Parameters ---------- size: sequence of int, int The size of the neuron group. ENa: float, ArrayType, Initializer, callable The reversal potential of sodium. Default is 50 mV. gNa: float, ArrayType, Initializer, callable The maximum conductance of sodium channel. Default is 120 msiemens. EK: float, ArrayType, Initializer, callable The reversal potential of potassium. Default is -77 mV. gK: float, ArrayType, Initializer, callable The maximum conductance of potassium channel. Default is 36 msiemens. EL: float, ArrayType, Initializer, callable The reversal potential of learky channel. Default is -54.387 mV. gL: float, ArrayType, Initializer, callable The conductance of learky channel. Default is 0.03 msiemens. V_th: float, ArrayType, Initializer, callable The threshold of the membrane spike. Default is 20 mV. C: float, ArrayType, Initializer, callable The membrane capacitance. Default is 1 ufarad. V_initializer: ArrayType, Initializer, callable The initializer of membrane potential. m_initializer: ArrayType, Initializer, callable The initializer of m channel. h_initializer: ArrayType, Initializer, callable The initializer of h channel. n_initializer: ArrayType, Initializer, callable The initializer of n channel. method: str The numerical integration method. name: str The group name. """ def dV(self, V, t, m, h, n, I): I_Na = (self.gNa * m * m * m * h) * (V - self.ENa) n2 = n * n I_K = (self.gK * n2 * n2) * (V - self.EK) I_leak = self.gL * (V - self.EL) dVdt = (- I_Na - I_K - I_leak + I) / self.C return dVdt @property def derivative(self): return JointEq(self.dV, self.dm, self.dh, self.dn)
[docs] def update(self, x=None): x = 0. if x is None else x x = self.sum_current_inputs(self.V.value, init=x) return super().update(x)
[docs] class MorrisLecarLTC(NeuDyn): r"""The Morris-Lecar neuron model with liquid time constant. **Model Descriptions** The Morris-Lecar model [4]_ (Also known as :math:`I_{Ca}+I_K`-model) is a two-dimensional "reduced" excitation model applicable to systems having two non-inactivating voltage-sensitive conductances. This model was named after Cathy Morris and Harold Lecar, who derived it in 1981. Because it is two-dimensional, the Morris-Lecar model is one of the favorite conductance-based models in computational neuroscience. The original form of the model employed an instantaneously responding voltage-sensitive Ca2+ conductance for excitation and a delayed voltage-dependent K+ conductance for recovery. The equations of the model are: .. math:: \begin{aligned} C\frac{dV}{dt} =& - g_{Ca} M_{\infty} (V - V_{Ca})- g_{K} W(V - V_{K}) - g_{Leak} (V - V_{Leak}) + I_{ext} \\ \frac{dW}{dt} =& \frac{W_{\infty}(V) - W}{ \tau_W(V)} \end{aligned} Here, :math:`V` is the membrane potential, :math:`W` is the "recovery variable", which is almost invariably the normalized :math:`K^+`-ion conductance, and :math:`I_{ext}` is the applied current stimulus. **Model Parameters** ============= ============== ======== ======================================================= **Parameter** **Init Value** **Unit** **Explanation** ------------- -------------- -------- ------------------------------------------------------- V_Ca 130 mV Equilibrium potentials of Ca+.(mV) g_Ca 4.4 \ Maximum conductance of corresponding Ca+.(mS/cm2) V_K -84 mV Equilibrium potentials of K+.(mV) g_K 8 \ Maximum conductance of corresponding K+.(mS/cm2) V_Leak -60 mV Equilibrium potentials of leak current.(mV) g_Leak 2 \ Maximum conductance of leak current.(mS/cm2) C 20 \ Membrane capacitance.(uF/cm2) V1 -1.2 \ Potential at which M_inf = 0.5.(mV) V2 18 \ Reciprocal of slope of voltage dependence of M_inf.(mV) V3 2 \ Potential at which W_inf = 0.5.(mV) V4 30 \ Reciprocal of slope of voltage dependence of W_inf.(mV) phi 0.04 \ A temperature factor. (1/s) V_th 10 mV The spike threshold. ============= ============== ======== ======================================================= References ---------- .. [4] Lecar, Harold. "Morris-lecar model." Scholarpedia 2.10 (2007): 1333. .. [5] http://www.scholarpedia.org/article/Morris-Lecar_model .. [6] https://en.wikipedia.org/wiki/Morris%E2%80%93Lecar_model """ supported_modes = (bm.NonBatchingMode, bm.BatchingMode) def __init__( self, size: Union[int, Sequence[int]], sharding: Any = None, keep_size: bool = False, mode: bm.Mode = None, name: str = None, method: str = 'exp_auto', init_var: bool = True, # neuron parameters V_Ca: Union[float, ArrayType, Callable] = 130., g_Ca: Union[float, ArrayType, Callable] = 4.4, V_K: Union[float, ArrayType, Callable] = -84., g_K: Union[float, ArrayType, Callable] = 8., V_leak: Union[float, ArrayType, Callable] = -60., g_leak: Union[float, ArrayType, Callable] = 2., C: Union[float, ArrayType, Callable] = 20., V1: Union[float, ArrayType, Callable] = -1.2, V2: Union[float, ArrayType, Callable] = 18., V3: Union[float, ArrayType, Callable] = 2., V4: Union[float, ArrayType, Callable] = 30., phi: Union[float, ArrayType, Callable] = 0.04, V_th: Union[float, ArrayType, Callable] = 10., W_initializer: Union[Callable, ArrayType] = OneInit(0.02), V_initializer: Union[Callable, ArrayType] = Uniform(-70., -60.), ): # initialization super().__init__(size=size, sharding=sharding, keep_size=keep_size, mode=mode, name=name, method=method) # parameters self.V_Ca = self.init_param(V_Ca) self.g_Ca = self.init_param(g_Ca) self.V_K = self.init_param(V_K) self.g_K = self.init_param(g_K) self.V_leak = self.init_param(V_leak) self.g_leak = self.init_param(g_leak) self.C = self.init_param(C) self.V1 = self.init_param(V1) self.V2 = self.init_param(V2) self.V3 = self.init_param(V3) self.V4 = self.init_param(V4) self.phi = self.init_param(phi) self.V_th = self.init_param(V_th) # initializers self._W_initializer = is_initializer(W_initializer) self._V_initializer = is_initializer(V_initializer) # integral self.integral = odeint(method=method, f=self.derivative) # model if init_var: self.reset_state(self.mode) def reset_state(self, batch_or_mode=None, **kwargs): self.V = self.init_variable(self._V_initializer, batch_or_mode) self.W = self.init_variable(self._W_initializer, batch_or_mode) self.spike = self.init_variable(partial(bm.zeros, dtype=bool), batch_or_mode) def dV(self, V, t, W, I): I = self.sum_current_inputs(V, init=I) M_inf = (1 / 2) * (1 + bm.tanh((V - self.V1) / self.V2)) I_Ca = self.g_Ca * M_inf * (V - self.V_Ca) I_K = self.g_K * W * (V - self.V_K) I_Leak = self.g_leak * (V - self.V_leak) dVdt = (- I_Ca - I_K - I_Leak + I) / self.C return dVdt def dW(self, W, t, V): tau_W = 1 / (self.phi * bm.cosh((V - self.V3) / (2 * self.V4))) W_inf = (1 / 2) * (1 + bm.tanh((V - self.V3) / self.V4)) dWdt = (W_inf - W) / tau_W return dWdt @property def derivative(self): return JointEq(self.dV, self.dW)
[docs] def update(self, x=None): t = share.load('t') dt = share.load('dt') x = 0. if x is None else x V, W = self.integral(self.V, self.W, t, x, dt) V += self.sum_delta_inputs() spike = bm.logical_and(self.V < self.V_th, V >= self.V_th) self.V.value = V self.W.value = W self.spike.value = spike return spike
def return_info(self): return self.spike
[docs] class MorrisLecar(MorrisLecarLTC): r"""The Morris-Lecar neuron model. **Model Descriptions** The Morris-Lecar model [4]_ (Also known as :math:`I_{Ca}+I_K`-model) is a two-dimensional "reduced" excitation model applicable to systems having two non-inactivating voltage-sensitive conductances. This model was named after Cathy Morris and Harold Lecar, who derived it in 1981. Because it is two-dimensional, the Morris-Lecar model is one of the favorite conductance-based models in computational neuroscience. The original form of the model employed an instantaneously responding voltage-sensitive Ca2+ conductance for excitation and a delayed voltage-dependent K+ conductance for recovery. The equations of the model are: .. math:: \begin{aligned} C\frac{dV}{dt} =& - g_{Ca} M_{\infty} (V - V_{Ca})- g_{K} W(V - V_{K}) - g_{Leak} (V - V_{Leak}) + I_{ext} \\ \frac{dW}{dt} =& \frac{W_{\infty}(V) - W}{ \tau_W(V)} \end{aligned} Here, :math:`V` is the membrane potential, :math:`W` is the "recovery variable", which is almost invariably the normalized :math:`K^+`-ion conductance, and :math:`I_{ext}` is the applied current stimulus. **Model Parameters** ============= ============== ======== ======================================================= **Parameter** **Init Value** **Unit** **Explanation** ------------- -------------- -------- ------------------------------------------------------- V_Ca 130 mV Equilibrium potentials of Ca+.(mV) g_Ca 4.4 \ Maximum conductance of corresponding Ca+.(mS/cm2) V_K -84 mV Equilibrium potentials of K+.(mV) g_K 8 \ Maximum conductance of corresponding K+.(mS/cm2) V_Leak -60 mV Equilibrium potentials of leak current.(mV) g_Leak 2 \ Maximum conductance of leak current.(mS/cm2) C 20 \ Membrane capacitance.(uF/cm2) V1 -1.2 \ Potential at which M_inf = 0.5.(mV) V2 18 \ Reciprocal of slope of voltage dependence of M_inf.(mV) V3 2 \ Potential at which W_inf = 0.5.(mV) V4 30 \ Reciprocal of slope of voltage dependence of W_inf.(mV) phi 0.04 \ A temperature factor. (1/s) V_th 10 mV The spike threshold. ============= ============== ======== ======================================================= References ---------- .. [4] Lecar, Harold. "Morris-lecar model." Scholarpedia 2.10 (2007): 1333. .. [5] http://www.scholarpedia.org/article/Morris-Lecar_model .. [6] https://en.wikipedia.org/wiki/Morris%E2%80%93Lecar_model """ def dV(self, V, t, W, I): M_inf = (1 / 2) * (1 + bm.tanh((V - self.V1) / self.V2)) I_Ca = self.g_Ca * M_inf * (V - self.V_Ca) I_K = self.g_K * W * (V - self.V_K) I_Leak = self.g_leak * (V - self.V_leak) dVdt = (- I_Ca - I_K - I_Leak + I) / self.C return dVdt
[docs] def update(self, x=None): x = 0. if x is None else x x = self.sum_current_inputs(self.V.value, init=x) return super().update(x)
[docs] class WangBuzsakiHHLTC(NeuDyn): r"""Wang-Buzsaki model [9]_, an implementation of a modified Hodgkin-Huxley model with liquid time constant. Each model is described by a single compartment and obeys the current balance equation: .. math:: C_{m} \frac{d V}{d t}=-I_{\mathrm{Na}}-I_{\mathrm{K}}-I_{\mathrm{L}}-I_{\mathrm{syn}}+I_{\mathrm{app}} where :math:`C_{m}=1 \mu \mathrm{F} / \mathrm{cm}^{2}` and :math:`I_{\mathrm{app}}` is the injected current (in :math:`\mu \mathrm{A} / \mathrm{cm}^{2}` ). The leak current :math:`I_{\mathrm{L}}=g_{\mathrm{L}}\left(V-E_{\mathrm{L}}\right)` has a conductance :math:`g_{\mathrm{L}}=0.1 \mathrm{mS} / \mathrm{cm}^{2}`, so that the passive time constant :math:`\tau_{0}=C_{m} / g_{\mathrm{L}}=10 \mathrm{msec} ; E_{\mathrm{L}}=-65 \mathrm{mV}`. The spike-generating :math:`\mathrm{Na}^{+}` and :math:`\mathrm{K}^{+}` voltage-dependent ion currents :math:`\left(I_{\mathrm{Na}}\right.` and :math:`I_{\mathrm{K}}` ) are of the Hodgkin-Huxley type (Hodgkin and Huxley, 1952). The transient sodium current :math:`I_{\mathrm{Na}}=g_{\mathrm{Na}} m_{\infty}^{3} h\left(V-E_{\mathrm{Na}}\right)`, where the activation variable :math:`m` is assumed fast and substituted by its steady-state function :math:`m_{\infty}=\alpha_{m} /\left(\alpha_{m}+\beta_{m}\right)` ; :math:`\alpha_{m}(V)=-0.1(V+35) /(\exp (-0.1(V+35))-1), \beta_{m}(V)=4 \exp (-(V+60) / 18)`. The inactivation variable :math:`h` obeys a first-order kinetics: .. math:: \frac{d h}{d t}=\phi\left(\alpha_{h}(1-h)-\beta_{h} h\right) where :math:`\alpha_{h}(V)=0.07 \exp (-(V+58) / 20)` and :math:`\beta_{h}(V)=1 /(\exp (-0.1(V+28)) +1) \cdot g_{\mathrm{Na}}=35 \mathrm{mS} / \mathrm{cm}^{2}` ; :math:`E_{\mathrm{Na}}=55 \mathrm{mV}, \phi=5 .` The delayed rectifier :math:`I_{\mathrm{K}}=g_{\mathrm{K}} n^{4}\left(V-E_{\mathrm{K}}\right)`, where the activation variable :math:`n` obeys the following equation: .. math:: \frac{d n}{d t}=\phi\left(\alpha_{n}(1-n)-\beta_{n} n\right) with :math:`\alpha_{n}(V)=-0.01(V+34) /(\exp (-0.1(V+34))-1)` and :math:`\beta_{n}(V)=0.125\exp (-(V+44) / 80)` ; :math:`g_{\mathrm{K}}=9 \mathrm{mS} / \mathrm{cm}^{2}`, and :math:`E_{\mathrm{K}}=-90 \mathrm{mV}`. References ---------- .. [9] Wang, X.J. and Buzsaki, G., (1996) Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. Journal of neuroscience, 16(20), pp.6402-6413. **Examples** Here is a simple usage example: .. code-block:: python import brainpy as bp import matplotlib.pyplot as plt neu = bp.dyn.WangBuzsakiHHLTC(1, ) inputs = bp.inputs.ramp_input(.1, 1, 700, 100, 600, ) runner = bp.DSRunner(neu, monitors=['V']) runner.run(inputs=inputs) plt.plot(runner.mon['ts'], runner.mon['V']) plt.legend(['Membrane potential/mA', loc='upper right') plt.tight_layout() plt.show() Parameters ---------- size: sequence of int, int The size of the neuron group. ENa: float, ArrayType, Initializer, callable The reversal potential of sodium. Default is 50 mV. gNa: float, ArrayType, Initializer, callable The maximum conductance of sodium channel. Default is 120 msiemens. EK: float, ArrayType, Initializer, callable The reversal potential of potassium. Default is -77 mV. gK: float, ArrayType, Initializer, callable The maximum conductance of potassium channel. Default is 36 msiemens. EL: float, ArrayType, Initializer, callable The reversal potential of learky channel. Default is -54.387 mV. gL: float, ArrayType, Initializer, callable The conductance of learky channel. Default is 0.03 msiemens. V_th: float, ArrayType, Initializer, callable The threshold of the membrane spike. Default is 20 mV. C: float, ArrayType, Initializer, callable The membrane capacitance. Default is 1 ufarad. phi: float, ArrayType, Initializer, callable The temperature regulator constant. V_initializer: ArrayType, Initializer, callable The initializer of membrane potential. h_initializer: ArrayType, Initializer, callable The initializer of h channel. n_initializer: ArrayType, Initializer, callable The initializer of n channel. method: str The numerical integration method. name: str The group name. """ def __init__( self, size: Union[int, Sequence[int]], sharding: Any = None, keep_size: bool = False, mode: bm.Mode = None, name: str = None, method: str = 'exp_auto', init_var: bool = True, # neuron parameters ENa: Union[float, ArrayType, Callable] = 55., gNa: Union[float, ArrayType, Callable] = 35., EK: Union[float, ArrayType, Callable] = -90., gK: Union[float, ArrayType, Callable] = 9., EL: Union[float, ArrayType, Callable] = -65, gL: Union[float, ArrayType, Callable] = 0.1, V_th: Union[float, ArrayType, Callable] = 20., phi: Union[float, ArrayType, Callable] = 5.0, C: Union[float, ArrayType, Callable] = 1.0, V_initializer: Union[Callable, ArrayType] = OneInit(-65.), h_initializer: Union[Callable, ArrayType] = OneInit(0.6), n_initializer: Union[Callable, ArrayType] = OneInit(0.32), ): # initialization super().__init__(size=size, sharding=sharding, keep_size=keep_size, mode=mode, name=name, method=method) # parameters self.ENa = self.init_param(ENa) self.EK = self.init_param(EK) self.EL = self.init_param(EL) self.gNa = self.init_param(gNa) self.gK = self.init_param(gK) self.gL = self.init_param(gL) self.phi = self.init_param(phi) self.C = self.init_param(C) self.V_th = self.init_param(V_th) # initializers self._h_initializer = is_initializer(h_initializer) self._n_initializer = is_initializer(n_initializer) self._V_initializer = is_initializer(V_initializer) # integral self.integral = odeint(method=method, f=self.derivative) # model if init_var: self.reset_state(self.mode) def reset_state(self, batch_size=None): self.V = self.init_variable(self._V_initializer, batch_size) self.h = self.init_variable(self._h_initializer, batch_size) self.n = self.init_variable(self._n_initializer, batch_size) self.spike = self.init_variable(partial(bm.zeros, dtype=bool), batch_size) def m_inf(self, V): # alpha = -0.1 * (V + 35) / (bm.exp(-0.1 * (V + 35)) - 1) alpha = 1. / bm.exprel(-0.1 * (V + 35)) beta = 4. * bm.exp(-(V + 60.) / 18.) return alpha / (alpha + beta) def dh(self, h, t, V): alpha = 0.07 * bm.exp(-(V + 58) / 20) beta = 1 / (bm.exp(-0.1 * (V + 28)) + 1) dhdt = alpha * (1 - h) - beta * h return self.phi * dhdt def dn(self, n, t, V): # alpha = -0.01 * (V + 34) / (bm.exp(-0.1 * (V + 34)) - 1) alpha = 1. / bm.exprel(-0.1 * (V + 34)) beta = 0.125 * bm.exp(-(V + 44) / 80) dndt = alpha * (1 - n) - beta * n return self.phi * dndt def dV(self, V, t, h, n, I): I = self.sum_current_inputs(V, init=I) INa = self.gNa * self.m_inf(V) ** 3 * h * (V - self.ENa) IK = self.gK * n ** 4 * (V - self.EK) IL = self.gL * (V - self.EL) dVdt = (- INa - IK - IL + I) / self.C return dVdt @property def derivative(self): return JointEq(self.dV, self.dh, self.dn)
[docs] def update(self, x=None): t = share.load('t') dt = share.load('dt') x = 0. if x is None else x V, h, n = self.integral(self.V, self.h, self.n, t, x, dt) V += self.sum_delta_inputs() self.spike.value = bm.logical_and(self.V < self.V_th, V >= self.V_th) self.V.value = V self.h.value = h self.n.value = n return self.spike.value
def return_info(self): return self.spike
[docs] class WangBuzsakiHH(WangBuzsakiHHLTC): r"""Wang-Buzsaki model [9]_, an implementation of a modified Hodgkin-Huxley model. Each model is described by a single compartment and obeys the current balance equation: .. math:: C_{m} \frac{d V}{d t}=-I_{\mathrm{Na}}-I_{\mathrm{K}}-I_{\mathrm{L}}-I_{\mathrm{syn}}+I_{\mathrm{app}} where :math:`C_{m}=1 \mu \mathrm{F} / \mathrm{cm}^{2}` and :math:`I_{\mathrm{app}}` is the injected current (in :math:`\mu \mathrm{A} / \mathrm{cm}^{2}` ). The leak current :math:`I_{\mathrm{L}}=g_{\mathrm{L}}\left(V-E_{\mathrm{L}}\right)` has a conductance :math:`g_{\mathrm{L}}=0.1 \mathrm{mS} / \mathrm{cm}^{2}`, so that the passive time constant :math:`\tau_{0}=C_{m} / g_{\mathrm{L}}=10 \mathrm{msec} ; E_{\mathrm{L}}=-65 \mathrm{mV}`. The spike-generating :math:`\mathrm{Na}^{+}` and :math:`\mathrm{K}^{+}` voltage-dependent ion currents :math:`\left(I_{\mathrm{Na}}\right.` and :math:`I_{\mathrm{K}}` ) are of the Hodgkin-Huxley type (Hodgkin and Huxley, 1952). The transient sodium current :math:`I_{\mathrm{Na}}=g_{\mathrm{Na}} m_{\infty}^{3} h\left(V-E_{\mathrm{Na}}\right)`, where the activation variable :math:`m` is assumed fast and substituted by its steady-state function :math:`m_{\infty}=\alpha_{m} /\left(\alpha_{m}+\beta_{m}\right)` ; :math:`\alpha_{m}(V)=-0.1(V+35) /(\exp (-0.1(V+35))-1), \beta_{m}(V)=4 \exp (-(V+60) / 18)`. The inactivation variable :math:`h` obeys a first-order kinetics: .. math:: \frac{d h}{d t}=\phi\left(\alpha_{h}(1-h)-\beta_{h} h\right) where :math:`\alpha_{h}(V)=0.07 \exp (-(V+58) / 20)` and :math:`\beta_{h}(V)=1 /(\exp (-0.1(V+28)) +1) \cdot g_{\mathrm{Na}}=35 \mathrm{mS} / \mathrm{cm}^{2}` ; :math:`E_{\mathrm{Na}}=55 \mathrm{mV}, \phi=5 .` The delayed rectifier :math:`I_{\mathrm{K}}=g_{\mathrm{K}} n^{4}\left(V-E_{\mathrm{K}}\right)`, where the activation variable :math:`n` obeys the following equation: .. math:: \frac{d n}{d t}=\phi\left(\alpha_{n}(1-n)-\beta_{n} n\right) with :math:`\alpha_{n}(V)=-0.01(V+34) /(\exp (-0.1(V+34))-1)` and :math:`\beta_{n}(V)=0.125\exp (-(V+44) / 80)` ; :math:`g_{\mathrm{K}}=9 \mathrm{mS} / \mathrm{cm}^{2}`, and :math:`E_{\mathrm{K}}=-90 \mathrm{mV}`. References ---------- .. [9] Wang, X.J. and Buzsaki, G., (1996) Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. Journal of neuroscience, 16(20), pp.6402-6413. **Examples** Here is an example: .. code-block:: python import brainpy as bp import matplotlib.pyplot as plt neu = bp.dyn.WangBuzsakiHH(1, ) inputs = bp.inputs.ramp_input(.1, 1, 700, 100, 600, ) runner = bp.DSRunner(neu, monitors=['V']) runner.run(inputs=inputs) plt.plot(runner.mon['ts'], runner.mon['V']) plt.legend(['Membrane potential/mA', loc='upper right') plt.tight_layout() plt.show() Parameters ---------- size: sequence of int, int The size of the neuron group. ENa: float, ArrayType, Initializer, callable The reversal potential of sodium. Default is 50 mV. gNa: float, ArrayType, Initializer, callable The maximum conductance of sodium channel. Default is 120 msiemens. EK: float, ArrayType, Initializer, callable The reversal potential of potassium. Default is -77 mV. gK: float, ArrayType, Initializer, callable The maximum conductance of potassium channel. Default is 36 msiemens. EL: float, ArrayType, Initializer, callable The reversal potential of learky channel. Default is -54.387 mV. gL: float, ArrayType, Initializer, callable The conductance of learky channel. Default is 0.03 msiemens. V_th: float, ArrayType, Initializer, callable The threshold of the membrane spike. Default is 20 mV. C: float, ArrayType, Initializer, callable The membrane capacitance. Default is 1 ufarad. phi: float, ArrayType, Initializer, callable The temperature regulator constant. V_initializer: ArrayType, Initializer, callable The initializer of membrane potential. h_initializer: ArrayType, Initializer, callable The initializer of h channel. n_initializer: ArrayType, Initializer, callable The initializer of n channel. method: str The numerical integration method. name: str The group name. """ def dV(self, V, t, h, n, I): INa = self.gNa * self.m_inf(V) ** 3 * h * (V - self.ENa) IK = self.gK * n ** 4 * (V - self.EK) IL = self.gL * (V - self.EL) dVdt = (- INa - IK - IL + I) / self.C return dVdt
[docs] def update(self, x=None): x = 0. if x is None else x x = self.sum_current_inputs(self.V.value, init=x) return super().update(x)