Data Explorer
Contents
Data Explorer¶
Table of Contents¶
Computer Models¶
Modelling passive and active properties of neurons based on electric circuit models from Week 1.
Setup¶
Import and define functions
#@title {display-mode: "form"}
#@markdown Run this code cell to import packages and define functions
from pathlib import Path
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import pandas as pd
from datetime import datetime,timezone,timedelta
pal = sns.color_palette(n_colors=15)
pal = pal.as_hex()
from ipywidgets import interactive, HBox, VBox, widgets, interact
import ipywidgets as widgets # interactive 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 hyper_fit(t,r,c):
return r + ((r*c)/t)
print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
Computer implementation of computational models¶
To implement this computational model in computer code, you need to define a function.
A function in computer programming is executable code (it will do things when you “run” it). Functions take information that they need to do what they need to do.
As an analogy, your computational model is a function (an equation in this case). You executed that function when you were given values for ion equilibrium potential and resistance of a neuron and calculated the neuron’s resting membrane potential.
To define a function in python so that the computer can execute it, we need to write down instructions in a specific way:
def function_name( parameters ):
variable = equation
return [variable]
def indicates that we are defining a function
function_name is the name we give the function
in parentheses after the function name, we list the information that will be provided to the function when it is executed (these are info that will be needed to do what we want the function to do… in this case calculate the resting membrane potential of a neuron)
:transitions the code into the instructions for the function to follow when it is executedthe instructions for the function to execute are on the rest of the written lines, with each line indented
return terminates the function and provides information that can be accessed after the function executes
Let’s consider a basic function that calculates the square of a number (the number times itself).
Examine the contents of the following code cell and then run the code cell to execute the definition of the function (this step actually creates the function that is defined with the text)
def square(x): # this sets the name of the function and any variables that it will use
result = x*x # this is the equation that this function uses to transform 'x' input into 'result' output
return result # this line makes the variable 'result' available as an output of the function
Nothing obvious will happen when you run a code cell that defines a function (other than an indication that the code cell ran). But… you can now use the function that you defined.
In the code cell below, assign a number to the varable my_number. Then, run the code cell to execute that variable assignment.
my_number = 2
Now, in the next code cell, we will execute that function by typing the name of the function followed by parentheses containing the variable my_number (to provide the function with that information).
Examine the code cell below and then run it.
square(my_number)
The function executes and returns the result of the instructions that it followed. In this case, the result is printed below the code cell, but can be assigned to another variable name by typing result = square(my_number).
Steady State Model¶
You created a computational model of the steady state membrane potential that depends on the equilibrium potential and resistance of two ion conductances. To implement this computational model in computer code, you need to define a function that computes Vin based on your computational model equation. In the following code cell, replace ... with the equation from your computational model that defines \(V_{in}\).
# define the function
def simulate_Vin_steady_state(E1,R1,E2,R2):
Vin = ...
return Vin
Click here for solution
#define the function
def simulate_Vin_steady_state(E1,R1,E2,R2):
Vin = ... #( (E1/R1) + (E2/R2) ) / ( (1/R1) + (1/R2) )
return Vin
In the code cell below, assign the values that you used in your circuit membrane model (Passive Membrane Properties lab: Part I) to each variable needed by the function. Then, run the code cell to assign those values to the variables.
R1 = ...
R2 = ...
E1 = ...
E2 = ...
Now, examine the contents of the code cell below and replace the ... with the model parameters you just defined (separated by commas). Then run the code cell to simulate your computational model given the values you just provided. Note that the order of information provided to the function matters - it must match the order you defined when defining the function.
simulate_Vin_steady_state( ... )
Click here for solution
simulate_Vin_steady_state(E1,R1,E2,R2)
How do the results of your computational model compare to your manually measured electric circuit model results from that lab?
Interactive Simulation 1: Steady-State Membrane Potential Model¶
Let’s implement your computational model into a more interactive simulation that you can use to efficiently explore the predictions of your model. There is more to learn about computer programming in Python than we will cover here to be able to make a more interactive simulation, so I have done the coding for you. But you can click Show Code if you want to take a look if you want.
Run the code cell below to run the interactive simulation. Toggle the sliders to explore the behavior of the model across the parameter space.
#@title {display-mode: 'form'}
#@markdown Run this code cell to generate an interactive widget using the scripted
#@markdown implementation (of your passive membrane model) that you just created. <br>
#@markdown The results of the simulation (given the selected paramaters) is printed on the top line.
#@markdown (note that parameter selection actually happens when you release the mouse click from the slider).
slider_E1 = widgets.FloatSlider(
min=-100,
max=100,
value=0,
step= 1,
readout=True,
continuous_update=False,
description='E1 (mV)')
slider_E1.layout.width = '600px'
slider_R1 = widgets.FloatSlider(
min=1,
max=100,
value=1,
step= 1,
readout=True,
continuous_update=False,
description='R1 (MOhm)')
slider_R1.layout.width = '600px'
slider_E2 = widgets.FloatSlider(
min=-100,
max=100,
value=0,
step= 1,
readout=True,
continuous_update=False,
description='E2 (mV)')
slider_E2.layout.width = '600px'
slider_R2 = widgets.FloatSlider(
min=1,
max=100,
value=1,
step= 1,
readout=True,
continuous_update=False,
description='R2 (MOhm)')
slider_R2.layout.width = '600px'
label_result = widgets.Label(
value='The resting membrane potential of the neuron will be: '
)
label_result.layout.width = '600px'
display(label_result)
# a function that will modify the xaxis range
def update_result(E1,R1,E2,R2):
Vin = simulate_Vin_steady_state(E1,R1,E2,R2)
label_result.value= f'The resting membrane potential of the neuron will be: {Vin:.2f} mV'
w = interact(update_result, E1=slider_E1,R1 = slider_R1,E2=slider_E2,R2=slider_R2);
Dynamic Model¶
The dynamic computatinoal model that we formulated in class is a differential equation. We can simulate a differential equation by iterating through each time step across an array of time.
Remember that your computational model is in a format useful for calculating the change in voltage (\(\partial{V}\)) across some duration of time:
where \(Iapp\) is the current applied across the membrane, \(V\) is the current membrane potential (\(V_{in}\)), \(E\) is the net equilibrium potential of the membrane at that moment in time, \(R\) is the net resistance of the membrane at that moment in time, and \(C\) is the membrane capacitance. This equation can be evaluated over any value of \(\partial{t}\).
Now let’s look at how to script a function that iterates through sequential time steps and uses the computational model to determine the membrane potential of the neuron at each time step. In this way, the script defines a function that will implements a simulation of the computational model.
def simulate_Vin_dynamic(E, R, C, Iapp, duration, dt):
'''
The simulation starts at time = 0 with the voltage given by the input to the function (V)
E = the equilibrium potential of the net ionic conductance
R = the resistance of the membrane in Ohm
C = the capacitance of the membrane in Farads
Iapp = the current (microAmps) applied to the system throughout the simulation
duration = the total time (seconds) that the current is applied to the system (in this case also the duration of the simulation)
dt = the time step (seconds) for analysing the simulation (how much time elapses between each analysis of dV)
'''
timesteps = int(duration/dt)
'''
we will just use the equilibrium membrane potential of the neuron
as the initial membrane potential
'''
V_record = [E] # initialize an array of membrane potential across time with the initial membrane potential
''' then, we will iterate through time using a "for loop" '''
for i in range(timesteps): # for each time step, iterate through the simulation
V = V_record[-1] # what is the current membrane potential?
dV = (Iapp - (1/R)*(V-E))*(dt/C) # how much does the voltage change by on this time step?
V = (V+dV) # add the change in voltage to the current voltage
V_record.append(V) # append the new voltage to the array of membrane voltage across time.
return V_record # make the array "V_record" available as an output of the scripted function
# Make sure to "run" this code cell before moving on,
# so that the function "simulate_Vin_dynamic" can be used by later code cells.
To run a simulation of the model, we need to define: 1) the initial conditions of the membrane, and 2) the amplitude of the applied current.
Before using the scripted version of your computational model to simulate the membrane potential response to an applied current, it is helpful to manually calculate the response for just a few timesteps after a current is applied.
Use the following initial conditions: E = -60 millivolts, R = 100 MOhm, C = 100 pF. Apply a 100 pA current across the membrane.
Time Step 1
What is the inital value of \(V_{in}\)?
This is the value of \(V_{in}\) at the moment when the current is applied.
At the first time step, what is the value of (\(\partial{V}\))?
Update \(V_{in}\) according to the initial \(V_{in}\) + (\(\partial{V}\)) on the first time step.
Time Step 2
What is the initial value of \(V_{in}\) on this time step?
At this time step, what is the value of (\(\partial{V}\))?
WHY is the value of (\(\partial{V}\)) different for this time step versus the last time step?
Note there are several ways to address why. There is an answer in terms of the value of \(V_{in}\). There is also an answer in terms of the current accross the capacitor. Remember, voltage across a capacitor only changes when there is current across the capacitor (and the magnitude of the voltage change is proportional to the amount of current).
Update \(V_{in}\) according to the inital \(V_{in}\) + (\(\partial{V}\)) on the first time step.
Time Step 3
Repeat the steps for time step 2 and notice how/why (\(\partial{V}\)) changes again on this time step compared to the first two time steps.
In the code cell below, we are defining the net equilibrium potential of the resting membrane conductances (\(E\)), the net resistance across the membrane (\(R\)), and the capacitance of the membrane (\(C\)). We are then defining the magnitude of the applied current (\(Iapp\)). We will also specify a duration for the square pulse stimulus (in this case equal to the duration of our simulation).
Importantly, we need to set a timestep value (\(\partial{t}\)) the time between iterations of the simulation across time).
After all of the variables required by the scripted function are defined, the script run a simulation of the computational model.
After you run the code cell below, you will see a plot of the simulated membrane potential across time in response to the applied square-wave current step. The calculated \(\tau\) for the model membrane will be printed above the plot.
What happens when you increase or decrease the value of \(\partial{t}\)? (try it out)
# define the initial conditions
E = -60*1e-3 # millivolts
R = 100*1e6 # MOhm
C = 100*1e-12 # pF
Iapp = 100*1e-12 #pA
duration = 0.1
dt = 0.1*1e-3 #msec
'''
Now we run the simulation
- V_record is the output of the function (the array of membrane potential values across time steps)
'''
V_record = simulate_Vin_dynamic(E, R, C, Iapp, duration, dt)
'''
Now we plot the results of the simulation (and label the axes)
'''
hfig, ax = plt.subplots(figsize=(10,5))
ax.plot(np.linspace(0,duration,int(duration/dt)+1),V_record,color='black');
ax.set_xlabel('seconds');
ax.set_ylabel('Volts');
'''
Finally, we can calculate the time constant of this model membrane and print the value above the figure
'''
print(f'tau = {(R*C):0.4f} seconds')
Interactive Simulation 2: Dynamic Membrane Potential Model¶
Let’s implement your computational model into a more interactive simulation that you can use to efficiently explore predictions of your model.
#@title {display-mode: "form"}
#@markdown Run this code cell to genrate an interactive widget using the model you just created. <br>
#@markdown The stimulus is a square-wave current pulse starting at t=0 with a duration of 1 second
# set up the applied current
delay = 0.1 #seconds
duration = 1 #seconds
# define the simulation conditions
total_duration = 1.2 #seconds
dt = 0.00001 #seconds
# Create sliders for the figure widget inputs
E_slider = widgets.FloatSlider(
min=-100,
max=100,
value=0,
step= 1,
readout=True,
continuous_update=False,
description='Equilibrium Potential (mV)',
style = {'description_width': '200px'})
E_slider.layout.width = '600px'
R_slider = widgets.FloatSlider(
min=0,
max=1000,
value=100,
step= 1,
readout=True,
continuous_update=False,
description='R (MOhm)')
R_slider.layout.width = '600px'
C_slider = widgets.FloatSlider(
min=0.001,
max=1,
value=0.1,
step= 0.001,
readout=True,
readout_format='.3f',
continuous_update=False,
description='C (nF)')
C_slider.layout.width = '600px'
thresh_slider = widgets.FloatSlider(
min=0,
max=30,
value=10,
step= 0.5,
readout=True,
continuous_update=True,
description='Relative Spike Threshold (mV)',
style = {'description_width': '200px'})
thresh_slider.layout.width = '600px'
amplitude_slider = widgets.FloatSlider(
min=-500,
max=500,
value=100,
step= 1,
readout=True,
continuous_update=True,
description='I (pA)')
amplitude_slider.layout.width = '600px'
label_tau = widgets.Label(
value='tau = '
)
label_tau.layout.width = '600px'
display(label_tau)
def update_plot(E,R,C,thresh,amplitude):# set up the simulation items
E = E*1e-3 # convert from millivolts to volts
# V = E # start the simulation at equilibrium
R = R*1e6 # convert from MegaOhm to Ohm
# tau = tau*1e-3 # convert from milliseconds to seconds
amplitude = amplitude*1e-12 # Convert from picoAmps to Amps
C = C*1e-9 # convert from microfarad to farad
# print(R,tau,C,amplitude)
time = np.linspace(0,total_duration,int(total_duration/dt)+1)-delay
stimulus = np.zeros((int(total_duration/dt)))
stimulus[int(delay/dt):int((delay+duration)/dt)]=amplitude
V_record = np.empty((int(total_duration/dt)+1))
# write the "for" loop
# simulate_V(E, R, C, Iapp, duration, dt)
# for i,current in enumerate(stimulus):
# dV = dt*((current - (V-E)/R)/C)
# V = (V+dV)
# V_record[i]=V
V_record[0]=E
for i,Iapp in enumerate(stimulus): # for each time step, iterate through the simulation
V = V_record[i]
dV = (Iapp - (1/R)*(V-E))*(dt/C) # how much does the voltage change by on this time step?
V = (V+dV) # add the change in voltage to the current voltage
V_record[i+1]=V
# V_record = [V_initial]
# for i in range(timesteps): # for each time step, iterate through the simulation
# V = V_record[-1]
# dV = (Iapp - (1/R)*(V-0))*(dt/C) # how much does the voltage change by on this time step?
# V = (V+dV) # add the change in voltage to the current voltage
# V_record.append(V)
# return V_record
hfig,ax = plt.subplots(figsize=(10,4))
# add simulation vectors to the figure widget
ax.plot(time,V_record*1e3) # convert from Volts to milliVolts
ax.axhline((V_record[0]+thresh), time[0], time[-1], color='k', ls='--')
ax.set_xlabel('seconds')
ax.set_ylabel('mV')
label_tau.value=f'tau = {(R*C):0.4f} seconds'
w = interact(update_plot,
E=E_slider,
R=R_slider,
C=C_slider,
thresh=thresh_slider,
amplitude=amplitude_slider);
Note that \( 1mV = 100pA * 10MOhm \). Why would you want to measure and report in terms of these units and not V, A, Ohm?
Spiking¶
Simulating spiking output of your model membrane does not require mathematically solving for the voltage-gated conductances. Many computational models of neuron spiking activity simply use knowledge about spike thresholds and refractory periods to implement conditional statements in the model that simulate spikes. These are often called Leaky Integrate and Fire models (LIF).
The following gives a hybrid english-python description of the implementation:
spike_treshold = x (in mV)
spike_amplitude = x (in mV)
spike_refractory = x (in ms)
after_hyperpolarization = x (in mV)
start at t = 0
dt = x (in ms)
V(t = 0) = Equilibrium Potential
repeat the following steps until the end of time is reached:
Determine the current membrane potential
> Vinitial = V(t)
If the current membrane potential is greater than or equal to spike_treshold:
Instead of calculating dV for this timestep,
set the voltage at the next timestep equal to the amplitude of a spike
> V(t + dt) = spike_amplitude
Move on to the next time step after the refractory period has ended
> t = t + spike_refractory
Set the membrane potential at this time to the after hyperpolarization value
> V(t) = after_hyperpolarization
If the current membrane potential is less than spike_treshold:
Calculate dV according to the RC circuit model of the membrane
Set the voltage at the next timestep equal to Vinitial + dV
> V(t + dt) = Vinitial + dV
Move on to the next time step
> t = t + dt
Interactive Simulation 3: LIF model neuron¶
This simulation incorporates “spiking” into your RC circuit model using conditional statements. The stimulation protocol in the simulation is a 1 second duration (Long Square) current step.
Use this simulation to explore the relationships among membrane resistance, capacitance, equilibrium potential, spike threshold, and amplitude of an applied current.
Parameters are resported in the following units:
Applied current : pA
Membrane Capacitance : nF
Membrane Resistance : MOhm
Membrane Potentials : mV
Note, that the units of membrane capacitance and resistance have been selected to make the calculation of the Time constant (in milliseconds) easy because: Tau (in ms) = C (in nF) * R (in MOhm)
#@title {display-mode: "form"}
#@markdown Run this cell to enable the interactive plot.
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['C'] = 0.500 # membrane capacitance [nF]
pars['R_L'] = 50. # leak Resistance [MOhm]
pars['V_init'] = -75. # initial potential [mV]
pars['E_L'] = -75. # leak reversal potential [mV]
pars['tref'] = 1. # refractory time (ms)
# simulation parameters #
pars['T'] = 1200. # Total duration of simulation [ms]
pars['dt'] = .001 # Simulation time step [ms]
pars['delay'] = 100 # delay from trial onset to stimulation start [ms]
# external parameters if any #
for k in kwargs:
pars[k] = kwargs[k]
pars['range_t'] = np.arange(0, pars['T'], pars['dt'])-pars['delay'] # Vector of discretized time points [ms]
return pars
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
"""
hfig,ax = plt.subplots(figsize=(10,5),num=1)
V_th = pars['V_th']*1e-3
dt, range_t = pars['dt']*1e-3, pars['range_t']*1e-3
if sp.size:
sp_num = (sp / dt).astype(int) - 1
v[sp_num] += (40*1e-3) # draw nicer spikes
ax.plot(range_t, v, 'b')
ax.axhline(V_th, 0, pars['T'], color='k', ls='--')
ax.set_xlabel('Time (s)',fontsize=14)
ax.set_ylabel('V',fontsize=14)
# ax.set_xticks(fontsize=14)
# ax.set_yticks(fontsize=14)
ax.legend(['Membrane\npotential', r'Threshold V$_{\mathrm{th}}$'],
loc=[1.05, 0.75])
ax.set_ylim([-0.090, -0.020])
def run_LIF(pars, Iinj):
"""
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']*1e-3, pars['V_reset']*1e-3 # convert to V
R_L = pars['R_L']*1e6 #convert to Ohm from MegaOhm
C = pars['C']*1e-9 # convert to Farad from nanoF
V_init, E_L = pars['V_init']*1e-3, pars['E_L']*1e-3 # convert to V
dt, range_t = pars['dt']*1e-3, pars['range_t']*1e-3 # convert to Sec
fs = 1/dt
Lt = range_t.size
tref = pars['tref']*1e-3
# Initialize voltage
v = np.zeros(Lt)
v[0] = V_init
# Set current time course
Iinj = (Iinj*1e-12) * np.ones(Lt) # convert to Amp from picoamp
durstep = int(((1)*fs)) # hard code for now as 1 second.... was int(((durstep*1e-3)*fs)/2)
# If current pulse, set beginning and end to 0
# if stop: (for this one, just do a step
# Iinj[:int(len(Iinj) / 2) - durstep] = 0
# Iinj[int(len(Iinj) / 2) + durstep:] = 0
delay = int((pars['delay']*1e-3)*fs)
Iinj[:delay] = 0
# Iinj[:int(len(Iinj) / 2) - durstep] = 0
Iinj[(delay+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)
dv = (Iinj[it] - ((v[it] - E_L)/R_L))*(dt/C)
# 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
my_layout.width = '700px'
# my_layout.description_width = 'initial'
style = {'description_width': 'initial'}
@widgets.interact(
current_injection=widgets.FloatSlider(50., min=0., max=800., step=0.1,
layout=my_layout,style=style,continuous_update=False),
# Simulation_Duration=widgets.FloatSlider(200., min=0.02, max=5., step=0.01,
# 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=0.01, max=1000., step=0.01,
# layout=my_layout,style=style),
Equilibrium_Potential=widgets.FloatSlider(-75., min=-90, max=-30., step=0.2,
layout=my_layout,style=style,continuous_update=False),
Resistance=widgets.FloatSlider(50., min=0, max=1000., step=1.,
layout=my_layout,style=style,continuous_update=False),
# AHP_Voltage=widgets.FloatSlider(-70., min=-100., max=-30.,step=1.,
# layout=my_layout,style=style,continuous_update=False),
Spike_Threshold=widgets.FloatSlider(-55., min=-100., max=-20., step=0.2,
layout=my_layout,style=style,continuous_update=True),
Capacitance=widgets.FloatSlider(0.200, min=0.001, max=1, step=0.001,readout_format='.3f',
layout=my_layout,style=style,continuous_update=False)
)
def diff_DC(current_injection,
# Simulation_Duration,
# Injection_Step,
# Current_Step_Duration,
Equilibrium_Potential,
Resistance,
Capacitance,
# AHP_Voltage,
Spike_Threshold): #I_dc=200., tau_m=10.):
# pars = default_pars(T=100.)
# pars['tau_m'] = tau_m
pars = default_pars(
E_L = Equilibrium_Potential,
R_L = Resistance,
V_init = Equilibrium_Potential,
V_reset = Spike_Threshold - 10,
V_th = Spike_Threshold,
# tau_m = Membrane_Tau,
C = Capacitance) #,
# T = Simulation_Duration)
# pars=default_pars()
v, sp = run_LIF(pars, Iinj=current_injection) #,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)))))
Interactive Simulation 4: LIF model neuron¶
The stimulation protocol in the simulation is a 0.01-10 millisecond duration (Short Square) current step.
Use this simulation to perform a ‘strenght-duration’ stimulus threshold experiment on model neurons.
Parameters are reported in the following units:
Applied current : nA (note that this is different from the Allen database and simulation 3)
Membrane Capacitance : nF
Membrane Resistance : MOhm
Membrane Potentials : mV
Note, that the units of membrane capacitance and resistance have been selected to make the calculation of the Time constant (in milliseconds) easy because: Tau (in ms) = C (in nF) * R (in MOhm)
#@title {display-mode: "form"}
#@markdown Run this cell to enable the interactive plot.
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['C'] = 500. # membrane capacitance [pF]
pars['R_L'] = 50. # leak Resistance [MOhm]
pars['V_init'] = -75. # initial potential [mV]
pars['E_L'] = -75. # leak reversal potential [mV]
pars['tref'] = 1. # refractory time (ms)
# simulation parameters #
pars['T'] = 11. # Total duration of simulation [ms]
pars['dt'] = .001 # Simulation time step [ms]
pars['delay'] = 0.5 #delay from trial onset to stimulus onset [ms]
# external parameters if any #
for k in kwargs:
pars[k] = kwargs[k]
pars['range_t'] = np.arange(0, pars['T'], pars['dt'])-pars['delay'] # Vector of discretized time points [ms]
return pars
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
"""
hfig,ax = plt.subplots(figsize=(10,5),num=1)
V_th = pars['V_th']*1e-3
dt, range_t = pars['dt']*1e-3, pars['range_t']*1e-3
if sp.size:
sp_num = (sp / dt).astype(int) - 1
v[sp_num] += (40*1e-3) # draw nicer spikes
ax.plot(pars['range_t'], v, 'b')
ax.axhline(V_th, 0, pars['T'], color='k', ls='--')
ax.set_xlabel('Time (ms)',fontsize=14)
ax.set_ylabel('V)',fontsize=14)
# ax.set_xticks(fontsize=14)
# ax.set_yticks(fontsize=14)
ax.legend(['Membrane\npotential', r'Threshold V$_{\mathrm{th}}$'],
loc=[1.05, 0.75])
ax.set_ylim([-0.090, -0.020])
def run_LIF(pars, Iinj, durstep):
"""
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']*1e-3, pars['V_reset']*1e-3 # convert to V
R_L = pars['R_L']*1e6 #convert to Ohm
C = pars['C']*1e-9 # convert to Farad
V_init, E_L = pars['V_init']*1e-3, pars['E_L']*1e-3 # convert to V
dt, range_t = pars['dt']*1e-3, pars['range_t']*1e-3 # convert to Sec
fs = 1/dt
Lt = range_t.size
tref = pars['tref']*1e-3
# Initialize voltage
v = np.zeros(Lt)
v[0] = V_init
# Set current time course
Iinj = (Iinj*1e-9) * np.ones(Lt)
durstep = int(((durstep*1e-3)*fs))
# If current pulse, set beginning and end to 0
delay = int((0.5*1e-3)*fs)
Iinj[:delay] = 0
# Iinj[:int(len(Iinj) / 2) - durstep] = 0
Iinj[(delay+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)
dv = (Iinj[it] - ((v[it] - E_L)/R_L))*(dt/C)
# 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
my_layout.width = '700px'
# my_layout.description_width = 'initial'
style = {'description_width': 'initial'}
@widgets.interact(
current_injection=widgets.FloatSlider(0.500, min=0., max=100., step=0.05,
layout=my_layout,style=style),
# Simulation_Duration=widgets.FloatSlider(10., min=0.02, max=20., step=0.01,
# 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(10., min=0.01, max=10., step=0.01,
layout=my_layout,style=style),
Equilibrium_Potential=widgets.FloatSlider(-75., min=-90, max=-30., step=0.2,
layout=my_layout,style=style),
Resistance=widgets.FloatSlider(50., min=0, max=1000., step=1.,
layout=my_layout,style=style),
# AHP_Voltage=widgets.FloatSlider(-70., min=-100., max=-30.,step=2.,
# layout=my_layout,style=style),
Spike_Threshold=widgets.FloatSlider(-55., min=-100., max=-20., step=0.2,
layout=my_layout,style=style),
Capacitance=widgets.FloatSlider(0.200, min=0.001, max=1, step=0.001,readout_format='.3f',
layout=my_layout,style=style)
)
def diff_DC(current_injection,
# Simulation_Duration,
Current_Step_Duration,
Equilibrium_Potential,
Resistance,
Capacitance,
# AHP_Voltage,
Spike_Threshold): #I_dc=200., tau_m=10.):
# pars = default_pars(T=100.)
# pars['tau_m'] = tau_m
pars = default_pars(
E_L = Equilibrium_Potential,
R_L = Resistance,
V_init = Equilibrium_Potential,
# V_reset = AHP_Voltage,
V_th = Spike_Threshold,
# tau_m = Membrane_Tau,
C = Capacitance)#,
# T = Simulation_Duration)
# pars=default_pars()
v, sp = run_LIF(pars, Iinj=current_injection,durstep=Current_Step_Duration)
plot_volt_trace(pars, v, sp)
print('Interactive demo initiated at ' + str(datetime.now(timezone(-timedelta(hours=5)))));
Strength-Duration Tuning¶
First, enter your stimulus duration and amplitude data into a .csv file. Then, you will import that .csv file so that you can fit a hyperbolic estimate to the data. Each row in the table should be a different stimulus condition. You should have a column in the table for stimulus duration and a column for each neuron’s amplitude data that you collected. The first row of the spreadsheet needs to contain a unique identifying header/name for each column.
Read, fit, and plot strength-duration threshold data¶
The following code cell will import the data from your .csv file and plot it as a scatterplot. Additionally, the data will be fit to the following hyperbolic equation:
where \(a\) is the stimulus amplitude, \(r\) is the rheobase, \(c\) is the chronaxie, and \(t\) is the stimulus duration. The fit \(r\) and \(c\) parameter values will be printed above the figure each time you run the code cell.
#@title {display-mode:"form"}
#@markdown Speficy the filepath to a csv file
filepath = 'full filepath goes here' #@param
# filepath = '/Users/kperks/Documents/Teaching/Neurophysiology-Lab/modules/computational-model/Results-strength-duration.csv' #@param
#@markdown Speficy the header name of one column you want for your x points
x_column = 'x_header' #@param
#@markdown Speficy the header name of one column you want for your y points
y_column = 'y_header' #@param
#@markdown Then run this code cell to plot the raw data and the estimated strength-duration curve for this one sample.
#@markdown The 'r' and 'c' values will be printed out above the data plot.
df = pd.read_csv(filepath)
hfig,ax = plt.subplots(figsize=(6,5))
sns.scatterplot(x=x_column,y=y_column,data=df,color='black');
ax.set_xscale('log')
df_ = df[[x_column, y_column]].dropna()
params, covs = curve_fit(hyper_fit, df_[x_column],df_[y_column],maxfev = 5000)
t_ = np.arange(0.01,10,0.01)
a_ = hyper_fit(t_,params[0],params[1])
ax.plot(t_,a_,color = 'gray',linestyle='--')
# ax.hlines(params[0],0,10,color = 'green')
# ax.hlines(2*params[0],0,10,color = 'green',linestyle='--')
# ax.vlines(params[1],0,2*params[0],color='purple',linestyle='--')
ax.set_ylabel('amplitude (v)')
ax.set_xlabel('duration (ms)');
ax.set_ylim(-2,100)
print(f'r = {params[0]:0.2f}; c = {params[1]:0.2f}')
Compare strength-duration tuning across neurons¶
#@title {display-mode:"form"}
#@markdown Specify the rheobase for each neuron
r = [r1,r2] #@param
#@markdown Specify the chronaxie for each neuron
c = [c1,c2] #@param
#@markdown Specify a label for each sample
label = ['neuron 1 ID', 'neuron 2 ID'] #@param
#@markdown Specify whether you want to plot the x-axis as 'log' or 'linear'
x_axis_scale = 'log'
#@markdown Then run this code cell to plot the estimated strength-duration curve for each sample.
#@markdown The chronaxie will be marked with scatter dot and a gray line connected to the x-axis for each sample.
t_ = np.arange(0.01,10,0.01)
hfig,ax = plt.subplots(figsize=(6,5))
for r_,c_,label_ in zip(r,c,label):
a_ = hyper_fit(t_,r_,c_)
ax.plot(t_,a_,label=label_)
# ax.scatter(c_,2*r_,zorder=3)
# ax.vlines(c_,0,2*r_,color='gray')
ax.set_ylabel('amplitude (v)')
ax.set_xlabel('duration (ms)');
ax.set_ylim(-2,100)
ax.set_xscale(x_axis_scale)
plt.legend();
Written by Dr. Krista Perks for courses taught at Wesleyan University.