Open In Colab

Neuron Coding

Neurons are able to coordinate movement of sets of muscles because they are synaptically connected to each other into circuits. The way a neural circuit functions depends on two main things:

  1. How neurons are connected (types of synapses and identity of connections).

  2. Intrinsic physiology of individual neurons.

Neurons are quite diverse. Differences in basic intrinsic properties effect how they are recruited by the circuits in which they are integrated. Intrinsic neuron properties therefore effect the patterned output generated by neural circuits.

We can explore the effects of some basic intrinsic properties on neural recruitment and spike coding by using neuron “models.” In this notebook you will work with simulations of the classic and very basic Leaky Integrate-and-Fire (LIF) neuron model. You will study response dynamics to various types of applied current (ie. input).

We will especially emphasize identifying conditions (input statistics) under which a neuron can either: 1. spike at low firing rates and in an irregular manner, or 2. spike at a precise time or in a regular manner. The reason for focusing on this is that different parts of the motor system need to meet different coding requirements.

The Leaky Integrate-and-Fire (LIF) model

This supplementary video (https://youtube.com/watch?v=rSExvwCVRYg) introduces the reduction of a biological neuron to a simple leaky-integrate-fire (LIF) neuron model. It is more detailed than you need to know to complete this notebook, but helpful if you are interested in understanding the process more in depth.

The basic idea of LIF neuron was proposed in 1907 by Louis Édouard Lapicque, long before we understood the electrophysiology of a neuron (see a translation of Lapicque’s paper ). More details of the model can be found in the book Theoretical neuroscience by Peter Dayan and Laurence F. Abbott.

The subthreshold membrane potential dynamics of a LIF neuron is described by

(1)CmdVdt=gL(VEL)+I,(1)

where V is the membrane potential, Cm is the membrane capacitance, gL is the leak conductance (gL=1/R, the inverse of the leak resistance R), EL is the resting potential (the equilibrium potential of the leak conductance), and I is an applied input current. Dividing both sides of the above equation by gL gives

(2)τmdVdt=(VEL)+IgL,(2)

where the τm is membrane time constant and is defined as τm=Cm/gL.

Note that dividing capacitance by conductance gives units of time!

The LIF neuron model uses Eqn.(2) to simulate membrane potential dynamics. The model implements a shortcut to simulating spiking activity:

If I is sufficiently strong such that V reaches a certain threshold value Vth, V is reset to a reset potential Vreset<Vth, and voltage is clamped to Vreset for τref ms, mimicking the refractoriness of the neuron during an action potential:

The LIF model captures the facts that a neuron:

  • performs spatial and temporal integration of synaptic inputs

  • generates a spike when the voltage reaches a certain threshold

  • goes refractory during the action potential

  • has a leaky membrane

However, the LIF model assumes that the spatial and temporal integration of inputs is linear (which is not strictly true in real neurons). Also, membrane potential dynamics close to the spike threshold are much slower in LIF neurons than in real neurons.

Setup

Click to show
#@title {display-mode: "form"}

#@markdown Execute this cell to initialize the notebook
#@markdown and define functions required to simulate an LIF model neuron. 

# Imports
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets  # interactive display
from IPython.display import display
from datetime import datetime,timezone,timedelta
%config InlineBackend.figure_format = 'retina'
# use NMA plot style
plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle")
my_layout = widgets.Layout()

def plot_volt_trace(pars, v, sp):
    """
    Plot trajetory of membrane potential for a single neuron

    Expects:
    pars   : parameter dictionary
    v      : volt trajetory
    sp     : spike train

    Returns:
    figure of the membrane potential trajetory for a single neuron
    """

    V_th = pars['V_th']
    dt, range_t = pars['dt'], pars['range_t']
    if sp.size:
        sp_num = (sp / dt).astype(int) - 1
        v[sp_num] += 40  # draw nicer spikes
    plt.clf()
    plt.plot(pars['range_t'], v, 'b')
    plt.axhline(V_th, 0, 1, color='k', ls='--')
    plt.xlabel('Time (ms)',fontsize=14)
    plt.ylabel('V (mV)',fontsize=14)
    plt.xticks(fontsize=14)
    plt.yticks(fontsize=14)
    plt.legend(['Membrane\npotential', r'Threshold V$_{\mathrm{th}}$'],
             loc=[1.05, 0.75])
    plt.ylim([-90, -20])
    plt.show()


def plot_GWN(pars, I_GWN):
    """
    Args:
    pars  : parameter dictionary
    I_GWN : Gaussian white noise input

    Returns:
    figure of the gaussian white noise input
    """

    plt.figure(figsize=(12, 4))
    plt.subplot(121)
    plt.plot(pars['range_t'][::3], I_GWN[::3], 'b')
    plt.xlabel('Time (ms)',fontsize=14)
    plt.ylabel(r'Current (pA)',fontsize=14)
    plt.xticks(fontsize=14)
    plt.yticks(fontsize=14)
    plt.subplot(122)
    plot_volt_trace(pars, v, sp)
    plt.tight_layout()
    plt.show()

def plot_IF_curve(spk_count, spk_count_dc, I_mean, sig_gwn):
    """
    Plot trajetory of membrane potential for a single neuron

    Expects:
    pars   : parameter dictionary
    v      : volt trajetory
    sp     : spike train

    Returns:
    figure of the membrane potential trajetory for a single neuron
    """
#     hfig,ax = plt.subplots(1)
    plt.clf()
    plt.plot(I_mean, spk_count, 'k',
           label=r'noise=%.2f' % sig_gwn)
    plt.plot(I_mean, spk_count_dc, 'k--', alpha=0.5, lw=4, dashes=(2, 2),
           label='constant current')
    plt.ylabel('Spike count',fontsize=14)
    plt.xlabel('Average injected current (pA)',fontsize=14)
    plt.xticks(fontsize=14)
    plt.yticks(fontsize=14)
    plt.legend(loc='best',fontsize=14)
    plt.show()

def my_hists(isi1, isi2, cv1, cv2, sigma1, sigma2):
    """
    Args:
    isi1 : vector with inter-spike intervals
    isi2 : vector with inter-spike intervals
    cv1  : coefficient of variation for isi1
    cv2  : coefficient of variation for isi2

    Returns:
    figure with two histograms, isi1, isi2

    """
    plt.figure()
    my_bins = np.linspace(10, 30, 20)
    plt.subplot(121)
    plt.hist(isi1, bins=my_bins, color='b', alpha=0.5)
    plt.xlabel('ISI (ms)',fontsize=14)
    plt.ylabel('count',fontsize=14)
    plt.xticks(fontsize=14)
    plt.yticks(fontsize=14)
    plt.title(r'$\sigma_{GWN}=$%.1f, CV$_{\mathrm{isi}}$=%.3f' % (sigma1, cv1),fontsize=14)

    plt.subplot(122)
    plt.hist(isi2, bins=my_bins, color='b', alpha=0.5)
    plt.xlabel('ISI (ms)')
    plt.ylabel('count')
    plt.xticks(fontsize=14)
    plt.yticks(fontsize=14)
    plt.title(r'$\sigma_{GWN}=$%.1f, CV$_{\mathrm{isi}}$=%.3f' % (sigma2, cv2))
    plt.tight_layout()
    plt.show()

def default_pars(**kwargs):
    pars = {}

    # typical neuron parameters#
    pars['V_th'] = -55.     # spike threshold [mV]
    pars['V_reset'] = -75.  # reset potential [mV]
    pars['tau_m'] = 10.     # membrane time constant [ms]
    pars['g_L'] = 10.       # leak conductance [nS]
    pars['V_init'] = -75.   # initial potential [mV]
    pars['E_L'] = -75.      # leak reversal potential [mV]
    pars['tref'] = 2.       # refractory time (ms)

    # simulation parameters #
    pars['T'] = 400.  # Total duration of simulation [ms]
    pars['dt'] = .1   # Simulation time step [ms]

    # external parameters if any #
    for k in kwargs:
        pars[k] = kwargs[k]

    pars['range_t'] = np.arange(0, pars['T'], pars['dt'])  # Vector of discretized time points [ms]

    return pars

def run_LIF(pars, Iinj, stop=False, durstep=100):
    """
    Simulate the LIF dynamics with external input current
    Args:
    pars       : parameter dictionary
    Iinj       : input current [pA]. The injected current here can be a value
                 or an array
    stop       : boolean. If True, use a current pulse
    Returns:
    rec_v      : membrane potential
    rec_sp     : spike times
    """

    # Set parameters
    V_th, V_reset = pars['V_th'], pars['V_reset']
    tau_m, g_L = pars['tau_m'], pars['g_L']
    V_init, E_L = pars['V_init'], pars['E_L']
    dt, range_t = pars['dt'], pars['range_t']
    fs = 1/dt
    Lt = range_t.size
    tref = pars['tref']

    # Initialize voltage
    v = np.zeros(Lt)
    v[0] = V_init

    # Set current time course
    Iinj = Iinj * np.ones(Lt)
    durstep = int((durstep*fs)/2)
    # If current pulse, set beginning and end to 0
    if stop:
        # Iinj[:int(len(Iinj) / 2) - (int(durstep/pars['dt'])/2)] = 0
        # Iinj[int(len(Iinj) / 2) + (int(durstep/pars['dt'])/2):] = 0

        Iinj[:int(len(Iinj) / 2) - durstep] = 0
        Iinj[int(len(Iinj) / 2) + durstep:] = 0

    # Loop over time
    rec_spikes = []  # record spike times
    tr = 0.  # the count for refractory duration

    for it in range(Lt - 1):

        if tr > 0:  # check if in refractory period
            v[it] = V_reset  # set voltage to reset
            tr = tr - 1 # reduce running counter of refractory period

        elif v[it] >= V_th:  # if voltage over threshold
            rec_spikes.append(it)  # record spike event
            v[it] = V_reset  # reset voltage
            tr = tref / dt  # set refractory time

        # Calculate the increment of the membrane potential
        dv = (-(v[it] - E_L) + Iinj[it] / g_L) * (dt / tau_m)

        # Update the membrane potential
        v[it + 1] = v[it] + dv

    # Get spike times in ms
    rec_spikes = np.array(rec_spikes) * dt

    return v, rec_spikes

def my_GWN(pars, mu, sig, myseed=False):
    """
    Function that generates Gaussian white noise input

    Args:
    pars       : parameter dictionary
    mu         : noise baseline (mean)
    sig        : noise amplitute (standard deviation)
    myseed     : random seed. int or boolean
                 the same seed will give the same
                 random number sequence

    Returns:
    I          : Gaussian white noise input
    """

    # Retrieve simulation parameters
    dt, range_t = pars['dt'], pars['range_t']
    Lt = range_t.size

    # Set random seed
    if myseed:
        np.random.seed(seed=myseed)
    else:
        np.random.seed()

    # Generate GWN
    # we divide here by 1000 to convert units to sec.
    I_gwn = mu + sig * np.random.randn(Lt) / np.sqrt(dt / 1000.)

    return I_gwn

# help(my_GWN)

def isi_cv_LIF(spike_times):
  """
  Calculates the inter-spike intervals (isi) and
  the coefficient of variation (cv) for a given spike_train
  Args:
    spike_times : (n, ) vector with the spike times (ndarray)
  Returns:
    isi         : (n-1,) vector with the inter-spike intervals (ms)
    cv          : coefficient of variation of isi (float)
  """

  if len(spike_times) >= 2:
    # Compute isi
    isi = np.diff(spike_times)
    # Compute cv
    cv = isi.std()/isi.mean()
  else:
    isi = np.nan
    cv = np.nan

  return isi, cv

print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

Intrinsic physiology.

Interactive Demo 1

Running the following code cell initiates a demo that enables you to explore how intrinsic physiology effects neuron recruitment and spike rate. In the simulation, a constant current (DC) input is provided to the model neuron.

Make sure you understand all of the parameters in the simulation while exploring.

You can specify the duration of the simulation, the amount of current injected, and whether the current injection is continuous throughout the simulation or a “step” in the middle.

To understand the effects of each parameter, it is helpful to change just one at a time. How many different ways can you change the response of the neuron (whether it spikes or not and at what rate)?

Click to show
#@title {display-mode: "form"}

#@markdown Run this cell to enable the interactive demo.
#@markdown Use the sliders to set the value of each parameter of the simulation. 


my_layout.width = '700px'
# my_layout.description_width = 'initial'
style = {'description_width': 'initial'}
@widgets.interact(
    applied_current=widgets.FloatSlider(50., min=0., max=1000., step=2.,
                               layout=my_layout,style=style),
    Simulation_Duration=widgets.FloatSlider(400., min=0., max=1000., step=10.,
                               layout=my_layout,style=style),
    Injection_Step=widgets.Checkbox(value=True,
                               description='current step (vs continuous)',
                               layout=my_layout,style=style),
    Current_Step_Duration=widgets.FloatSlider(100., min=10, max=200., step=2.,
                               layout=my_layout,style=style),
    Leak_Reversal=widgets.FloatSlider(-75., min=-90, max=-30., step=2.,
                               layout=my_layout,style=style),
    Leak_Conductance=widgets.FloatSlider(10., min=1, max=50., step=2.,
                               layout=my_layout,style=style),
    AHP_Voltage=widgets.FloatSlider(-80., min=-90., max=-30.,step=2.,
                               layout=my_layout,style=style),
    Spike_Threshold=widgets.FloatSlider(-55., min=-100., max=-30., step=2.,
                               layout=my_layout,style=style),
    Membrane_Tau=widgets.FloatSlider(5., min=1., max=15., step=1.,
                               layout=my_layout,style=style)
)


def diff_DC(applied_current,
            Simulation_Duration,
            Injection_Step,
            Current_Step_Duration,
            Leak_Reversal,
            Leak_Conductance,
            AHP_Voltage,
            Spike_Threshold,
            Membrane_Tau): #I_dc=200., tau_m=10.):
    # pars = default_pars(T=100.)
    # pars['tau_m'] = tau_m
    pars = default_pars(
      E_L = Leak_Reversal,
      g_L = Leak_Conductance,
      V_init = Leak_Reversal,
      V_reset = AHP_Voltage,
      V_th = Spike_Threshold,
      tau_m = Membrane_Tau,
      T=Simulation_Duration)
    # pars=default_pars()
    v, sp = run_LIF(pars, Iinj=applied_current,stop=Injection_Step,durstep=Current_Step_Duration)
    plot_volt_trace(pars, v, sp)

print('Interactive demo initiated at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

If you are a bit stuck starting your exploration, here are some ideas to get you started. Start with “current step (versus continuous)” checked (True/Yes).

  • How much DC input current is needed to reach the threshold (this is called the rheobase current)?

  • Once you apply enough current to get the neuron just above spike threshold, what is the effect of changing the:

    • membrane time constant?

    • AHP voltage?

    • leak reversal?

  • Find an input current amplitude that drives the neuron to a voltage about halfway between the leak reversal and the spike threshold. What is the effect of changing the:

    • leak conductance?

Input “noise”

Given the noisy nature of neuronal activity in vivo, neurons usually receive complex, time-varying inputs rather than steady “DC” inputs. In this section, you will use the model neuron to explore the effect of input current noise on the spiking response.

We will mimic membrane current noise as gaussian white noise (random fluctuations with a standard deviation equal to the noise magnitude).

Interactive Demo 2

The noise describes only the fluctuations of the input received by a neuron. The input current also has a mean amplitude (as it did before). Indeed, when the noise = 0, the input current is just constant as it was before.

So the question arises how does membrane current noise affect the spiking behavior of the neuron. For instance we may want to know:

  1. How does the minimum input current needed to make a neuron spike change with increase in noise? And how does this effect rate coding?

  2. How does the spike regularity change with increase in noise? And how does this effect temporal coding across repetitions of brief (~50msec) excitatory input?

This demo enables you to explore these questions by changing the amplitude and fluctuation sizes (noise) of the membrane current. The input current is applied at time=0 for the simulation duration.

Click to show
#@title {display-mode: "form"}

#@markdown Run this cell to run the interactive demo.  
#@markdown >NOTE: With any given set of parameters, you can "repeat stimulation" 
#@markdown by hitting that button. Because noise is "random" the spiking can be
#@markdown different on different trials even with the same settings. 


my_layout.width = '450px'

@widgets.interact(
    amplitude=widgets.FloatSlider(200., min=0., max=300., step=5., layout=my_layout),
    noise=widgets.FloatSlider(2.5, min=0., max=5., step=.1, layout=my_layout),
    repeat_trial=widgets.ToggleButton(value=False, description='Repeat Simulation')
    )

def diff_GWN_to_LIF(amplitude,noise,repeat_trial):
    pars = default_pars(T=200.)
    I_GWN = my_GWN(pars, mu=amplitude, sig=noise)
#     pars['tau_m'] = tau_m
    v, sp = run_LIF(pars, Iinj=I_GWN,stop=False,durstep=100)
    plot_volt_trace(pars, v, sp)

# @widgets.interact(
#     amplitude=widgets.FloatSlider(200., min=0., max=300., step=5.,
#                                layout=my_layout),
#     noise=widgets.FloatSlider(2.5, min=0., max=5., step=.1,
#                                 layout=my_layout),
#     step=widgets.Checkbox(value=False, description='input a discrete step',
#                                 layout=my_layout),
#     Current_Step_Duration=widgets.FloatSlider(100., min=10, max=200., step=2.,
#                                layout=my_layout,style=style),
#     repeat_trial=widgets.ToggleButton(value=False,
#                                       description='Repeat Simulation')
#     )

# def diff_GWN_to_LIF(amplitude, noise,step,Current_Step_Duration,repeat_trial):
#     pars = default_pars(T=400.)
#     I_GWN = my_GWN(pars, mu=amplitude, sig=noise)
# #     pars['tau_m'] = tau_m
#     v, sp = run_LIF(pars, Iinj=I_GWN,stop=step,durstep=Current_Step_Duration)
#     plot_volt_trace(pars, v, sp)

print('Interactive demo initiated at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

Interactive Demo 3

As you saw in the Interactive Demo #2, the noisiness of the membrane current effects the spike count/rate.

When we plot spike rate as a function of membrane current, it is called the input-output transfer function of the neuron (ie F-I curve: “firing rate”-“input current”). In this simulation, the number of spikes evoked by a 1-second long current step are plotted as a function of current amplitude.

This demo enables you to explore the following questions by changing the magnitude of the input noise.

  1. Does the membrane current noise always effect the spike count, or does its effects depend on the amplitude of the current?

  2. How does the F-I curve of the LIF neuron change as we increase the membrane current noise?

Click to show
#@title {display-mode: "form"}

#@markdown Run this cell to run the interactive demo.  

my_layout.width = '450px'
@widgets.interact(
    noise=widgets.FloatSlider(3.0, min=0., max=6., step=0.5,
                                layout=my_layout)
)

def diff_std_affect_fI(noise):
    pars = default_pars(T=1000.)
    I_mean = np.arange(100., 400., 10.)
    spk_count = np.zeros(len(I_mean))
    spk_count_dc = np.zeros(len(I_mean))

    for idx in range(len(I_mean)):
        I_GWN = my_GWN(pars, mu=I_mean[idx], sig=noise, myseed=2020)
        v, rec_spikes = run_LIF(pars, Iinj=I_GWN)
        v_dc, rec_sp_dc = run_LIF(pars, Iinj=I_mean[idx])
        spk_count[idx] = len(rec_spikes)
        spk_count_dc[idx] = len(rec_sp_dc)
    
    plot_IF_curve(spk_count, spk_count_dc, I_mean, noise)

print('Interactive demo initiated at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

Interactive Demo 4

As you saw in the Interactive Demo #2, the noisiness of the membrane current also effects the temporal precision (and regularity) of spiking. Spiking regularity is often quantified by the Coefficient of Variation of the inter-spike-interval (CVISI). A neurons’s “inter-spike-interval” (ISI) is the time between spikes (and the inverse of its spike rate). The CVISI is calculated by:

(3)CVISI=std(ISI)mean(ISI)

For a clocklike (regular) ISI, CVISI= 0 because std(ISI)=0.

This demo enables you to explore the following questions by changing the magnitude of the input noise.

  1. Is there a difference between the effect of small injected current amplitude and high injected current amplitude and the CV? How does this difference depend on the magnitude of membrane current noise? Why do you think that is?

  2. Why does increasing the amplitude of current injection decrease the CVISI?

  3. If you plot spike count (or rate) vs. CVISI, should there be a relationship between the two?

Click to show
#@title {display-mode: "form"}

#@markdown Run this cell to run the interactive demo. 

my_layout.width = '450px'
@widgets.interact(
    noise=widgets.FloatSlider(0.0, min=0., max=5.,
                                step=0.1, layout=my_layout)
)


def diff_std_affect_fI(noise):
    pars = default_pars(T=1000.)
    I_mean = np.arange(100., 400., 20)
    spk_count = np.zeros(len(I_mean))
    cv_isi = np.empty(len(I_mean))

    for idx in range(len(I_mean)):
        I_GWN = my_GWN(pars, mu=I_mean[idx], sig=noise)
        v, rec_spikes = run_LIF(pars, Iinj=I_GWN)
        spk_count[idx] = len(rec_spikes)
        if len(rec_spikes) > 3:
            isi = np.diff(rec_spikes)
            cv_isi[idx] = np.std(isi) / np.mean(isi)

    # Plot the F-I curve i.e. Output firing rate as a function of input mean.
    plt.figure()
    plt.clf()
    plt.plot(I_mean[spk_count > 5], cv_isi[spk_count > 5], 'bo', alpha=0.5)
    plt.xlabel('Average injected current (pA)',fontsize=14)
    plt.xticks(fontsize=14)
    plt.ylabel(r'Spike irregularity ($\mathrm{CV}_\mathrm{ISI}$)',fontsize=14)
    plt.yticks(fontsize=14)
    plt.ylim(-0.1, 1.5)
    plt.xlim(50, 400)
    plt.grid(True)
    plt.show()

print('Interactive demo initiated at ' + str(datetime.now(timezone(-timedelta(hours=5)))))

Summary

Congratulations! You’ve just used a classic integrate-and-fire (LIF) neuron model simulation to study the dynamics of neural activity. You drove the LIF neuron with external inputs, such as direct current and Gaussian white noise. You examined how different inputs affect the LIF neuron’s output (firing rate and spike time irregularity). Throughout the course, think about what types of neural coding are required by various neural circuits for motor control. You can now start to relate neural coding properties back to the basic intrinsic physiology of neurons. There are so many other anatomical and physiological properties that effect neural coding that the ones you examined here. This is just scratching the surface on neural complexity in circuits, coding and motor control.


Created by Dr. Krista Perks for a Motor Systems course (BIOL358) at Wesleyan University based on ‘Neuromatch Academy” Week 2, Day 3: Biological Neuron Models

Original Notebook Tutorial Information:

Content creators: Qinglong Gu, Songtin Li, John Murray, Richard Naud, Arvind Kumar

Content reviewers: Maryam Vaziri-Pashkam, Ella Batty, Lorenzo Fontolan, Richard Gao, Matthew Krause, Spiros Chavlis, Michael Waskom, Ethan Cheng