"""
(C) Asaph Zylbertal 01.03.2015, HUJI, Jerusalem, Israel
Run the mitral cell used in the article and produce example figures (Fig 4, Fig 5A-C, Fig 6C)
"""
import neuron
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import filtfilt
def draw_model (mitral_mod):
inst_dat={}
inst_dat['v']=[]
inst_dat['t']=[]
# run model until steady state is reached, save in mitral_mod.steady
mitral_mod.init_steady_state(mitral_mod.soma(0.5), min_slope=5e-9, init_run_chunk=2000.)
# run model for short current injections
for amp in [-0.06, -0.03, 0.03, 0.06, 0.1, 0.15, 0.2, 0.35]:
mitral_mod.steady.restore()
mitral_mod.init_square_stim(amp)
mitral_mod.init_recording(mitral_mod.soma(0.5))
mitral_mod.run_model()
outv=np.array(mitral_mod.rec_v)
outt=np.array(mitral_mod.rec_t)
mitral_mod.stop_recording()
inst_dat['v'].append(outv)
inst_dat['t'].append(outt)
del mitral_mod.steady
del mitral_mod.stim
######### short current injections figure (article Fig 4) ##########
fig=plt.figure(figsize=(5.5, 12))
pas=fig.add_subplot(7, 1, 1)
pas.plot(inst_dat['t'][0]/1000., inst_dat['v'][0],'g')
pas.plot(inst_dat['t'][1]/1000., inst_dat['v'][1],'g')
plt.tick_params(axis='x', which='both', bottom='on', labelbottom='off')
plt.ylim(-75., -30.)
plt.title('Passive response')
plt.ylabel('Vm (mV)')
mod2=fig.add_subplot(7, 1, 2)
mod2.plot(inst_dat['t'][2], inst_dat['v'][2],'g')
plt.ylim(-75., 60.)
plt.title('I=30pA')
plt.ylabel('Vm (mV)')
plt.tick_params(axis='x', which='both', bottom='on', labelbottom='off')
mod3=fig.add_subplot(7, 1, 3)
mod3.plot(inst_dat['t'][3], inst_dat['v'][3],'g')
plt.ylim(-75., 60.)
plt.title('I=60pA')
plt.ylabel('Vm (mV)')
plt.tick_params(axis='x', which='both', bottom='on', labelbottom='off')
mod4=fig.add_subplot(7, 1, 4)
mod4.plot(inst_dat['t'][4], inst_dat['v'][4],'g')
plt.ylim(-75., 60.)
plt.title('I=100pA')
plt.ylabel('Vm (mV)')
plt.tick_params(axis='x', which='both', bottom='on', labelbottom='off')
mod5=fig.add_subplot(7, 1, 5)
mod5.plot(inst_dat['t'][5], inst_dat['v'][5],'g')
plt.ylim(-75., 60.)
plt.title('I=150pA')
plt.ylabel('Vm (mV)')
plt.tick_params(axis='x', which='both', bottom='on', labelbottom='off')
mod6=fig.add_subplot(7, 1, 6)
mod6.plot(inst_dat['t'][6], inst_dat['v'][6],'g')
plt.ylim(-75., 60.)
plt.title('I=200pA')
plt.ylabel('Vm (mV)')
plt.tick_params(axis='x', which='both', bottom='on', labelbottom='off')
mod7=fig.add_subplot(7, 1, 7)
mod7.plot(inst_dat['t'][7], inst_dat['v'][7],'g')
plt.ylim(-75., 60.)
plt.title('I=350pA')
plt.xlabel('t (s)')
plt.ylabel('Vm (mV)')
################
four_spikes_frame_time=20.
trains_frame_time=80.
four_spikes_sim_time=15000
trains_sim_time=60000
# record "fluorescence" from tuft #1
recording_comp=mitral_mod.tuft1(0.5)
inst_dat={}
inst_dat['f']=[]
inst_dat['v']=[]
inst_dat['t']=[]
inst_dat['i']=[]
init_run_time=1500000
######### pump and exchanger current figure (article figure 6C) #########
plt.figure()
ca_pump_current=mitral_mod.init_vec_recording(mitral_mod.tuft1(0.5)._ref_ica_pmp_cadp)
ca_ncx_rate=mitral_mod.init_vec_recording(mitral_mod.tuft1(0.5)._ref_rate_ncx)
for freq in [1., 15., 30.]:
if freq==1.:
# when the frequency is 1Hz (only four spikes) clamp to -55mV and inject -0.03nA during IC epoch (like the real experiment)
(t, ih, fl, v)=hybrid_clamp(mitral_mod,init_run_time, 3800-140,4400,4000,4000,freq/1000, 1,7.,four_spikes_sim_time,-55,6,-0.03, return_fluor=True, fluor_comp=recording_comp)
else:
# when the frequency is 15Hz or 30Hz clamp to -70mV and inject -0.06nA during IC epoch (like the real experiment)
(t, ih, fl, v)=hybrid_clamp(mitral_mod,init_run_time, 3800-140,4400,4000,4000,freq/1000, 1,7.,trains_sim_time,-70,6,-0.06, return_fluor=True, fluor_comp=recording_comp)
outt=np.array(t)
outf=np.array(fl)
outv=np.array(v)
if not freq==1:
outi=np.array(ih)
inst_dat['i'].append(outi)
inst_dat['t'].append(outt)
inst_dat['f'].append(outf)
inst_dat['v'].append(outv)
measure_frame_time=trains_frame_time/100
if freq>1:
interp_block_ca_pump_current=np.interp(np.arange(0,trains_sim_time+init_run_time,measure_frame_time), outt, np.array(ca_pump_current))[(init_run_time/measure_frame_time):]
interp_block_ca_ncx_current=-2*np.interp(np.arange(0,trains_sim_time+init_run_time,measure_frame_time), outt, np.array(ca_ncx_rate))[(init_run_time/measure_frame_time):]
measures_t=np.arange(0, trains_sim_time, measure_frame_time)
if freq==30:
plt.plot(measures_t*1e-3 , interp_block_ca_pump_current*1e3, 'r')
plt.plot(measures_t*1e-3 , interp_block_ca_ncx_current*1e3, 'g')
plt.plot(measures_t*1e-3 , (interp_block_ca_ncx_current+interp_block_ca_pump_current)*1e3, 'b')
plt.plot([0, 50], [0, 0], '--k')
plt.xlim(0, 50)
plt.ylim(-1., 1.)
plt.title('Ca2+ pump and exchanger currents')
plt.legend(('Ca2+ pump', 'Na+ - Ca2+ exchanger', 'Total'))
plt.xlabel('t (sec)')
plt.ylabel(r'$\mu A/cm^2$')
########## four spikes fluorescence figure (article figure 5A) #############
plt.figure()
interp_block=interp_f_result(inst_dat, four_spikes_sim_time, init_run_time, four_spikes_frame_time, 0, 0, mitral_mod.params['filt_order'], mitral_mod.params['time_shift'])
plt.plot(np.arange(0,four_spikes_sim_time,four_spikes_frame_time)*1e-3, interp_block, 'b', linewidth=2.0)
plt.xlabel('t (sec)')
plt.ylabel('df/f')
plt.title('Tuft fluorescence - 1Hz spiking')
########### 15Hz and 30Hz fluorescence and current figure (article figure 5B-C) ##########
fig2=plt.figure(figsize=(14, 5))
subfig=[0]*2
subfig[0]=fig2.add_subplot(1, 2, 1)
interp_block=interp_f_result(inst_dat, trains_sim_time, init_run_time, trains_frame_time, 1, 1, mitral_mod.params['filt_order'], mitral_mod.params['time_shift'])
subfig[0].plot(np.arange(0, trains_sim_time, trains_frame_time)/1000., filtfilt(np.ones(5)/5, [1], interp_block), 'g', linewidth=2.0)
interp_block=interp_f_result(inst_dat, trains_sim_time, init_run_time, trains_frame_time, 2, 2, mitral_mod.params['filt_order'], mitral_mod.params['time_shift'])
subfig[0].plot(np.arange(0, trains_sim_time, trains_frame_time)/1000., filtfilt(np.ones(5)/5, [1], interp_block), 'r', linewidth=2.0)
subfig[0].legend(('15Hz', '30Hz'))
plt.xlabel('t (sec)')
plt.ylabel('df/f')
plt.title('Tuft fluorescence')
subfig[1]=fig2.add_subplot(1, 2, 2)
interp_block= interp_i_result(inst_dat, trains_sim_time, init_run_time, trains_frame_time, 1, 0)
subfig[1].plot(np.arange(0, trains_sim_time, trains_frame_time)/1000., interp_block, 'g', linewidth=2.0)
interp_block=interp_i_result(inst_dat, trains_sim_time, init_run_time, trains_frame_time, 2, 1)
subfig[1].plot(np.arange(0, trains_sim_time, trains_frame_time)/1000., interp_block, 'r', linewidth=2.0)
subfig[1].legend(('15Hz', '30Hz'))
plt.xlabel('t (sec)')
plt.ylabel('I (nC)')
plt.title('Current recorded in the soma')
plt.show()
def hybrid_clamp(inst, init_run_time, ic_delay, ic_duration, train_delay, train_duration, freq, amp, pulse_duration, sim_time, v_clamp, rs, dc=0.0, return_fluor=False, clamp_vec=None, clamp_t=None, fluor_comp=None):
"""
run hybrid clamp experiment in the model:
----------------------------------------------
inst - model cell instance
init_run_time - initialization run time (how long to run to achieve steady state values in all state variables)
ic_delay - delay before switch to from VC to IC
ic_duration - duration of IC epoch
train_delay - delay before pulse train injection
train_duration - duration of the injected pulse train
freq - pulse frequency
amp - pulse amplitude
pulse_duration - duration of each pulse in the train
sim_tim - simulation duration (after initialization run)
v_clamp - voltage to clamp to during VC
rs - series resistance
dc - DC current injection during IC epoch
return_fluor (boolean) - should the function return simulated fluorescence data?
clamp_vec - vector of command voltage (if not constant)
clamp_t - time stamps for clamp_vec
fluo_comp - compartment to read simulated fluorescence from (default=soma)
"""
if fluor_comp==None:
fluor_comp=inst.soma(0.5)
t=inst.init_vec_recording(neuron.h._ref_t)
hc=neuron.h.hybrid(inst.soma(0.5))
i_hybrid=inst.init_vec_recording(hc._ref_i)
v=inst.init_vec_recording(inst.tuft1(0.5)._ref_v)
if return_fluor:
fl=inst.init_vec_recording(fluor_comp.cadp._ref_f)
inst.init_train_stim(train_delay+init_run_time, train_duration, freq, pulse_duration, amp, 0.0, limit_dc=True)
dcstim=neuron.h.IClamp(inst.soma(0.5))
dcstim.delay=ic_delay+init_run_time
dcstim.dur=ic_duration
dcstim.amp=dc
hc.rs=rs
hc.delay=ic_delay+init_run_time
hc.dur=ic_duration
hc.tot_dur=sim_time+init_run_time
if clamp_vec==None:
hc.vc_amp=v_clamp
else:
stimv_vec=neuron.h.Vector(clamp_vec)
t_vec=neuron.h.Vector(clamp_t)
stimv_vec.play(hc._ref_vc_amp, t_vec)
if inst.cv.active()==1:
inst.cv.re_init()
neuron.h.finitialize(v_clamp)
neuron.h.fcurrent()
neuron.run(sim_time+init_run_time)
del inst.stim
del hc
del dcstim
if return_fluor:
return (t, i_hybrid, fl, v)
else:
return (t, i_hybrid,v)
def interp_f_result(inst_dat, sim_time, init_time, frame_time, t_num, f_num, filt_order, time_shift):
"""
Interpolate and process fluorescence result
"""
interp_block=np.interp(np.arange(0,sim_time+init_time,frame_time), inst_dat['t'][t_num], inst_dat['f'][f_num])[(init_time/frame_time):]
interp_block=df_f(interp_block)
interp_block=filtfilt(np.ones(filt_order)/filt_order, [1], interp_block)
origin_len=len(interp_block)
interp_block=np.concatenate([np.zeros(time_shift/frame_time), interp_block])
interp_block=interp_block[0:origin_len]
return interp_block
def interp_i_result(inst_dat, sim_time, init_time, frame_time, t_num, i_num):
"""
Interpolate current result
"""
interp_block=np.interp(np.arange(0,sim_time+init_time,frame_time), inst_dat['t'][t_num], inst_dat['i'][i_num])[(init_time/frame_time):]
return interp_block-np.mean(interp_block[0:3000/frame_time])
def df_f(f_vec):
f_rest=np.mean(f_vec[0:10])
return (f_vec-f_rest)/f_rest