Memristive Hodgkin-Huxley Spiking Neuron Model for Reproducing Neuron Behaviors

Abstract

The Hodgkin-Huxley (HH) spiking neuron model reproduces the dynamic characteristics of the neuron by mimicking the action potential, ionic channels, and spiking behaviors. The memristor is a nonlinear device with variable resistance. In this paper, the memristor is introduced to the HH spiking model, and the memristive Hodgkin-Huxley spiking neuron model (MHH) is presented. We experimentally compare the HH spiking model and the MHH spiking model by applying different stimuli. First, the individual current pulse is injected into the HH and MHH spiking models. The comparison between action potentials, current densities, and conductances is carried out. Second, the reverse single pulse stimulus and a series of pulse stimuli are applied to the two models. The effects of current density and action time on the production of the action potential are analyzed. Finally, the sinusoidal current stimulus acts on the two models. The various spiking behaviors are realized by adjusting the frequency of the sinusoidal stimulus. We experimentally demonstrate that the MHH spiking model generates more action potential than the HH spiking model and takes a short time to change the memductance. The reverse stimulus cannot activate the action potential in both models. The MHH spiking model performs smoother waveforms and a faster speed to return to the resting potential. The larger the external stimulus, the faster action potential generated, and the more noticeable change in conductances. Meanwhile, the MHH spiking model shows the various spiking patterns of neurons.

1. Introduction

Neurons with highly nonlinear characteristics act as the basic functional unit of receiving and propagating signals. The whole procedure of processing signals in the nerve system needs the cooperation of neurons. Some theoretical knowledge and research methods are beneficial to unveil the mechanism of information propagation in neurons. Italian scientist Camillo Golgi worked on the nervous system structure and earned the Nobel Prize for physiology and medicine in 1906 (Dröscher, 1998). In 1998, Ramon y Cajal pointed out that the neurons without directly connecting each other in the nerve system (Raviola and Mazzarello, 2011). To replicate the functions and mechanisms of neurons, we urgently need to construct the biophysical model. A variety of neuron models are emerging, and the Hodgkin-Huxley (HH) spiking neuron model is the original (Hodgkin and Huxley, 1989). Stochastic Hodgkin-Huxley Neuron Systems with the NEF is helpful to study neuron sensitivity (Chen and Li, 2010). The Hodgkin-Huxley Model with automatic parameter estimation is applied to the neuromimetic chips (Buhry et al., 2011). The space-clamped Hodgkin-Huxley model effectively inhibits the production of spikes under the injection of the noisy synaptic input (Tuckwell and Ditlevsen, 2016). The Langevin is combined with the Hodgkin-Huxley system performs accurate interspike interval (ISI) and realizes the accuracy minimal loss (Pu and Thomas, 2020). The Berger-Levy theory is introduced to the Hodgkin-Huxley model, demonstrate that the information communication between neurons is related to the presynaptic firing rate and the synchronization (Ghavami et al., 2018).

The memristor with the non-volatility and variable resistance characteristics is regarded as the fourth passive circuit element. Therefore, it becomes a hot topic in neural computing (Le et al., 2015), learning and memorizing (Sayyaparaju et al., 2018), micro-circuitry design (Berdan et al., 2014), biological synapse (Mandal and Saha, 2016), and neuron modeling (Maheshwar et al., 2014), and so on. The synaptic plasticity of biological neuronal systems can be realized by memristors and memristive crossbar in 3-D architecture to mimic the human brain (Truong et al., 2016). The memristor with hysteresis and memory characteristics is the most promising candidate for establishing the brain-like neuromorphic system (Mokhtar et al., 2017). The key features of biological neurons and synapses can be mimicked by memristors (Berdan et al., 2016; Mandal and Saha, 2016). The ion motion in neurons is represented by the electrical conductance change of a memristor (Xia and Yang, 2019). A memristor is used as a two-terminal resistor with memory (Chua, 1971; Strukov et al., 2008) performs well in storing information according to the physical laws (Yang et al., 2013). The memristor entirely avoids the data transformation bottleneck between the memory and computation (Li and Wang, 2019). The memristor crossbar array can be used to integrate the co-processor chip, which will realize machine learning algorithms and neuromorphic computing (James, 2019).

This work elaborates on the construction of the memristive Hodgkin-Huxley spiking neuron model. The mathematical expressions and the circuit of the HH spiking model are presented and analyzed in sections 2, 3. Section 4 describes the MHH spiking model and discusses the memristors used to mimic the ion channels. The comparison between two models under the different stimuli is conducted in section 5. Section 6 is the conclusion of the paper.

2. The Hodgkin-Huxley (HH) Spiking Neuron Model

The neuron cell membrane is a voltage-gated ion channel, which has high selectivity for the permeability of external and internal ions in body fluid. Only one type of ion can pass through specific channels. There involves four ionic components, sodium, potassium, calcium, and chloride. The transmembrane current depends on the rapid inward current caused by sodium and the slow outward current caused by potassium (Häusser, 2000). The ion concentration difference inside and outside of the cell is the primary driving force of neural activities. When the sodium channels are opened, the high concentration sodium flows from extracellular to intracellular, the depolarization is produced, the action potential is generated. And then, the sodium channels are closed, and the potassium channels are opened, the potassium permeates from intracellular to extracellular, the repolarization is performed. Finally, the membrane potential undergoes a hyperpolarization phase, the membrane potential shifts back to the resting potential. The above process is the generation mechanism of the action potential in a neuron.

The inside of the axon membrane is full of ionic fluids (cytoplasm), the outside of the axon membrane is filled with body fluids. The fluids (conductor) of intracellular and extracellular are separated by the axon membrane (insulator). When an insulator separates two conductors, the capacitor emerges to model the charge storage capacity. The part of the axon membrane without ion channels is equivalent to a capacitor (C_m). The axon membrane of the neuron consists of the lipid bilayer, the membrane protein, and ion channels (the upper image in Figure 1). The sodium ion channel is represented by a nonlinear conductance (g_Na), the potassium ion channel is denoted by a nonlinear conductance (g_K), and other ion channels are described as a linear conductance (g_L) (Beck et al., 2020). When the neuron is in the resting state, a potential difference is caused by the ionic concentration between the intracellular and extracellular fluids. The potential difference is called the equilibrium potential of each ion (E_Na, E_K, and E_L), which is equivalent to a driving power supply (the lower image in Figure 1).

FIGURE 1

Figure 1. The voltage-gated channel of the axon cell membrane.

When the neuron is in the resting state, there is a resting potential. Here, we choose v_rest = −65 mV as the resting potential in experiments (Hodgkin and Huxley, 1952). The V_m denotes the membrane potential, E_Na (50 mV), E_K (−70 mV), and E_L (−50 mV) represent the Nernst equilibrium potentials. When the potassium current passes through the potassium channel, the potassium current is proportional to the difference between the membrane potential and E_K (Hodgkin and Huxley, 1989; Börgers, 2017):

IK=gK(Vm−EK) (1) IK=gK(Vm-EK) (1)

Here, g_K is the potassium conductance, (V_m-E_K) is the potassium driving force. The sodium current and the leaky current are described as:

INa=gNa(Vm−ENa) (2) INa=gNa(Vm-ENa) (2)

IL=gL(Vm−EL) (3) IL=gL(Vm-EL) (3)

The ion channels are sensitive to membrane potential, which control the open and close states of channels.

In the Hodgkin-Huxley spiking model, the conductance value of each ion channel is decided by the gate-controlled variables m, n, h, and 0 ≤ m ≤ 1, 0 ≤ n ≤ 1, 0 ≤ h ≤ 1. The potassium channel depends on four active gate variables (n). The sodium channel is controlled by three active gate variables (m) and one inactive gate variable (h). The potassium conductance, the sodium conductance, and the leaky conductance are described as:

gK=gKmaxn4 (4) gK=gKmaxn4 (4)

gNa=gNamaxm3h (5) gNa=gNamaxm3h (5)

gL=gLmax (6) gL=gLmax (6)

Here, g_Kmax, g_Namax, and g_Lmax denote the maximum values of potassium, sodium, and leaky conductances, accordingly. Their values are 36, 120, 0.3 Ohm⁻¹cm⁻² (Hodgkin and Huxley, 1952, 1989). The expressions of gate-controlled variables of ion channels are written as follows:

dm/dt=1/τm(m∞−m) (7) dm/dt=1/τm(m∞-m) (7)

dn/dt=1/τn(n∞−n) (8) dn/dt=1/τn(n∞-n) (8)

dh/dt=1/τh(h∞−h) (9) dh/dt=1/τh(h∞-h) (9)

The time constants τ_m, τ_n, and τ_h change with m, n, and h, accordingly. The transition rate α characterizes the ion channels change from the close state to the open state. The transition rate β indicates the ion channels vary from the open state to the close state. m_∞, n_∞, and h_∞ are the steady-state values of the gate variables m, n, and h, accordingly (Saïgai et al., 2011). They are all the functions of the membrane potential. Their expressions are:

m∞=αm/(αm+βm) (10) m∞=αm/(αm+βm) (10)

n∞=αn/(αn+βn) (11) n∞=αn/(αn+βn) (11)

h∞=αh/(αh+βh) (12) h∞=αh/(αh+βh) (12)

τm=1/(αm+βm) (13) τm=1/(αm+βm) (13)

τn=1/(αn+βn) (14) τn=1/(αn+βn) (14)

τh=1/(αh+βh) (15) τh=1/(αh+βh) (15)

αm=φ(2.5−0.1(Vm−Vrest))/(e(2.5−0.1(Vm−Vrest))−1) (16) αm=φ(2.5-0.1(Vm-Vrest))/(e(2.5-0.1(Vm-Vrest))-1) (16)

αn=φ(0.1−0.01(Vm−Vrest))/(e(1−0.1(Vm−Vrest))−1) (17) αn=φ(0.1-0.01(Vm-Vrest))/(e(1-0.1(Vm-Vrest))-1) (17)

αh=0.07φe(−(Vm−Vrest))/20 (18) αh=0.07φe(-(Vm-Vrest))/20 (18)

βm=4φe(−(Vm−Vrest))/20 (19) βm=4φe(-(Vm-Vrest))/20 (19)

βn=0.125φe(−(Vm−Vrest))/80 (20) βn=0.125φe(-(Vm-Vrest))/80 (20)

βh=φ/(e(3.0−0.1(Vm−Vrest))+1) (21) βh=φ/(e(3.0-0.1(Vm-Vrest))+1) (21)

Here, φ=3^(T−6.3)/10. The relationship between the transition state and the membrane potential is shown in Figure 2 (Hodgkin and Huxley, 1952, 1989; Börgers, 2017).

FIGURE 2

Figure 2. The relationship between transition state and membrane potential. (A) The evolution of the transition rate [alpha] at the temperature of 6.3°C. (B) The variation of the transition rate [beta] at the temperature of 6.3°C. (C) The evolution of the transition rate [alpha] at the temperature of 18.5°C. (D) The change of the transition rate [beta] at the temperature of 18.5°C.

The HH spiking neuron model is strongly dependent on the temperature, and the early experiments were carried out under the temperatures T = 6.3°C and T = 18.5°C. When the temperature is 6.3°C, the transition rates of the active gates α_n and α_m (Figure 2A), the inactive rate β_h (Figure 2B) increase with the rise of the membrane potential. The inactive transition rate α_h (Figure 2A), the active transition rates β_n and β_m (Figure 2B) decrease with the increase of the membrane potential. When the temperature is increased to 18.5°C, the transition rates α and β show the same experimental phenomena (Figures 2C,D) as above. We compare the transition rates at different temperatures, and the difference is performed in the light blue ellipse. When the temperature is 6.3°C, α_n varies from 0 to 10, α_m alters from 0 to 1, α_h changes from 0.5 to 0 (the enlarged plot in Figure 2A). When the temperature is 18.5°C, α_n varies from 0 to 36, α_m adjusts from 0 to 3.5, α_h changes from 2 to 0 (the enlarged plot in Figure 2C). When the temperature is set to 6.3°C, β_n varies from 37 to 0, β_m adjusts from 0.2 to 0, β_h changes from 0 to 1 (the enlarged plot in Figure 2B). When the temperature is increased to 18.5°C, β_n varies from 140 to 0, β_m adjusts from 0.8 to 0, β_h changes from 0 to 4 (the enlarged plot in Figure 2D). The higher the temperature, the greater the range of conversion rates, the longer time needed to return to the critical value of the transition rate.

When the temperatures are T = 6.3°C and T = 18.5°C, the simulation plots between the steady values of gate variables (m_∞, n_∞, and h_∞) and the membrane potential, the relationship between the time constant (τ_m, τ_n, and τ_h) and the membrane potential, as shown in Figure 3.

FIGURE 3

Figure 3. The relationship between gate-controlled variables, the time constant, and the membrane potential. (A) The evolution of the gate-controlled variables (m, n, and h) at the temperature of 6.3°C. (B) The change of the time constant ([tau]) at the temperature of 6.3°C. (C) The change of the gate-controlled variables (m, n, and h) at the temperature of 18.5°C. (D) The evolution of the time constant ([tau]) at the temperature of 18.5°C.

The steady-state values (m_∞ and n_∞) of activation gate variables (m and n) change from 0 to 1 with the increase of the membrane potential. The steady-state value (h_∞) of the inactivation gate variable (h) decreases with the increase of the membrane potential (Figures 3A,C). The steady-state values are not affected by the change of temperature. When the temperature is 6.3°C, τ_n varies from 5.8 to 1, τ_m adjusts from 0.8 to 0, τ_h changes from 9 to 1. When the temperature is increased to 18.5°C, τ_n varies from 1.5 to 0.25, τ_m adjusts from 0.2 to 0, τ_h changes from 2.25 to 0.25 (Figures 3B,D). The higher temperature, the smaller the range of τ.

3. The Electrical Circuit of the Hodgkin-Huxley Spiking Neuron

The significant electrical properties of a neuron can be precisely replicated by the HH circuit model, as shown in Figure 4A (Hodgkin and Huxley, 1989).

FIGURE 4

Figure 4. The electrical circuit of the axon cell membrane. (A) The HH circuit model. (B) The MHH circuit model.

Here, C is the membrane capacitor. g_Na is the sodium conductance, g_K is the potassium conductance, and g_L is the leaky conductance. V_m is the membrane potential. I_C is the capacitor current, I_Na is the sodium current, I_K is the potassium current, and I_L is the leaky current. I_ext is the external stimulus. E_Na, E_K, and E_L are ion concentration differences of sodium, potassium, and leakage [namely, the equilibrium potentials (Emili et al., 2003) are calculated by the Nernst equation (Hill, 1992)]. The arrow directions of currents are pointing from inside to outside of the membrane. The value of the extracellular potential is set to zero (V_out= 0, namely, the extracellular is grounded) (Hodgkin and Huxley, 1989).

According to Kirchhoff's voltage-current law, the circuit equations are described as:

Vm=Vin−Vout (22) Vm=Vin-Vout (22)

IC=dQ/dt (23) IC=dQ/dt (23)

Q=CVm (24) Q=CVm (24)

Im=INa+IK+IL (25) Im=INa+IK+IL (25)

Iext=IC+INa+IK+IL=IC+Im (26) Iext=IC+INa+IK+IL=IC+Im (26)

In the giant squid axon experiment, the current through the axon membrane is expressed as the current density J(t, x). It represents the amount of the electric current per square centimeter, and its unit is mAcm⁻². Based on the mathematical analysis of the RC equivalent circuit (Figure 4A), the following voltage-current equations are obtained.

C∂Vm(t,x)/∂t=−Jm(t,x)+Jext(t,x)+1/(2rin)∂2Vm(t)/∂x2 (27) C∂Vm(t,x)/∂t=-Jm(t,x)+Jext(t,x)+1/(2rin)∂2Vm(t)/∂x2 (27)

Jm=JNa+JK+JL (28) Jm=JNa+JK+JL (28)

JNa=gNa(Vm−ENa) (29) JNa=gNa(Vm-ENa) (29)

JK=gK(Vm−EK) (30) JK=gK(Vm-EK) (30)

JL=gL(Vm−EL) (31) JL=gL(Vm-EL) (31)

The left side of (27) is the charging or discharging rate per unit area for the capacitor. J_m(t, x) is the total current density that flows through the membrane. J_Na is the current density passing through sodium conductance. J_K is the current density of potassium. V_m is the membrane potential. J_ext(t, x) is the external stimulus. The last term is the charge rate of longitudinal current along the inside membrane surface. It depends only on the time t rather than the location x, so the quadratic partial differential term equals zero, (27) can be rewritten as:

C∂Vm(t,x)/∂t=−Jm(t,x)+Jext(t,x) (32) C∂Vm(t,x)/∂t=-Jm(t,x)+Jext(t,x) (32)

The propagated action potential is performed by (32). The action potential is sensitive to the temperature. The action potential of the cell membrane shows distinct firing behaviors under various temperatures.

When the temperature is 6.3°C, the HH spiking model generates three action potentials in 20 ms, the duration of a spike is 7.65 ms (Figure 5A). When the temperature becomes 15°C, the HH spiking model generates six action potentials in 20 ms, the duration of a spike decreases to 3.35 ms (Figure 5B). When the temperature is increased to 20°C, the HH spiking model generates nine action potentials in 20 ms, the duration of a spike reduces to 1.95 ms (Figure 5C). We increase the temperature to 35°C, and there is no action potential produced after one action potential is generated (Figure 5D). We decrease the temperature to −20°C, and the action potential cannot be obtained (Figure 5E). The temperature affects the time duration of the spike, the generation of action potentials, and the firing frequency of a neuron. It is hard to achieve the action potential when the temperature is too high or low. The increase of temperature has significantly decreased the time duration of the spike and remarkably produced a higher firing frequency.

FIGURE 5

Figure 5. The action potentials under the different temperatures (the external stimulus is 0.08 mAcm−2). (A) The temperature is 6.3°C. (B) The temperature is 15°C. (C) The temperature is 20°C. (D) The temperature is 35°C. (E) The temperature is −20°C.

The external stimuli with various intensities act on the HH spiking model, which performs different action potentials. When the current density is 0.001mAcm⁻², the HH spiking model cannot produce the action potential (Figure 6A). When the current densities are increased to 0.01 and 0.09mAcm⁻², the action potentials are obtained (Figures 6B,C). However, when the current density becomes 0.2mAcm⁻², the HH spiking model generates one action potential. After that, it cannot produce the action potentials (Figure 6D). The external stimulus is related to the generation of the action potential. The larger the external stimulus, the higher the firing frequency. If the external stimulus is too larger or small, the HH spiking model cannot reproduce the action potential.

FIGURE 6

Figure 6. The distinct stimuli are applied to the HH model. (A) The external stimulus is 0.001 mAcm−2. (B) The external stimulus is 0.01 mAcm−2. (C) The external stimulus is 0.09 mAcm−2. (D) The external stimulus is 0.2 mAcm−2.

When the action time of the external stimulus is 1 ms, there is not enough time to show the complete firing process (Figure 7A). Therefore, the action time is increased to 10 ms, and the action potential is generated (Figure 7B). When the action time becomes 20 or 50 ms, the HH spiking model produces more action potentials (Figures 7C,D). Thus, the action time of the external stimulus has a strong influence on the generation of the action potential. The longer the action time, the more action potentials generated. But when the action time is too long or short, the HH spiking model cannot perform the firing process.

FIGURE 7

Figure 7. The firing behaviors under the various action time of the external stimulus. (A) The action time is 1 ms. (B) The action time is 10 ms. (C) The action time is 20 ms. (D) The action time is 50 ms.

4. The Memristive Hodgkin-Huxley (MHH) Spiking Neuron Model

In the HH circuit model, the potassium conductance and the sodium conductance are voltage-gated channels, which can be described by time and membrane potential. The flux-controlled memristor with the nonvolatile property is the function of time and voltage, which can be used in a nonlinear circuit system (Petras, 2010; Corinto and Forti, 2017; Corinto et al., 2018). Based on the HH spiking model, we replace the sodium and potassium conductances with the flux-controlled memristors (Wang et al., 2012), and the memristive Hodgkin-Huxley spiking neuron model is constructed (Figure 4B).

Some of the mathematical expressions in the HH spiking model need to be modified. g_Na and g_K in (4) and (5) are replaced by the memristance and rewritten as:

gMK=1/MKn4 (33) gMK=1/MKn4 (33)

gMNa=1/MNam3h (34) gMNa=1/MNam3h (34)

The conductance values of the sodium and potassium ion channels become the function of time, and the membrane potential will change with the evolution of the memristance.

The flux-controlled memristor is described as (Wang et al., 2012):

M(ϕ(t))=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪20000−3.98×108ϕ(t)+108−−−−−−−−−−−−−−−−−−√ 100ϕ(t)<−0.75ϕ(t)≥−0.75 andϕ(t)<0.25ϕ(t))≥0.25 (35) M(ϕ(t))={20000ϕ(t)<−0.75−3.98×108ϕ(t)+108ϕ(t)≥−0.75 andϕ(t)<0.25100ϕ(t))≥0.25 (35)

Where M_K=M_Na=M is the function of time. The potassium memristance (g_MK) and the sodium memristance (g_MNa) are functions involved with time and membrane potential. When the various external stimuli act on the MHH spiking neuron model, changes in g_MK and g_MNa are performed in Figure 8.

FIGURE 8

Figure 8. Changes of the memristance under the distinct input currents. (A) The variations of the sodium memristance. (B) The variations of the potassium memristance.

The initial values of memductances and reductance g_MK= 0.5×10⁻⁴ Ohm⁻¹cm⁻², g_MNa= 0.5×10⁻⁴ Ohm⁻¹cm⁻², and g_L= 0.3×10⁻³ Ohm⁻¹cm⁻² [0.5×10⁻⁴ is the reciprocal of the maximum value (20,000 Ohmcm⁻²) of a memristor]. The temperature is 6.3°C, C is 1 μF m⁻². E_Na is 50 mV, E_K is −70 mV, and E_L is −50 mV.

When the external stimulus [0.008 mA cm⁻² (g_MNa)] is applied to the MHH spiking model, the sodium memductance (the coral color curve) does not change in the time range from 0 to 1.025 ms (the enlarged plot in Figure 8A). Then, the sodium memductance increases to 0.029 Ohm⁻¹cm⁻² and then decreases to zero. When the MHH spiking model receives the external stimulus [0.08 mAcm⁻² (g_MNa1)], the sodium memductance (the dark red curve) remains the same in the time range from 0ms to 1.38ms (the enlarged plot in Figure 8A). And the maximum value of the sodium memductance is 0.031 Ohm⁻¹cm⁻². Likewise, when the external stimulus [0.8 mAcm⁻² (g_MNa2)] acts on the MHH spiking model, the sodium memductance (the purple curve) does not change in the time range from 0 to 0.97 ms (the enlarged plot in Figure 8A). And the maximum value of the sodium memductance is 0.038 Ohm⁻¹cm⁻².

When the external stimulus [0.04 mAcm⁻² (g_MK)] is injected into the MHH spiking model, the potassium memductance (the coral color curve) does not change from 0 to 1.5 ms (the enlarged plot in Figure 8B). Then, the potassium memductance increases and attains 0.0324 Ohm⁻¹cm⁻². Likewise, the MHH spiking model receives the external stimuli [(0.08 mAcm⁻² (g_MK1) and 0.16 mAcm⁻² (g_MK2)], the potassium memductance (the dark red curve reaches 0.0348 Ohm⁻¹cm⁻² and the purple curve attains 0.0359 Ohm⁻¹cm⁻² (the enlarged plot in Figure 8B) are stable at constant values (Figure 8B).

The sodium memductance and the potassium memductance are associated with the external stimulus. The stronger the external input, the faster the memductance changes, the larger the memductance value. The change curves of sodium and potassium memductance are similar to the theoretical curves (refer to Hodgkin and Huxley, 1989). Therefore, the memristors can mimic the sodium ion channel and the potassium ion channel.

The temperature is selected as 6.3°C, and the external current is 0.08 mAcm⁻². The transition rate parameters (α and β), gate variables (m_∞, n_∞, and h_∞), and the time constant (τ) in the MHH spiking model are shown in Figure 9.

FIGURE 9

Figure 9. The transition rate, gate variables, and time constant of the MHH spiking model. (A) The variation process of the transition [alpha]. (B) The variation process of the transition rate [beta]. (C) The variation process of gate variables. (D) The variation of the time constant [tau].

The transition rates of the active gates (α_n and α_m, Figure 9A), the inactive transiton rate (β_h, Figure 9B) enhance with the increase of the membrane potential. The inactive transition rate (α_h, Figure 9A), the active transition rates (β_n and β_m, Figure 9B) decrease with the rise in the membrane potential. The steady-state values (m_∞ and n_∞) of activation gate variables (m and n) change from 0 to 1 with the increase of the membrane potential. The steady-state value (h_∞) of the inactivation gate variable (h) decreases with the increase of the membrane potential (Figure 9C). The time constant τ_n changes from 4.52 to 0, τ_m adjusts from 0.5 to 0, and τ_h varies from 8.57 to 0 (Figure 9D). The changing processes of the transition rate, gate variables, and the time constant in the MHH spiking model have high similarities with those of the HH spiking model in Figures 2, 3. Therefore, the memristors can be utilized as the sodium ion channel and the potassium ion channel.

When the current density J_m in (28) is replaced by J_M, conductances g_Na and g_K in (29) and (30) are replaced by g_MNa and g_MK, and the current equations are rewritten as:

JM=JMNa+JMK+JL (36) JM=JMNa+JMK+JL (36)

JMNa=gMNa(V−ENa) (37) JMNa=gMNa(V-ENa) (37)

JMK=gMK(V−EK) (38) JMK=gMK(V-EK) (38)

The membrane potential V_m in (32) is replaced by V_M, and the membrane potential of the MHH spiking neuron model is described as:

C∂VM(t,x)/∂t=−JM(t,x)+Jext(t,x) (39) C∂VM(t,x)/∂t=-JM(t,x)+Jext(t,x) (39)

The electrical equivalent circuit of the HH spiking model is based on the voltage-clamp experimental method. When the voltage-clamp values are distinct, the variables perform various variations in the HH and MHH spiking models. Here, the temperature T = 6.3°C. The clamp voltage is denoted by V_clamp, and its value is selected as +20 or +80 mV. The resting potential V_rest = −65 mV. The membrane potential V_m=V_clamp +V_rest.

When the clamp-voltage value is 20 mV, the membrane potential becomes −45 mV. Changes of Na⁺ and K⁺ gate variables in the MHH spiking model (the plots on the left in Figure 10B) are the same as those in the HH spiking model (the plots on the left in Figure 10A). The HH spiking model generates the reverse curves of J_Na and J_m, and their maxima are −0.17 and −0.21 mAcm⁻². The maximum of the forward curve J_K is 0.14 mAcm⁻², and the forward curve J_L reaches 0.009 mAcm⁻². The peak values of g_K and g_Na are 4.54 and 2.25 mOhm⁻¹cm⁻² (the plots on the right in Figure 10A). The MHH spiking model produces the reverse curves of J_MNa and J_M, and their maxima are −0.18 and −0.16 mAcm⁻². The forward curves of J_MK and J_L attain their maxima 0.04 and 0.009 mAcm⁻². The maxima of g_MK and g_MNa are 1.26 and 1.88 mOhm⁻¹cm⁻² (the plots on the right in Figure 10B).

FIGURE 10

Figure 10. The distinct clamp voltages are applied to the HH spiking model and the MHH spiking model. (A) The HH model with Vclamp = 20 mV. (B) The MHH model with Vclamp = 20 mV. (C) The HH model with Vclamp = 80 mV. (D) The MHH model with Vclamp = 80 mV.

The variable values of the HH spiking model are more significant than those of the MHH spiking model (because the memristance is large, its initial value is 10,000 Ohmcm⁻²). When the clamp-voltage value is 20 mV, both spiking models cannot generate the action potential.

A transient increase of sodium ions in the cell leads to the depolarization of the action potential. The waveforms of the two models change in the same way when the clamp voltage is 80 mV (the membrane potential is 15 mV). We take the MHH model as an example and make a vertical comparison (Figures 10B,D). With the increase of clamp voltage, the current densities of sodium and potassium increase significantly. The value of gate variable n changes from 0.5 to 1, and the value of gate variable m varies from 0.4 to 1. The potassium memductance changes from 1.26 to 8 mOhm⁻¹cm⁻², and the sodium memductance changes from 1.88 to 30 mOhm⁻¹cm⁻².

When the clamp-voltage value is 80 mV, the HH and MHH spiking models can produce the action potential. The gate variables n and m change with the identical waveforms. The current densities, the potassium conductance, and the sodium conductance are different. The maxima of J_MNa, J_MK, J_L, and J_M are −1.059, 0.74, 0.0297, and −0.97 mAcm⁻² (the right-upper plot in Figure 10D), which are larger than those of the HH spiking model (Figure 10C). The variation ranges of potassium conductance and sodium conductance for the MHH spiking model are [0 8], [0 30] less than those [0 29], [0 37] in the HH spiking model. The higher the voltage-clamp value, the larger the variable values, the smaller the conductance variation range.

5. The Comparison Between Two Models Under the Different Stimuli

5.1. The Individual Current Pulse Stimulus

The forward stimulus J_ext = 0.1mAcm⁻² (the pulse width is 0.1 ms) is applied to the HH spiking model and the MHH spiking model, the temperature is selected as 18.5°C, and the response time of the model is 5 ms. The initial value of the membrane potential is the resting potential, V_rest = −65 mV.

Here, J_ext is the external stimulus, J_Na (J_MNa) is the sodium current (the coral color curve), J_K (J_MK) is the potassium current (the blue curve), J_L (J_ML) is the leaky current (the green curve), and J_m (J_M) is the total current (the purple curve) flowing through the cell membrane in the HH (MHH) spiking model. V (V_M) is the action potential generated by the HH (MHH) spiking model. g_Na (g_MNa) is the sodium conductance (the sodium memductance), and g_K (g_MK) is the potassium conductance (the potassium memductance) in the HH (MHH) spiking model.

The HH and MHH spiking models receive the external stimuli and produce the corresponding current densities of the ion channels. The sodium current is negative because the sodium ions move from the outside to the inside of the cell. In contrast, the potassium current is positive because the potassium ions flow from intracellular to extracellular. The potassium and total current densities (the peak values: J_K = 0.82 mAcm⁻², J_m= −0.51 mAcm⁻²) generated by the HH spiking model are larger than those (the peak values: J_MK = 0.4 mAcm⁻², J_M = −0.53 mAcm⁻²) in the MHH spiking model. The sodium and leaky current densities (the peak values: J_Na= −0.7 mAcm⁻², J_L= 0.024 mAcm⁻²) generated by the HH spiking model are smaller than those (the peak values: J_MNa= −0.6 mAcm⁻², J_L= 0.026 mAcm⁻²) in the MHH spiking model. The sodium current of the MHH model has a smooth perturbation at around t = 1.072 s, and the sodium current of the HH model has an obvious perturbation at around t = 1.279 s. The perturbation is caused by the rapid variation of potassium conductance (potassium memductance). The curves formed by the MHH model (the left side plot in Figure 11B) are smoother than those in the HH model (the left side plot in Figure 11A) because the memristor has a unique time-varying property.

FIGURE 11

Figure 11. The single pulse and the reverse pules are applied to two models. (A) The single current pulse is injected into the HH spiking model. (B) The single current pulse is injected into the MHH spiking model. (C) The reverse current pulse is injected into the HH spiking model. (D) The reverse current pulse is injected into the MHH spiking model.

The HH spiking model and the MHH spiking model can perform the action potential. The membrane potential peak value (V_M = 38.33 mV at 1.188 ms) of the MHH model (the middle plot in Figure 11B) is stronger than that (V_m = 28.31 mV at 1.366 ms) of the HH model (the central plot in Figure 11A). Meanwhile, the MHH spiking model takes a short time to produce the action potential. After generating the action potential, both models return to the equilibrium state (the resting state, V_rest = −65 mV).

The HH spiking model takes 1.354 ms to reach the maximum value of g_Na (23.53 mOhm⁻¹cm⁻²) and needs 1.715 ms to get the peak value of g_K (12.45 mOhm⁻¹cm⁻²; the right side plot in Figure 11A). Therefore, the MHH spiking model takes 1.134 ms to attain the maximum value of g_MNa (20.81 mOhm⁻¹cm⁻²) and needs 1.673 ms to reach the peak value of g_MK (5.196 mOhm⁻¹cm⁻²) (the right side plot in Figure 11B). The rise in sodium conductance (sodium memductance) is faster than potassium conductance (potassium memductance). The MHH spiking model utilizes less time than the HH model to activate the change of the memductance; however, the obtained memductane is small. Because the variation in the memductance is slight in a short time (5 ms), it maintains a large memristance.

5.2. The Reverse Single Current Pulse Stimulus

The reverse stimulus (J_ext = −0.1 mAcm⁻², the pulse width is 0.1 ms) acts on the HH spiking model and the MHH spiking model, the temperature is 18.5°C, and the response time of the model is 5 ms.

There are not enough ions to move from intracellular (extracellular) to extracellular (intracellular); therefore, the sodium current and the potassium current cannot be produced (the left-side plots in Figures 11C,D). The significant variation of the conductance causes the generation of potassium and sodium currents. The sodium conductance (sodium memductance) is close to zero (the right-side plots in Figures 11C,D). The potassium conductance (potassium memductance) decreases from 0.37 mOhm⁻¹cm⁻² (0.36 mOhm⁻¹cm⁻²) to 0.17 mOhm⁻¹cm⁻² (0.14 mOhm⁻¹cm⁻²) and then increases to 0.35 mOhm⁻¹cm⁻² (0.26 mOhm⁻¹cm⁻²). The HH and MHH spiking models are unable to generate the action potential, and the membrane potentials become hyperpolarization before returning to their resting states (the middle plots in Figures 11C,D.

5.3. The Three External Stimuli With Different Intensity

The external stimuli J_ext1 = 0.5 mAcm⁻², J_ext2 =1 mAcm⁻², and J_ext3 =2 mAcm⁻² are injected into the HH spiking model and the MHH spiking model, the temperature is 18.5°C, the response time is 5 ms.

When the small external stimulus (J_ext1 = 0.5 mA.cm⁻²) is applied to the HH spiking model, the action potential cannot be produced. The membrane potential has a slight rise (V_m = −60 mv) and then returns to the resting potential (−65 mv) at 3 ms (the second plot in Figure 12A). The current density is zero (the first plot in Figure 12A). There is only a slight change in the conductance, which can be ignored (the third plot in Figure 12A). However, when the MHH spiking receives the stimulus J_ext1 = 0.5 mA.cm⁻², the action potential is obtained (the second plot in Figure 12B). The changes in current densities and the memductance are noticeable. When the external stimuli increase to J_ext2 = 1 mA.cm⁻² and J_ext1 = 2 mA.cm⁻², the values in current density, membrane potential, and conductances strengthen gradually (Figure 12).

FIGURE 12

Figure 12. The three different current pulses are applied to the HH and MHH spiking models. (A) The HH spiking model. (B) The MHH spiking model.

The larger the external stimulus, the faster the action potential is produced, the higher the peak value is generated, the more significant change in conductances, and the greater the current density. The smaller the external stimulus, the longer time it takes to produce the action potential. The peak value of membrane potential in the MHH model (the middle plot in Figure 12B) is greater than that of the HH model (the middle plot in Figure 12A). The maximum values of current densities and conductances in the MHH spiking model (the first and third plots in Figure 12B) are lower than those in the HH spiking model (the first and third plots in Figure 12A).

5.4. A Series of Pulse Stimuli

When a series of pulses (J_ext(n)= 1mAcm⁻², n = 1,2,......,18, the temperature is 18.5°C.) act on the HH and MHH spiking models, the action potentials are achieved. However, not every single pulse can cause the generation of the action potential (the first plots in Figures 13A,B). Only when the action potential generated by the previous pulse has enough time to return to its resting state, another action potential will be generated. The MHH spiking model [six action potentials (the second plot in Figure 13B)] generates more action potentials than the HH spiking model (five action potentials (the second plot in Figure 13A)). Meanwhile, the action potential performs two oscillation behaviors in the MHH spiking model (inside the blue ellipse in Figure 13B), and the action potential shows three oscillation behaviors in the HH spiking model (inside the blue ellipse in Figure 13A). The memductances in the MHH model (the third plot in Figure 13B) are smaller than those in the HH model (the third plot in Figure 13A), which causes the current density produced by the MHH model (the first plot in Figure 13B) to be lower than the HH model (the first plot in Figure 13A).

FIGURE 13

Figure 13. The distinct action time of the external stimulus is set for the two models. (A) The HH spiking model with 20 ms action time. (B) The MHH spiking model with 20 ms action time. (C) The HH spiking model with 100 ms action time. (D) The MHH spiking model with 100 ms action time. (E) The HH spiking model with 200 ms action time. (F) The MHH spiking model with 200 ms action time.

The action time of the external stimulus is extended to 100 ms, and two models can produce more action potentials than Figures 13A,B. The MHH spiking model generates more action potentials (the middle plot in Figure 13C) than the HH spiking model (the middle plot in Figure 13D).

The action time is increased to 200 ms, the doublet currents (Shigaki et al., 2020) are generated in the MHH spiking model, one is large, the other is small (the enlarged plot inside the left ellipse in Figure 13F). Meanwhile, the action potential is produced before the current pulse comes in the MHH model because the memristor has an initial charge even though it is very small (the enlarged plot inside the right ellipse in Figure 13F). The current intensity, the voltage peak value, and conductances in the HH spiking model (Figure 13E) are larger than the simulation results in Figure 13F.

With the increasing of time length, the conductance (or memductance) and the current density of sodium and potassium increase dramatically. The more time we give, the more action potentials are generated, the larger the peak values of current densities, conductances (or memductances), and action potentials. However, the action time length should not be too long; otherwise, the function of neurons cannot be replicated effectively (Chen et al., 2019).

5.5. The Sinusoidal Current Stimulus

The sinusoidal stimulus [J_ext = J_extm× sin(2 t/T_in), J_extm=0.01 mA.cm⁻²] is a positive-negative periodic signal with a single-frequency component. T_in is the time period of input signals, and the temperature is 18.5°C.

When T_in = 0.01 ms and T_in = 1 ms, the sinusoidal stimuli are applied to the HH spiking model. The action potential cannot be obtained because there is not enough time for the neuron to depolarize. But the MHH model generates action potentials under the same conditions. The frequency of the sinusoidal stimulus affects the generation of the action potential. When the frequency is low, there is sufficient time to depolarize, and the action potential occurs (Figure 14). When T_in = 5 ms, the HH and MHH spiking models produce the action potentials, their spiking patterns belong to tonic spikes in pyramidal neurons. When T_in = 20 ms, the MHH model generates the repetitive bursts with doublet spikes, and the HH model performs the tonic spiking. When the value of T_in is increased to 60 ms, the action potential cannot be produced in the HH spiking model but can be obtained in the MHH model. The frequency range of the sinusoidal stimulus in the MHH spiking model is wider than that of the HH spiking model. The various spiking patterns can be obtained by appropriately adjusting the frequency of the sinusoidal signal.

FIGURE 14

Figure 14. The time periods of the sinusoidal signal. Tin = 0.01 ms, Tin = 1 ms, Tin = 5 ms, Tin = 20 ms, Tin = 60 ms.

6. Conclusion

The biological neuron is expressed adequately by the classic HH spiking model. It is sensitive to the temperature, the strength of the external stimulus, and the action time of the stimulus. The MHH spiking model successfully simulates the generation of the action potential in a neuron. When the different external stimuli are applied to the HH and MHH spiking models, the action potential is produced, and various spiking patterns are achieved. The MHH spiking model has advantages in generating the action potential through the comparison with the HH spiking model. The waveforms with smaller perturbations formed by the MHH spiking model are smooth. The higher frequency of the external stimulus, the more action potentials generated. The response speed of the MHH spiking model is faster than that of the HH spiking model. The various spiking behaviors are obtained by adjusting the signal frequency in the MHH spiking model. And meanwhile, the combination between neuron models and a memristor provides the possibility to scale down the neuron circuit and gives a novel way to replicate the functions of the biological neuron.

Data Availability Statement

The original contributions generated for the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

XF built models and simulations, carried out the experimental analysis, and prepared the manuscript in this work. SD and LW supervised the content of the article and the results of the simulations. All authors contributed to the article and approved the submitted version.

Funding

Project supported by the National Key R and D Program of China (Grant No. 2018YFB1306600), the National Natural Science Foundation of China (Grant Nos. 62076207, 62076208, U20A20227), the Fundamental Science and Advanced Technology Research Foundation of Chongqing, China (Grant No. cstc2017jcyjBX0050).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Beck, M. E., Shylendra, A., Sangwan, V. K., Guo, S., Gaviria Rojas, W. A., Yoo, H., et al. (2020). Spiking neurons from tunable gaussian heterojunction transistors. Nat. Commun. 11:1565. doi: 10.1038/s41467-020-15378-7

Berdan, R., Khat, A., Papavassiliou, C., and Prodromakis, T. (2014). “Qualitative spice modeling accounting for volatile dynamics of TIO2 memristors,” in IEEE International Symposium on Circuits and Systems (ISCAS) (Melbourne, VIC), 2033–2036.

Berdan, R., Vasilaki, E., Khiat, A., Indiveri, G., Serb, A., and Prodromakis, T. (2016). Emulating short-term synaptic dynamics with memristive devices. Sci. Rep. Nat. 6:18639. doi: 10.1038/srep18639

Börgers (2017). An introduction to modeling neuronal dynamics. Texts Appl. Math. 66, 1–42. doi: 10.1007/978-3-319-51171-9

Buhry, L., Grassia, F., Giremus, A., Grivel, E., Renaud, S., and Saghi, S. (2011). Automated parameter estimation of the Hodgkin-Huxley model using the differential evolution algorithm: application to neuromimetic analog integrated circuits. Neural Comput. 23, 2599–2625. doi: 10.1162/NECO_a_00170

Chen, B. S., and Li, C. W. (2010). On the noise-enhancing ability of stochastic hodgkin-huxley neuron systems. Neural Comput. 22, 1737–1763. doi: 10.1162/neco.2010.07-09-1057

Chen, S., Zou, Y., and Zhang, X. (2019). An efficient method for hopf bifurcation control in fractional-order neuron model. IEEE Access 7, 77490–77498. doi: 10.1109/ACCESS.2019.2920007

Chua, L. (1971). Memristor–the missing circuit element. IEEE Trans. Circuit Theory 18, 507–519. doi: 10.1109/TCT.1971.1083337

Corinto, F., and Forti, M. (2017). Memristor circuits: bifurcations without parameters. IEEE Trans. Circuits Syst. I 64, 1540–1551. doi: 10.1109/TCSI.2016.2642112

Corinto, F., Gilli, M., and Forti, M. (2018). Flux-charge description of circuits with non-volatile switching memristor devices. IEEE Trans. Circuits Syst. II 65, 642–646. doi: 10.1109/TCSII.2018.2825447

Dröscher, C. (1998). Camillo golgi and the discovery of the Golgi apparatus. Histochem. Cell Biol. 109, 425–430. doi: 10.1007/s004180050245

Emili, Schiavoni, Francavilla, Roselli, and Sorrentino (2003). Computation of electromagnetic field inside a tissue at mobile communications frequencies. IEEE Trans. Microw. Theory Techn. 51, 178–186. doi: 10.1109/TMTT.2002.806899

Ghavami, S., Rahmati, V., Lahouti, F., and Schwabe, L. (2018). Neuronal synchronization can control the energy efficiency of inter-spike interval coding. IEEE Trans. Mol. Biol. Multi-Scale Commun. 4, 221–236. doi: 10.1109/TMBMC.2019.2937291

Häusser, M. (2000). The Hodgkin-Huxley theory of the action potential. Nat. Neurosci. 3, 1165–1165. doi: 10.1038/81426

Hill, B. (1992). Ionic Channels of Excitable Membranes. Sunderlan, MA: Sinauer Assoociates.

Hodgkin, A. L., and Huxley, A. F. (1952). The components of membrane conductance in the giant axon of loligo. J. Physiol. 116, 473–496. doi: 10.1113/jphysiol.1952.sp004718

Hodgkin, A. L., and Huxley, A. F. (1989). A quantitative description of membrane current and its application to conduction and excitation in nerve. Bull. Math. Biol. 52, 25–71. doi: 10.1016/S0092-8240(05)80004-7

James, A. P. (2019). A hybrid memristor-CMOS chip for AI. Nat. Electron. 2, 268–269. doi: 10.1038/s41928-019-0274-6

Le, Z., Sangho, S., and Sung-Mo, K. (2015). “Memristor-based synapses and neurons for neuromorphic computing,” in IEEE International Symposium on Circuits and Systems (ISCAS) (Lisbon), 1150–1153.

Li, C., Wang, Z., Rao, M., Belkin, D., Song, W., Jiang, H., et al. (2019). Long short-term memory networks in memristor crossbar arrays. Nat. Mach. Intell. 1, 49–57. doi: 10.1038/s42256-018-0001-4

Maheshwar, P. S., Kim, H., and Chua, L. O. (2014). Brains are made of memristors. Circuits Syst. Mag. IEEE 14, 12–36. doi: 10.1109/MCAS.2013.2296414

Mandal, S. and Saha, A. (2016). “Memristors act as synapses in neuromorphic architectures,” in IEEE International Conference on Communication and Electronics Systems (ICCES) (Coimbatore), 1–4. doi: 10.1109/CESYS.2016.7889862

Mokhtar, S. M. A., Abdullah, W. F. H., Kadiran, K. A., Rifin, R., and Omar, M. (2017). “Programmable delay element using memristor and case study in delay lock loop,” in Control & System Graduate Research Colloquium (Shah Alam), 17–21. doi: 10.1109/ICSGRC.2017.8070560

Petras, I. (2010). Fractional-order memristor-based Chua's circuit. IEEE Trans. Circuits Syst. II 57, 975–979. doi: 10.1109/TCSII.2010.2083150

Pu, S., and Thomas, P. J. (2020). Fast and accurate Langevin simulations of stochastic Hodgkin-Huxley dynamics. Neural Comput. 32, 1775–1835. doi: 10.1162/neco_a_01312

Raviola, E., and Mazzarello, P. (2011). The diffuse nervous network of camillo golgi: facts and fiction. Brain Res. Rev. 66, 425–430. doi: 10.1016/j.brainresrev.2010.09.005

Saïgai, S., Bornat, Y., Tomas, J., Le Masson, G., and Renaud, S. (2011). A library of analog operators based on the hodgkin-huxley formalism for the design of tunable, real-time, silicon neurons. IEEE Trans. Biomed. Circuits Syst. 5, 3–19. doi: 10.1109/TBCAS.2010.2078816

Sayyaparaju, S., Amer, S., and Rose, G. S. (2018). “A bi-memristor synapse with spike-timing-dependent plasticity for on-chip learning in memristive neuromorphic systems,” in IEEE19th International Symposium on Quality Electronic Design (ISQED) (Santa Clara, CA), 69–74. doi: 10.1109/ISQED.2018.8357267

Shigaki, S., Ohashi, H., Sakurai, T., Shimizu, M., Hosoda, K., and Kurabayashi, D. (2020). Real-time odor discrimination using single antenna of insect. IEEE Sensors Lett. 4, 1–4. doi: 10.1109/LSENS.2020.3024606

Strukov, D. B., Snider, G. S., Stewart, D. R., and Williams, R. S. (2008). The missing memristor found. Nature 453, 80–83. doi: 10.1038/nature06932

Truong, S. C., Pham, K. V., and Yang, W. (2016). “Memristor circuits and systems for future computing and bio-inspired information processing,” in Biomedical Circuits & Systems Conference (Shanghai), 456–459. doi: 10.1109/BioCAS.2016.7833830

Tuckwell, H. C., and Ditlevsen, S. (2016). The space-clamped Hodgkin-Huxley system with random synaptic input: Inhibition of spiking by weak noise and analysis with moment equations. Neural Comput. 28, 2129–2161. doi: 10.1162/NECO_a_00881

Wang, L., Drakakis, E., and Duan, S. (2012). Memristor model and its application for chaos generation. Int. J. Bifurc. Chaos 22:1250205. doi: 10.1142/S0218127412502057

Xia, Q., and Yang, J. J. (2019). Memristive crossbar arrays for brain-inspired computing. Nat. Mater. 18, 309–323. doi: 10.1038/s41563-019-0291-x

Yang, J. J., Strukov, D. B., and Stewart, D. R. (2013). Memristive devices for computing. Nat. Nanotech. 8, 13–24. doi: 10.1038/nnano.2012.240