Biological neuron model

One of the earliest models of a neuron was first investigated in 1907 by Louis Lapicque.[2] A neuron is represented in time by

${\displaystyle I(t)=C_{\mathrm {m} }{\frac {dV_{\mathrm {m} }(t)}{dt}}}$

which is just the time derivative of the law of capacitance, Q = CV. When an input current is applied, the membrane voltage increases with time until it reaches a constant threshold V_th, at which point a delta function spike occurs and the voltage is reset to its resting potential, after which the model continues to run. The firing frequency of the model thus increases linearly without bound as input current increases.

The model can be made more accurate by introducing a refractory period t_ref that limits the firing frequency of a neuron by preventing it from firing during that period. Through some calculus involving a Fourier transform, the firing frequency as a function of a constant input current thus looks like

${\displaystyle \,\!f(I)={\frac {I}{C_{\mathrm {m} }V_{\mathrm {th} }+t_{\mathrm {ref} }I}}}$.

A remaining shortcoming of this model is that it implements no time-dependent memory. If the model receives a below-threshold signal at some time, it will retain that voltage boost forever until it fires again. This characteristic is clearly not in line with observed neuronal behavior.

The Hodgkin–Huxley model (H&H model)[3][4][5][6] is a model of the relationship between the flow of ionic currents across the neuronal cell membrane and the membrane voltage of the cell.[3][4][5][6] It consists of a set of nonlinear differential equations describing the behaviour of ion channels that permeate the cell membrane of the squid giant axon. Hodgkin and Huxley were awarded the 1963 Nobel Prize in Physiology or Medicine for this work.

We note as before our voltage-current relationship, this time generalized to include multiple voltage-dependent currents:

${\displaystyle C_{\mathrm {m} }{\frac {dV(t)}{dt}}=-\sum _{i}I_{i}(t,V)}$.

Each current is given by Ohm's law as

${\displaystyle I(t,V)=g(t,V)\cdot (V-V_{\mathrm {eq} })}$

where g(t,V) is the conductance, or inverse resistance, which can be expanded in terms of its constant average ḡ and the activation and inactivation fractions m and h, respectively, that determine how many ions can flow through available membrane channels. This expansion is given by

${\displaystyle g(t,V)={\bar {g}}\cdot m(t,V)^{p}\cdot h(t,V)^{q}}$

and our fractions follow the first-order kinetics

${\displaystyle {\frac {dm(t,V)}{dt}}={\frac {m_{\infty }(V)-m(t,V)}{\tau _{\mathrm {m} }(V)}}=\alpha _{\mathrm {m} }(V)\cdot (1-m)-\beta _{\mathrm {m} }(V)\cdot m}$

with similar dynamics for h, where we can use either τ and m_∞ or α and β to define our gate fractions.

The Hodgkin–Huxley model may be expanded to include additional ionic currents. Typically, these include inward Ca²⁺ and Na⁺ input currents, as well as several varieties of K⁺ outward currents, including a "leak" current.

The end result can be at the small end 20 parameters which one must estimate or measure for an accurate model, and for complex systems of neurons not easily tractable by computer. Careful simplifications of the Hodgkin–Huxley model are therefore needed.

In the leaky integrate-and-fire model, the memory problem is solved by adding a "leak" term to the membrane potential, reflecting the diffusion of ions that occurs through the membrane when some equilibrium is not reached in the cell. The model looks like

${\displaystyle I(t)-{\frac {V_{\mathrm {m} }(t)}{R_{\mathrm {m} }}}=C_{\mathrm {m} }{\frac {dV_{\mathrm {m} }(t)}{dt}}}$

where R_m is the membrane resistance, as we find it is not a perfect insulator as assumed previously. This forces the input current to exceed some threshold I_th = V_th / R_m in order to cause the cell to fire, else it will simply leak out any change in potential. The firing frequency thus looks like

${\displaystyle f(I)={\begin{cases}0,&I\leq I_{\mathrm {th} }\\{[}t_{\mathrm {ref} }-R_{\mathrm {m} }C_{\mathrm {m} }\log(1-{\tfrac {V_{\mathrm {th} }}{IR_{\mathrm {m} }}}){]}^{-1},&I>I_{\mathrm {th} }\end{cases}}}$

which converges for large input currents to the previous leak-free model with refractory period.[7] The model can also be extended to inhibitory neurons.[8]

Recent advances in computational and theoretical fractional calculus lead to a new form of model, called Fractional-order leaky integrate-and-fire developed by Teka et al. The great advantage of this model is that it can capture and integrate all the past activities and can reproduce the time dependent spiking adaptations observed on pyramidal neurons. The model has the following form - more details can be found in.[9]

${\displaystyle I(t)-{\frac {V_{\mathrm {m} }(t)}{R_{\mathrm {m} }}}=C_{\mathrm {m} }{\frac {d^{\alpha }V_{\mathrm {m} }(t)}{d^{\alpha }t}}}$

3D visualization of the Galves–Löcherbach model for biological neural nets. This visualization is set for 4,000 neurons (4 layers with one population of inhibitory neurons and one population of excitatory neurons each) at 180 intervals of time.

The Galves–Löcherbach model is a specific development of the leaky integrate-and-fire model. It is inherently stochastic. It was developed by mathematicians Antonio Galves and Eva Löcherbach.[10] Given the model specifications, the probability that a given neuron ${\displaystyle i}$ spikes in a time period ${\displaystyle t}$ may be described by

${\displaystyle \mathop {\mathrm {Prob} } (X_{t}(i)=1|{\mathcal {F}}_{t-1})=\phi _{i}{\Biggl (}\sum _{j\in I}W_{j\rightarrow i}\sum _{s=L_{t}^{i}}^{t-1}g_{j}(t-s)X_{s}(j),~~~t-L_{t}^{i}{\Biggl )},}$

where ${\displaystyle W_{j\rightarrow i}}$ is a synaptic weight, describing the influence of neuron ${\displaystyle j}$ on neuron ${\displaystyle i}$, ${\displaystyle g_{j}}$ expresses the leak, and ${\displaystyle L_{t}^{i}}$ provides the spiking history of neuron ${\displaystyle i}$ before ${\displaystyle t}$, according to

${\displaystyle L_{t}^{i}={\text{sup}}\{s<t:X_{s}(i)=1\}.}$

In the exponential integrate-and-fire model, spike generation is exponential, following the equation:

${\displaystyle {\frac {dV}{dt}}={\frac {1}{\tau _{m}}}[E_{m}-V+\Delta _{T}\exp \left({\frac {V-V_{T}}{\Delta _{T}}}\right)]}$.

where ${\displaystyle V}$ is the membrane potential, ${\displaystyle V_{T}}$ is the membrane potential threshold, ${\displaystyle \tau _{m}}$ is the membrane time constant, ${\displaystyle E_{m}}$is the resting potential, and ${\displaystyle \Delta _{T}}$ is the sharpness of action potential initiation, usually around 1 mV for cortical pyramidal neurons. Once the membrane potential crosses ${\displaystyle V_{T}}$, it diverges to infinity in finite time.[11][12][13]

Sweeping simplifications to Hodgkin–Huxley were introduced by FitzHugh and Nagumo in 1961 and 1962. Seeking to describe "regenerative self-excitation" by a nonlinear positive-feedback membrane voltage and recovery by a linear negative-feedback gate voltage, they developed the model described by

${\displaystyle {\begin{array}{rcl}{\dfrac {dV}{dt}}&=&V-V^{3}/3-w+I_{\mathrm {ext} }\\\tau {\dfrac {dw}{dt}}&=&V-a-bw\end{array}}}$

where we again have a membrane-like voltage and input current with a slower general gate voltage w and experimentally-determined parameters a = -0.7, b = 0.8, τ = 1/0.08. Although not clearly derivable from biology, the model allows for a simplified, immediately available dynamic, without being a trivial simplification.[14]

In 1981 Morris and Lecar combined the Hodgkin–Huxley and FitzHugh–Nagumo models into a voltage-gated calcium channel model with a delayed-rectifier potassium channel, represented by

${\displaystyle {\begin{array}{rcl}C{\dfrac {dV}{dt}}&=&-I_{\mathrm {ion} }(V,w)+I\\{\dfrac {dw}{dt}}&=&\phi \cdot {\dfrac {w_{\infty }-w}{\tau _{w}}}\end{array}}}$

where ${\displaystyle I_{\mathrm {ion} }(V,w)={\bar {g}}_{\mathrm {Ca} }m_{\infty }\cdot (V-V_{\mathrm {Ca} })+{\bar {g}}_{\mathrm {K} }w\cdot (V-V_{\mathrm {K} })+{\bar {g}}_{\mathrm {L} }\cdot (V-V_{\mathrm {L} })}$.[7]

Building upon the FitzHugh–Nagumo model, Hindmarsh and Rose proposed in 1984 a model of neuronal activity described by three coupled first order differential equations:

${\displaystyle {\begin{array}{rcl}{\dfrac {dx}{dt}}&=&y+3x^{2}-x^{3}-z+I\\{\dfrac {dy}{dt}}&=&1-5x^{2}-y\\{\dfrac {dz}{dt}}&=&r\cdot (4(x+{\tfrac {8}{5}})-z)\end{array}}}$

with r² = x² + y² + z², and r ≈ 10⁻² so that the z variable only changes very slowly. This extra mathematical complexity allows a great variety of dynamic behaviors for the membrane potential, described by the x variable of the model, which include chaotic dynamics. This makes the Hindmarsh–Rose neuron model very useful, because being still simple, allows a good qualitative description of the many different patterns of the action potential observed in experiments.

Cable theory describes the dendritic arbor as a cylindrical structure undergoing a regular pattern of bifurcation, like branches in a tree. For a single cylinder or an entire tree, the input conductance at the base (where the tree meets the cell body, or any such boundary) is defined as

${\displaystyle G_{in}={\frac {G_{\infty }\tanh(L)+G_{L}}{1+(G_{L}/G_{\infty })\tanh(L)}}}$,

where L is the electrotonic length of the cylinder which depends on its length, diameter, and resistance. A simple recursive algorithm scales linearly with the number of branches and can be used to calculate the effective conductance of the tree. This is given by

${\displaystyle \,\!G_{D}=G_{m}A_{D}\tanh(L_{D})/L_{D}}$

where A_D = πld is the total surface area of the tree of total length l, and L_D is its total electrotonic length. For an entire neuron in which the cell body conductance is G_S and the membrane conductance per unit area is G_md = G_m / A, we find the total neuron conductance G_N for n dendrite trees by adding up all tree and soma conductances, given by

${\displaystyle G_{N}=G_{S}+\sum _{j=1}^{n}A_{D_{j}}F_{dga_{j}}}$,

where we can find the general correction factor F_dga experimentally by noting G_D = G_mdA_DF_dga.

The cable model makes a number of simplifications to give closed analytic results, namely that the dendritic arbor must branch in diminishing pairs in a fixed pattern. A compartmental model allows for any desired tree topology with arbitrary branches and lengths, but makes simplifications in the interactions between branches to compensate. Thus, the two models give complementary results, neither of which is necessarily more accurate.

Each individual piece, or compartment, of a dendrite is modeled by a straight cylinder of arbitrary length l and diameter d which connects with fixed resistance to any number of branching cylinders. We define the conductance ratio of the ith cylinder as B_i = G_i / G_∞, where ${\displaystyle G_{\infty }={\tfrac {\pi d^{3/2}}{2{\sqrt {R_{i}R_{m}}}}}}$ and R_i is the resistance between the current compartment and the next. We obtain a series of equations for conductance ratios in and out of a compartment by making corrections to the normal dynamic B_out,i = B_in,i+1, as

• ${\displaystyle B_{\mathrm {out} ,i}={\frac {B_{\mathrm {in} ,i+1}(d_{i+1}/d_{i})^{3/2}}{\sqrt {R_{\mathrm {m} ,i+1}/R_{\mathrm {m} ,i}}}}}$
• ${\displaystyle B_{\mathrm {in} ,i}={\frac {B_{\mathrm {out} ,i}+\tanh X_{i}}{1+B_{\mathrm {out} ,i}\tanh X_{i}}}}$
• ${\displaystyle B_{\mathrm {out,par} }={\frac {B_{\mathrm {in,dau1} }(d_{\mathrm {dau1} }/d_{\mathrm {par} })^{3/2}}{\sqrt {R_{\mathrm {m,dau1} }/R_{\mathrm {m,par} }}}}+{\frac {B_{\mathrm {in,dau2} }(d_{\mathrm {dau2} }/d_{\mathrm {par} })^{3/2}}{\sqrt {R_{\mathrm {m,dau2} }/R_{\mathrm {m,par} }}}}+\ldots }$

where the last equation deals with parents and daughters at branches, and ${\displaystyle X_{i}={\tfrac {l_{i}{\sqrt {4R_{i}}}}{\sqrt {d_{i}R_{m}}}}}$. We can iterate these equations through the tree until we get the point where the dendrites connect to the cell body (soma), where the conductance ratio is B_in,stem. Then our total neuron conductance is given by

${\displaystyle G_{N}={\frac {A_{\mathrm {soma} }}{R_{\mathrm {m,soma} }}}+\sum _{j}B_{\mathrm {in,stem} ,j}G_{\infty ,j}}$.

An example of a compartmental model of a neuron, with an algorithm to reduce the number of compartments (increase the computational speed) and yet retain the salient electrical characteristics, can be found in.[15]

Siebert[16][17] modeled the neuron spike firing pattern using a non-homogeneous Poisson process model, following experiments involving the auditory system.[16][17] According to Siebert, the probability of a spiking event at the time interval ${\displaystyle [t,t+\Delta _{t}]}$ is proportional to a non negative function ${\displaystyle g[s(t)]}$, where ${\displaystyle s(t)}$ is the raw stimulus.:

${\displaystyle P_{spike}(t\in [t',t'+\Delta _{t}])=\Delta _{t}\cdot g[s(t)]}$

Siebert considered several functions as ${\displaystyle g[s(t)]}$, including ${\displaystyle g[s(t)]\propto s^{2}(t)}$ for low stimulus intensities.

The main advantage of Siebert's model is its simplicity. The shortcomings of the model is its inability to reflect properly the following phenomena:

• The edge emphasizing property of the neuron in response to a stimulus pulse.
• The saturation of the firing rate.
• The values of inter-spike-interval-histogram at short intervals values (close to zero).

These shortcoming are addressed by the two state Markov Model.[18][19][20]

The spiking neuron model by Nossenson & Messer[18][19][20] produces the probability of the neuron to fire a spike as a function of either an external or pharmacological stimulus.[18][19][20] The model consists of a cascade of a receptor layer model and a spiking neuron model, as shown in Fig 4. The connection between the external stimulus to the spiking probability is made in two steps: First, a receptor cell model translates the raw external stimulus to neurotransmitter concentration, then, a spiking neuron model connects between neurotransmitter concentration to the firing rate (spiking probability). Thus, the spiking neuron model by itself depends on neurotransmitter concentration at the input stage.[18][19][20]

Fig 4: High level block diagram of the receptor layer and neuron model by Nossenson & Messer.[18][20]

Fig 5. The prediction for the firing rate in response to a pulse stimulus as given by the model by Nossenson & Messer.[18][20]

An important feature of this model is the prediction for neurons firing rate pattern which captures, using a low number of free parameters, the characteristic edge emphasized response of neurons to a stimulus pulse, as shown in Fig. 5. The firing rate is identified both as a normalized probability for neural spike firing, and as a quantity proportional to the current of neurotransmitters released by the cell. The expression for the firing rate takes the following form:

${\displaystyle R_{fire}(t)={\frac {P_{spike}(t;\Delta _{t})}{\Delta _{t}}}=[y(t)+R_{0}]\cdot P_{0}(t)}$

where,

• P0 is the probability of the neuron to be "armed" and ready to fire. It is given by the following differential equation:

${\displaystyle {\dot {P}}_{0}=-[y(t)+R_{0}+R_{1}]\cdot P_{0}(t)+R_{1}}$

P0 could be generally calculated recursively using Euler method, but in the case of a pulse of stimulus it yields a simple closed form expression.[18][21]

• y(t) is the input of the model and is interpreted as the neurotransmitter concentration on the cell surrounding (in most cases glutamate) . For an external stimulus it can be estimated through the receptor layer model:

${\displaystyle y(t)\simeq g_{gain}\cdot \langle s^{2}(t)\rangle }$ , with ${\displaystyle \langle s^{2}(t)\rangle }$ being short temporal average of stimulus power (given in Watt or other energy per time unit).

• R0 corresponds to the intrinsic spontaneous firing rate of the neuron.
• R1 is the recovery rate of the neuron from refractory state.

Other predictions by this model include:

1) The averaged Evoked Response Potential (ERP) due to population of many neurons in unfiltered measurements resembles the firing rate.[20]

2) The voltage variance of activity due to multiple neuron activity resembles the firing rate (also known as Multi-Unit-Activity power or MUA).[19][20]

3) The inter-spike-interval probability distribution takes the form a gamma-distribution like function.[18][21]

The following is a list of published non-Markovian neuron models:

• Johnson, and Swami[39]
• Berry and Meister[40]
• Kass and Ventura[41]

The theta model, or Ermentrout–Kopell canonical model, is a model originally developed to model neurons in the animal Aplysia, and is particularly well suited to describe neuron parabolic bursting.[42]

According to the model by Koch and Segev,[7] the response of a neuron to individual neurotransmitters can be modeled as an extension of the classical Hodgkin–Huxley model with both standard and nonstandard kinetic currents. Four neurotransmitters primarily have influence in the CNS. AMPA/kainate receptors are fast excitatory mediators while NMDA receptors mediate considerably slower currents. Fast inhibitory currents go through GABA_A receptors, while GABA_B receptors mediate by secondary G-protein-activated potassium channels. This range of mediation produces the following current dynamics:

• ${\displaystyle I_{\mathrm {AMPA} }(t,V)={\bar {g}}_{\mathrm {AMPA} }\cdot [O]\cdot (V(t)-E_{\mathrm {AMPA} })}$
• ${\displaystyle I_{\mathrm {NMDA} }(t,V)={\bar {g}}_{\mathrm {NMDA} }\cdot B(V)\cdot [O]\cdot (V(t)-E_{\mathrm {NMDA} })}$
• ${\displaystyle I_{\mathrm {GABA_{A}} }(t,V)={\bar {g}}_{\mathrm {GABA_{A}} }\cdot ([O_{1}]+[O_{2}])\cdot (V(t)-E_{\mathrm {Cl} })}$
• ${\displaystyle I_{\mathrm {GABA_{B}} }(t,V)={\bar {g}}_{\mathrm {GABA_{B}} }\cdot {\tfrac {[G]^{n}}{[G]^{n}+K_{\mathrm {d} }}}\cdot (V(t)-E_{\mathrm {K} })}$

where ḡ is the maximal[3][7] conductance (around 1S) and E is the equilibrium potential of the given ion or transmitter (AMDA, NMDA, Cl, or K), while [O] describes the fraction of receptors that are open. For NMDA, there is a significant effect of magnesium block that depends sigmoidally on the concentration of intracellular magnesium by B(V). For GABA_B, [G] is the concentration of the G-protein, and K_d describes the dissociation of G in binding to the potassium gates.

The dynamics of this more complicated model have been well-studied experimentally and produce important results in terms of very quick synaptic potentiation and depression, that is, fast, short-term learning.

The model by Nossenson and Messer translates neurotransmitter concentration at the input stage to the probability of releasing neurotransmitter at the output stage.[18][19][20] For a more detailed description of this model, see the Two state Markov model section above.

Further reading on this subject:[1][43]

Further reading on this subject:[44][45]

Further reading on this subject see:[46][47][48]

The neurotransmitter-based energy detection scheme[20][21] suggests that the neural tissue chemically executes a Radar-like detection procedure.

Fig. 6 The biological neural detection scheme as suggested by Nossenson et al.[20][21]

As shown in Fig. 6, the key idea of the conjecture is to account neurotransmitter concentration, neurotransmitter generation and neurotransmitter removal rates as the important quantities in executing the detection task, while referring to the measured electrical potentials as a side effect that only in certain conditions coincide with the functional purpose of each step. The detection scheme is similar to a radar like "energy detection" because it includes signal squaring, temporal summation and a threshold switch mechanism, just like the energy detector, but it also includes a unit that emphasizes stimulus edges and a variable memory length (variable memory). According to this conjecture, the physiological equivalent of the energy test statistics is neurotransmitter concentration, and the firing rate corresponds to neurotransmitter current. The advantage of this interpretation is that it leads to a unit consistent explanation which allows to bridge between electrophysiological measurements, biochemical measurements and psychophysical results.

The evidence reviewed in[20][21] suggest the following association between functionality to histological classification:

1. Stimulus squaring is likely to be performed by receptor cells.
2. Stimulus edge emphasizing and signal transduction is performed by neurons.
3. Temporal accumulation of neurotransmitters is performed by glia cells. Short term neurotransmitter accumulation is likely to occur also in some types of neurons.
4. Logical switching is executed by glia cells, and it results from exceeding a threshold level of neurotransmitter concentration. This threshold crossing is also accompanied by a change in neurotransmitter leak rate.
5. Physical all-or-non movement switching is due to muscle cells and results from exceeding a certain neurotransmitter concentration threshold on muscle surroundings.

Note that although the electrophysiological signals in Fig.6 are often similar to the functional signal (signal power / neurotransmitter concentration / muscle force), there are some stages in which the electrical observation is different from the functional purpose of the corresponding step. In particular, Nossenson et al. suggested that glia threshold crossing has a completely different functional operation compared to the radiated electrophysiological signal, and that the latter might only be a side effect of glia break.