Charlotte - Data Explorer
Contents
Charlotte - Data Explorer¶
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()
# define the function
def simulate_Vin_steady_state(E1,R1,E2,R2):
Vin = ( (E1/R1) + (E2/R2) ) / ( (1/R1) + (1/R2) )
return Vin
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.
def hyper_fit(t,r,c):
return r + ((r*c)/t)
print('Task completed at ' + str(datetime.now(timezone(-timedelta(hours=5)))))
Steady State Model¶
Interactive Simulation 1: Steady-State Membrane Potential Model¶
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:
\[
\partial{V} = (Iapp - \frac{V-E}{R}) * (\frac{\partial{t}}{C})
\]
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.
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).
In shrimp (a type of aquatic crustacean like the earthworm), intracellular measurement of the medial giant fiber shows a resting membrane potential of −64.3 ±7.4 mV (mean ± S.E.M., N=7) (Xu and Terakawa, 1991). The action potential is monophasic (image from Fenestration nodes and the wide submyelinic space form the basis for the unusually fast impulse conduction of shrimp myelinated axons; originally from Xu and Terakawa, 1993).
Perhaps from this figure, we can approximate a relative spike threshold of 20mV.
Interactive Simulation 3: LIF model neuron with SHORT square stimulus¶
The stimulation protocol in the simulation is a 1 second 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'] = 5. # 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.02,
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(4., min=0.01, max=4., step=0.02,
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:
\[
a = r + \frac{(r*c)}{t}
\]
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 Specify 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.