diff --git a/README.md b/README.md index c95109027..fec353b71 100644 --- a/README.md +++ b/README.md @@ -9,9 +9,7 @@ LICENSE Documentation PyPI version - Linux CI - Windows CI - MacOS CI + Continuous Integration

diff --git a/brainpy/_src/dyn/neurons/hh.py b/brainpy/_src/dyn/neurons/hh.py index 22a528136..0643237b9 100644 --- a/brainpy/_src/dyn/neurons/hh.py +++ b/brainpy/_src/dyn/neurons/hh.py @@ -21,7 +21,7 @@ class HHLTC(HHTypeNeuLTC): - r"""Hodgkin–Huxley neuron model. + r"""Hodgkin–Huxley neuron model with liquid time constant. **Model Descriptions** @@ -311,6 +311,176 @@ def return_for_delay(self): 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}) + + 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. + + *1. Center manifold* + + Because there are four state variables, visualizing the path in phase space can + be difficult. Usually two variables are chosen, voltage :math:`V_{m}(t)` and the + potassium gating variable :math:`n(t)`, allowing one to visualize the limit cycle. + However, one must be careful because this is an ad-hoc method of visualizing the + 4-dimensional system. This does not prove the existence of the limit cycle. + + .. image:: ../../../_static/Hodgkin_Huxley_Limit_Cycle.png + :align: center + + A better projection can be constructed from a careful analysis of the Jacobian of + the system, evaluated at the equilibrium point. Specifically, the eigenvalues of + the Jacobian are indicative of the center manifold's existence. Likewise, the + eigenvectors of the Jacobian reveal the center manifold's orientation. The + Hodgkin–Huxley model has two negative eigenvalues and two complex eigenvalues + with slightly positive real parts. The eigenvectors associated with the two + negative eigenvalues will reduce to zero as time :math:`t` increases. The remaining + two complex eigenvectors define the center manifold. In other words, the + 4-dimensional system collapses onto a 2-dimensional plane. Any solution + starting off the center manifold will decay towards the *center manifold*. + Furthermore, the limit cycle is contained on the center manifold. + + *2. Bifurcations* + + If the injected current :math:`I` were used as a bifurcation parameter, then the + Hodgkin–Huxley model undergoes a Hopf bifurcation. As with most neuronal models, + increasing the injected current will increase the firing rate of the neuron. + One consequence of the Hopf bifurcation is that there is a minimum firing rate. + This means that either the neuron is not firing at all (corresponding to zero + frequency), or firing at the minimum firing rate. Because of the all-or-none + principle, there is no smooth increase in action potential amplitude, but + rather there is a sudden "jump" in amplitude. The resulting transition is + known as a `canard `_. + + .. image:: ../../../_static/Hodgkins_Huxley_bifurcation_by_I.gif + :align: center + + The following image shows the bifurcation diagram of the Hodgkin–Huxley model + as a function of the external drive :math:`I` [3]_. The green lines show the amplitude + of a stable limit cycle and the blue lines indicate unstable limit-cycle behaviour, + both born from Hopf bifurcations. The solid red line shows the stable fixed point + and the black line shows the unstable fixed point. + + .. image:: ../../../_static/Hodgkin_Huxley_bifurcation.png + :align: center + + **Model Examples** + + .. plot:: + :include-source: True + + >>> import brainpy as bp + >>> group = bp.neurons.HH(2) + >>> runner = bp.DSRunner(group, monitors=['V'], inputs=('input', 10.)) + >>> runner.run(200.) + >>> bp.visualize.line_plot(runner.mon.ts, runner.mon.V, show=True) + + .. plot:: + :include-source: True + + >>> import brainpy as bp + >>> import brainpy.math as bm + >>> import matplotlib.pyplot as plt + >>> + >>> group = bp.neurons.HH(2) + >>> + >>> I1 = bp.inputs.spike_input(sp_times=[500., 550., 1000, 1030, 1060, 1100, 1200], sp_lens=5, sp_sizes=5., duration=2000, ) + >>> I2 = bp.inputs.spike_input(sp_times=[600., 900, 950, 1500], sp_lens=5, sp_sizes=5., duration=2000, ) + >>> I1 += bp.math.random.normal(0, 3, size=I1.shape) + >>> I2 += bp.math.random.normal(0, 3, size=I2.shape) + >>> I = bm.stack((I1, I2), axis=-1) + >>> + >>> runner = bp.DSRunner(group, monitors=['V'], inputs=('input', I, 'iter')) + >>> runner.run(2000.) + >>> + >>> fig, gs = bp.visualize.get_figure(1, 1, 3, 8) + >>> fig.add_subplot(gs[0, 0]) + >>> plt.plot(runner.mon.ts, runner.mon.V[:, 0]) + >>> plt.plot(runner.mon.ts, runner.mon.V[:, 1] + 130) + >>> plt.xlim(10, 2000) + >>> plt.xticks([]) + >>> plt.yticks([]) + >>> 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. + 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. + + 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. + """ def dV(self, V, t, m, h, n, I): I_Na = (self.gNa * m ** 3.0 * h) * (V - self.ENa) I_K = (self.gK * n ** 4.0) * (V - self.EK) @@ -330,7 +500,7 @@ def update(self, x=None): class MorrisLecarLTC(HHTypeNeuLTC): - r"""The Morris-Lecar neuron model. + r"""The Morris-Lecar neuron model with liquid time constant. **Model Descriptions** @@ -506,6 +676,78 @@ def return_for_delay(self): 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 Examples** + + .. plot:: + :include-source: True + + >>> import brainpy as bp + >>> + >>> group = bp.neurons.MorrisLecar(1) + >>> runner = bp.DSRunner(group, monitors=['V', 'W'], inputs=('input', 100.)) + >>> runner.run(1000) + >>> + >>> fig, gs = bp.visualize.get_figure(2, 1, 3, 8) + >>> fig.add_subplot(gs[0, 0]) + >>> bp.visualize.line_plot(runner.mon.ts, runner.mon.W, ylabel='W') + >>> fig.add_subplot(gs[1, 0]) + >>> bp.visualize.line_plot(runner.mon.ts, runner.mon.V, ylabel='V', show=True) + + + **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) @@ -532,7 +774,7 @@ def update(self, x=None): class WangBuzsakiModelLTC(HHTypeNeuLTC): - r"""Wang-Buzsaki model [9]_, an implementation of a modified Hodgkin-Huxley model. + 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: @@ -722,6 +964,89 @@ def return_for_delay(self): return self.spike class WangBuzsakiModel(WangBuzsakiModelLTC): + 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}`. + + + 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. + + 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. + + """ def m_inf(self, V): alpha = -0.1 * (V + 35) / (bm.exp(-0.1 * (V + 35)) - 1) beta = 4. * bm.exp(-(V + 60.) / 18.) diff --git a/brainpy/_src/dyn/neurons/lif.py b/brainpy/_src/dyn/neurons/lif.py index ad999cd7f..9b6e8db35 100644 --- a/brainpy/_src/dyn/neurons/lif.py +++ b/brainpy/_src/dyn/neurons/lif.py @@ -39,6 +39,10 @@ 'GifLTC', 'GifRef', 'GifRefLTC', + 'Izhikevich', + 'IzhikevichLTC', + 'IzhikevichRef', + 'IzhikevichRefLTC', ] @@ -614,8 +618,7 @@ def update(self, x=None): x += out(self.V.value) super().update(x) -ExpIF.__doc__ = ExpIFLTC.__doc__ % ('') -ExpIFLTC.__doc__ = ExpIFLTC.__doc__ % (ltc_doc) + class ExpIFRefLTC(ExpIFLTC): def __init__( @@ -739,8 +742,86 @@ def update(self, x=None): x += out(self.V.value) super().update(x) +ExpIF.__doc__ = ExpIFLTC.__doc__ % ('') +ExpIFRefLTC.__doc__ = ExpIFLTC.__doc__ % (ltc_doc) +ExpIFRef.__doc__ = ExpIFLTC.__doc__ % ('') +ExpIFLTC.__doc__ = ExpIFLTC.__doc__ % (ltc_doc) + class AdExIFLTC(GradNeuDyn): + r"""Adaptive exponential integrate-and-fire neuron model %s. + + **Model Descriptions** + + The **adaptive exponential integrate-and-fire model**, also called AdEx, is a + spiking neuron model with two variables [1]_ [2]_. + + .. math:: + + \begin{aligned} + \tau_m\frac{d V}{d t} &= - (V-V_{rest}) + \Delta_T e^{\frac{V-V_T}{\Delta_T}} - Rw + RI(t), \\ + \tau_w \frac{d w}{d t} &=a(V-V_{rest}) - w + \end{aligned} + + once the membrane potential reaches the spike threshold, + + .. math:: + + V \rightarrow V_{reset}, \\ + w \rightarrow w+b. + + The first equation describes the dynamics of the membrane potential and includes + an activation term with an exponential voltage dependence. Voltage is coupled to + a second equation which describes adaptation. Both variables are reset if an action + potential has been triggered. The combination of adaptation and exponential voltage + dependence gives rise to the name Adaptive Exponential Integrate-and-Fire model. + + The adaptive exponential integrate-and-fire model is capable of describing known + neuronal firing patterns, e.g., adapting, bursting, delayed spike initiation, + initial bursting, fast spiking, and regular spiking. + + **Model Examples** + + - `Examples for different firing patterns `_ + + **Model Parameters** + + ============= ============== ======== ======================================================================================================================== + **Parameter** **Init Value** **Unit** **Explanation** + ------------- -------------- -------- ------------------------------------------------------------------------------------------------------------------------ + V_rest -65 mV Resting potential. + V_reset -68 mV Reset potential after spike. + V_th -30 mV Threshold potential of spike and reset. + V_T -59.9 mV Threshold potential of generating action potential. + delta_T 3.48 \ Spike slope factor. + a 1 \ The sensitivity of the recovery variable :math:`u` to the sub-threshold fluctuations of the membrane potential :math:`v` + b 1 \ The increment of :math:`w` produced by a spike. + R 1 \ Membrane resistance. + tau 10 ms Membrane time constant. Compute by R * C. + tau_w 30 ms Time constant of the adaptation current. + tau_ref 0. ms Refractory time. + ============= ============== ======== ======================================================================================================================== + + **Model Variables** + + ================== ================= ========================================================= + **Variables name** **Initial Value** **Explanation** + ------------------ ----------------- --------------------------------------------------------- + V 0 Membrane potential. + w 0 Adaptation current. + input 0 External and synaptic input current. + spike False Flag to mark whether the neuron is spiking. + refractory False Flag to mark whether the neuron is in refractory period. + t_last_spike -1e7 Last spike time stamp. + ================== ================= ========================================================= + + **References** + + .. [1] Fourcaud-Trocmé, Nicolas, et al. "How spike generation + mechanisms determine the neuronal response to fluctuating + inputs." Journal of Neuroscience 23.37 (2003): 11628-11640. + .. [2] http://www.scholarpedia.org/article/Adaptive_exponential_integrate-and-fire_model + """ def __init__( self, size: Shape, @@ -1013,8 +1094,77 @@ def update(self, x=None): x += out(self.V.value) super().update(x) +AdExIF.__doc__ = AdExIFLTC.__doc__ % ('') +AdExIFRefLTC.__doc__ = AdExIFLTC.__doc__ % (ltc_doc) +AdExIFRef.__doc__ = AdExIFLTC.__doc__ % ('') +AdExIFLTC.__doc__ = AdExIFLTC.__doc__ % (ltc_doc) class QuaIFLTC(GradNeuDyn): + r"""Quadratic Integrate-and-Fire neuron model %s. + + **Model Descriptions** + + In contrast to physiologically accurate but computationally expensive + neuron models like the Hodgkin–Huxley model, the QIF model [1]_ seeks only + to produce **action potential-like patterns** and ignores subtleties + like gating variables, which play an important role in generating action + potentials in a real neuron. However, the QIF model is incredibly easy + to implement and compute, and relatively straightforward to study and + understand, thus has found ubiquitous use in computational neuroscience. + + .. math:: + + \tau \frac{d V}{d t}=c(V-V_{rest})(V-V_c) + RI(t) + + where the parameters are taken to be :math:`c` =0.07, and :math:`V_c = -50 mV` (Latham et al., 2000). + + **Model Examples** + + .. plot:: + :include-source: True + + >>> import brainpy as bp + >>> + >>> group = bp.neurons.QuaIF(1,) + >>> + >>> runner = bp.DSRunner(group, monitors=['V'], inputs=('input', 20.)) + >>> runner.run(duration=200.) + >>> bp.visualize.line_plot(runner.mon.ts, runner.mon.V, show=True) + + + **Model Parameters** + + ============= ============== ======== ======================================================================================================================== + **Parameter** **Init Value** **Unit** **Explanation** + ------------- -------------- -------- ------------------------------------------------------------------------------------------------------------------------ + V_rest -65 mV Resting potential. + V_reset -68 mV Reset potential after spike. + V_th -30 mV Threshold potential of spike and reset. + V_c -50 mV Critical voltage for spike initiation. Must be larger than V_rest. + c .07 \ Coefficient describes membrane potential update. Larger than 0. + R 1 \ Membrane resistance. + tau 10 ms Membrane time constant. Compute by R * C. + tau_ref 0 ms Refractory period length. + ============= ============== ======== ======================================================================================================================== + + **Model Variables** + + ================== ================= ========================================================= + **Variables name** **Initial Value** **Explanation** + ------------------ ----------------- --------------------------------------------------------- + V 0 Membrane potential. + input 0 External and synaptic input current. + spike False Flag to mark whether the neuron is spiking. + refractory False Flag to mark whether the neuron is in refractory period. + t_last_spike -1e7 Last spike time stamp. + ================== ================= ========================================================= + + **References** + + .. [1] P. E. Latham, B.J. Richmond, P. Nelson and S. Nirenberg + (2000) Intrinsic dynamics in neuronal networks. I. Theory. + J. Neurophysiology 83, pp. 808–827. + """ def __init__( self, size: Shape, @@ -1237,7 +1387,88 @@ def update(self, x=None): super().update(x) +QuaIF.__doc__ = QuaIFLTC.__doc__ % ('') +QuaIFRefLTC.__doc__ = QuaIFLTC.__doc__ % (ltc_doc) +QuaIFRef.__doc__ = QuaIFLTC.__doc__ % ('') +QuaIFLTC.__doc__ = QuaIFLTC.__doc__ % (ltc_doc) + + class AdQuaIFLTC(GradNeuDyn): + r"""Adaptive quadratic integrate-and-fire neuron model %s. + + **Model Descriptions** + + The adaptive quadratic integrate-and-fire neuron model [1]_ is given by: + + .. math:: + + \begin{aligned} + \tau_m \frac{d V}{d t}&=c(V-V_{rest})(V-V_c) - w + I(t), \\ + \tau_w \frac{d w}{d t}&=a(V-V_{rest}) - w, + \end{aligned} + + once the membrane potential reaches the spike threshold, + + .. math:: + + V \rightarrow V_{reset}, \\ + w \rightarrow w+b. + + **Model Examples** + + .. plot:: + :include-source: True + + >>> import brainpy as bp + >>> group = bp.neurons.AdQuaIF(1, ) + >>> runner = bp.DSRunner(group, monitors=['V', 'w'], inputs=('input', 30.)) + >>> runner.run(300) + >>> fig, gs = bp.visualize.get_figure(2, 1, 3, 8) + >>> fig.add_subplot(gs[0, 0]) + >>> bp.visualize.line_plot(runner.mon.ts, runner.mon.V, ylabel='V') + >>> fig.add_subplot(gs[1, 0]) + >>> bp.visualize.line_plot(runner.mon.ts, runner.mon.w, ylabel='w', show=True) + + **Model Parameters** + + ============= ============== ======== ======================================================= + **Parameter** **Init Value** **Unit** **Explanation** + ------------- -------------- -------- ------------------------------------------------------- + V_rest -65 mV Resting potential. + V_reset -68 mV Reset potential after spike. + V_th -30 mV Threshold potential of spike and reset. + V_c -50 mV Critical voltage for spike initiation. Must be larger + than :math:`V_{rest}`. + a 1 \ The sensitivity of the recovery variable :math:`u` to + the sub-threshold fluctuations of the membrane + potential :math:`v` + b .1 \ The increment of :math:`w` produced by a spike. + c .07 \ Coefficient describes membrane potential update. + Larger than 0. + tau 10 ms Membrane time constant. + tau_w 10 ms Time constant of the adaptation current. + ============= ============== ======== ======================================================= + + **Model Variables** + + ================== ================= ========================================================== + **Variables name** **Initial Value** **Explanation** + ------------------ ----------------- ---------------------------------------------------------- + V 0 Membrane potential. + w 0 Adaptation current. + input 0 External and synaptic input current. + spike False Flag to mark whether the neuron is spiking. + t_last_spike -1e7 Last spike time stamp. + ================== ================= ========================================================== + + **References** + + .. [1] Izhikevich, E. M. (2004). Which model to use for cortical spiking + neurons?. IEEE transactions on neural networks, 15(5), 1063-1070. + .. [2] Touboul, Jonathan. "Bifurcation analysis of a general class of + nonlinear integrate-and-fire neurons." SIAM Journal on Applied + Mathematics 68, no. 4 (2008): 1045-1079. + """ def __init__( self, size: Shape, @@ -1504,7 +1735,93 @@ def update(self, x=None): super().update(x) +AdQuaIF.__doc__ = AdQuaIFLTC.__doc__ % ('') +AdQuaIFRefLTC.__doc__ = AdQuaIFLTC.__doc__ % (ltc_doc) +AdQuaIFRef.__doc__ = AdQuaIFLTC.__doc__ % ('') +AdQuaIFLTC.__doc__ = AdQuaIFLTC.__doc__ % (ltc_doc) + + class GifLTC(GradNeuDyn): + r"""Generalized Integrate-and-Fire model %s. + + **Model Descriptions** + + The generalized integrate-and-fire model [1]_ is given by + + .. math:: + + &\frac{d I_j}{d t} = - k_j I_j + + &\frac{d V}{d t} = ( - (V - V_{rest}) + R\sum_{j}I_j + RI) / \tau + + &\frac{d V_{th}}{d t} = a(V - V_{rest}) - b(V_{th} - V_{th\infty}) + + When :math:`V` meet :math:`V_{th}`, Generalized IF neuron fires: + + .. math:: + + &I_j \leftarrow R_j I_j + A_j + + &V \leftarrow V_{reset} + + &V_{th} \leftarrow max(V_{th_{reset}}, V_{th}) + + Note that :math:`I_j` refers to arbitrary number of internal currents. + + **Model Examples** + + - `Detailed examples to reproduce different firing patterns `_ + + **Model Parameters** + + ============= ============== ======== ==================================================================== + **Parameter** **Init Value** **Unit** **Explanation** + ------------- -------------- -------- -------------------------------------------------------------------- + V_rest -70 mV Resting potential. + V_reset -70 mV Reset potential after spike. + V_th_inf -50 mV Target value of threshold potential :math:`V_{th}` updating. + V_th_reset -60 mV Free parameter, should be larger than :math:`V_{reset}`. + R 20 \ Membrane resistance. + tau 20 ms Membrane time constant. Compute by :math:`R * C`. + a 0 \ Coefficient describes the dependence of + :math:`V_{th}` on membrane potential. + b 0.01 \ Coefficient describes :math:`V_{th}` update. + k1 0.2 \ Constant pf :math:`I1`. + k2 0.02 \ Constant of :math:`I2`. + R1 0 \ Free parameter. + Describes dependence of :math:`I_1` reset value on + :math:`I_1` value before spiking. + R2 1 \ Free parameter. + Describes dependence of :math:`I_2` reset value on + :math:`I_2` value before spiking. + A1 0 \ Free parameter. + A2 0 \ Free parameter. + ============= ============== ======== ==================================================================== + + **Model Variables** + + ================== ================= ========================================================= + **Variables name** **Initial Value** **Explanation** + ------------------ ----------------- --------------------------------------------------------- + V -70 Membrane potential. + input 0 External and synaptic input current. + spike False Flag to mark whether the neuron is spiking. + V_th -50 Spiking threshold potential. + I1 0 Internal current 1. + I2 0 Internal current 2. + t_last_spike -1e7 Last spike time stamp. + ================== ================= ========================================================= + + **References** + + .. [1] Mihalaş, Ştefan, and Ernst Niebur. "A generalized linear + integrate-and-fire neural model produces diverse spiking + behaviors." Neural computation 21.3 (2009): 704-718. + .. [2] Teeter, Corinne, Ramakrishnan Iyer, Vilas Menon, Nathan + Gouwens, David Feng, Jim Berg, Aaron Szafer et al. "Generalized + leaky integrate-and-fire models classify multiple neuron types." + Nature communications 9, no. 1 (2018): 1-15. + """ def __init__( self, size: Shape, @@ -1828,7 +2145,79 @@ def update(self, x=None): super().update(x) +Gif.__doc__ = GifLTC.__doc__ % ('') +GifRefLTC.__doc__ = GifLTC.__doc__ % (ltc_doc) +GifRef.__doc__ = GifLTC.__doc__ % ('') +GifLTC.__doc__ = GifLTC.__doc__ % (ltc_doc) + + class IzhikevichLTC(GradNeuDyn): + r"""The Izhikevich neuron model %s. + + **Model Descriptions** + + The dynamics of the Izhikevich neuron model [1]_ [2]_ is given by: + + .. math :: + + \frac{d V}{d t} &= 0.04 V^{2}+5 V+140-u+I + + \frac{d u}{d t} &=a(b V-u) + + .. math :: + + \text{if} v \geq 30 \text{mV}, \text{then} + \begin{cases} v \leftarrow c \\ + u \leftarrow u+d \end{cases} + + **Model Examples** + + - `Detailed examples to reproduce different firing patterns `_ + + **Model Parameters** + + ============= ============== ======== ================================================================================ + **Parameter** **Init Value** **Unit** **Explanation** + ------------- -------------- -------- -------------------------------------------------------------------------------- + a 0.02 \ It determines the time scale of + the recovery variable :math:`u`. + b 0.2 \ It describes the sensitivity of the + recovery variable :math:`u` to + the sub-threshold fluctuations of the + membrane potential :math:`v`. + c -65 \ It describes the after-spike reset value + of the membrane potential :math:`v` caused by + the fast high-threshold :math:`K^{+}` + conductance. + d 8 \ It describes after-spike reset of the + recovery variable :math:`u` + caused by slow high-threshold + :math:`Na^{+}` and :math:`K^{+}` conductance. + tau_ref 0 ms Refractory period length. [ms] + V_th 30 mV The membrane potential threshold. + ============= ============== ======== ================================================================================ + + **Model Variables** + + ================== ================= ========================================================= + **Variables name** **Initial Value** **Explanation** + ------------------ ----------------- --------------------------------------------------------- + V -65 Membrane potential. + u 1 Recovery variable. + input 0 External and synaptic input current. + spike False Flag to mark whether the neuron is spiking. + refractory False Flag to mark whether the neuron is in refractory period. + t_last_spike -1e7 Last spike time stamp. + ================== ================= ========================================================= + + **References** + + .. [1] Izhikevich, Eugene M. "Simple model of spiking neurons." IEEE + Transactions on neural networks 14.6 (2003): 1569-1572. + + .. [2] Izhikevich, Eugene M. "Which model to use for cortical spiking neurons?." + IEEE transactions on neural networks 15.5 (2004): 1063-1070. + """ def __init__( self, size: Shape, @@ -1874,7 +2263,7 @@ def __init__( # initializers self._V_initializer = is_initializer(V_initializer) - self._u_initializer = is_initializer(u_initializer) + self._u_initializer = is_initializer(u_initializer, allow_none=True) # integral self.integral = odeint(method=method, f=self.derivative) @@ -1942,6 +2331,7 @@ def du(self, u, t, V): dudt = self.a * (self.b * V - u) return dudt + @property def derivative(self): return JointEq([self.dV, self.du]) @@ -2012,7 +2402,7 @@ def __init__( # initializers self._V_initializer = is_initializer(V_initializer) - self._u_initializer = is_initializer(u_initializer) + self._u_initializer = is_initializer(u_initializer, allow_none=True) # integral self.integral = odeint(method=method, f=self.derivative) @@ -2077,6 +2467,7 @@ def du(self, u, t, V): dudt = self.a * (self.b * V - u) return dudt + @property def derivative(self): return JointEq([self.dV, self.du]) @@ -2084,4 +2475,10 @@ def update(self, x=None): x = 0. if x is None else x for out in self.cur_inputs.values(): x += out(self.V.value) - super().update(x) \ No newline at end of file + super().update(x) + + +Izhikevich.__doc__ = IzhikevichLTC.__doc__ % ('') +IzhikevichRefLTC.__doc__ = IzhikevichLTC.__doc__ % (ltc_doc) +IzhikevichRef.__doc__ = IzhikevichLTC.__doc__ % ('') +IzhikevichLTC.__doc__ = IzhikevichLTC.__doc__ % (ltc_doc) diff --git a/brainpy/_src/tests/test_brainpy_deprecations.py b/brainpy/_src/tests/test_brainpy_deprecations.py index 519e0e4e8..bd23632ce 100644 --- a/brainpy/_src/tests/test_brainpy_deprecations.py +++ b/brainpy/_src/tests/test_brainpy_deprecations.py @@ -5,7 +5,7 @@ mode_deprecated_names = list(brainpy.modes.__deprecations.keys()) tools_deprecated_names = list(brainpy.tools.__deprecations.keys()) train_deprecated_names = list(brainpy.train.__deprecations.keys()) -dyn_deprecated_names = list(brainpy.dyn.__deprecations.keys()) +# dyn_deprecated_names = list(brainpy.dyn.__deprecations.keys()) intg_deprecated_names = list(brainpy.integrators.__deprecations.keys()) io_deprecated_names = list(brainpy.base.io.__deprecations.keys()) @@ -40,13 +40,13 @@ def test_brainpy_train(self, name): with self.assertWarns(DeprecationWarning): getattr(brainpy.train, name) - @parameterized.product( - name=dyn_deprecated_names - ) - def test_brainpy_dyn(self, name): - with self.assertWarns(DeprecationWarning): - getattr(brainpy.dyn, name) - + # @parameterized.product( + # name=dyn_deprecated_names + # ) + # def test_brainpy_dyn(self, name): + # with self.assertWarns(DeprecationWarning): + # getattr(brainpy.dyn, name) + # @parameterized.product( name=intg_deprecated_names ) diff --git a/brainpy/dyn/channels.py b/brainpy/dyn/channels.py index e69de29bb..f4f0d0283 100644 --- a/brainpy/dyn/channels.py +++ b/brainpy/dyn/channels.py @@ -0,0 +1,54 @@ +from brainpy._src.dyn.channels.base import ( + Ion, + IonChannel, + Calcium, + IhChannel, + CalciumChannel, + SodiumChannel, + PotassiumChannel, + LeakyChannel, +) + +from brainpy._src.dyn.channels.Ca import ( + CalciumFixed, + CalciumChannel, + CalciumDetailed, + CalciumFirstOrder, + CalciumDyna, + ICaN_IS2008, + ICaT_HM1992, + ICaT_HP1992, + ICaHT_HM1992, + ICaL_IS2008, +) + +from brainpy._src.dyn.channels.K import ( + IKDR_Ba2002, + IK_TM1991, + IK_HH1952, + IKA1_HM1992, + IKA2_HM1992, + IKK2A_HM1992, + IKK2B_HM1992, + IKNI_Ya1989, +) + +from brainpy._src.dyn.channels.IH import ( + Ih_HM1992, + Ih_De1996, +) + +from brainpy._src.dyn.channels.KCa import ( + IAHP_De1994 +) + +from brainpy._src.dyn.channels.Na import ( + INa_Ba2002, + INa_TM1991, + INa_HH1952, +) + +from brainpy._src.dyn.channels.leaky import ( + IL, + IKL, +) diff --git a/brainpy/dyn/neurons.py b/brainpy/dyn/neurons.py index 79d6c1caa..61ab26852 100644 --- a/brainpy/dyn/neurons.py +++ b/brainpy/dyn/neurons.py @@ -1,6 +1,9 @@ from brainpy._src.dyn.base import ( NeuDyn, + GradNeuDyn, + HHTypeNeu, + HHTypeNeuLTC ) from brainpy._src.dyn.neurons.lif import ( @@ -8,7 +11,46 @@ LifLTC, LifRefLTC, LifRef, + ExpIF, + ExpIFLTC, + ExpIFRefLTC, + ExpIFRef, + AdExIF, + AdExIFLTC, + AdExIFRefLTC, + AdExIFRef, + QuaIF, + QuaIFLTC, + QuaIFRefLTC, + QuaIFRef, + AdQuaIF, + AdQuaIFLTC, + AdQuaIFRefLTC, + AdQuaIFRef, + Gif, + GifLTC, + GifRefLTC, + GifRef, + Izhikevich, + IzhikevichLTC, + IzhikevichRefLTC, + IzhikevichRef ) +from brainpy._src.dyn.neurons.hh import ( + HH, + HHLTC, + MorrisLecar, + MorrisLecarLTC, + WangBuzsakiModel, + WangBuzsakiModelLTC, +) + +from brainpy._src.dyn.neurons.input import ( + InputGroup, + OutputGroup, + SpikeTimeGroup, + PoissonGroup, +)