-
Notifications
You must be signed in to change notification settings - Fork 1
/
fI_curve.py
188 lines (162 loc) · 8.08 KB
/
fI_curve.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
import os
import sys
import json
import pickle
import argparse as arg
import numpy as np
import matplotlib.pyplot as plt
from random import randint
import dlutils.utils as dlu
import dlutils.graphics as dlg
import neuron
from current_injection import inject_current_step
use_scoop = True
if use_scoop:
try:
from scoop import futures
map_fun = futures.map
except:
map_fun = map
else:
map_fun = map
if __name__ == '__main__':
parser = arg.ArgumentParser(description='Compute the f-I curve of a neuron model.')
parser.add_argument('I', type=str, action='store', help='current values in pA, either comma separated or interval and steps, as in 100:300:50')
parser.add_argument('-f','--swc-file', type=str, help='SWC file defining the cell morphology', required=True)
parser.add_argument('-p','--params-files', type=str, default='', help='JSON file(s) containing the parameters of the model (comma separated)')
parser.add_argument('-m','--mech-file', type=str, default='', help='JSON file containing the mechanisms to be inserted into the cell')
parser.add_argument('-c','--config-file', type=str, default='', help='JSON file(s) containing the configuration')
parser.add_argument('-n','--cell-name', default='', type=str, help='cell name, if the mechanisms are stored in new style format')
parser.add_argument('-P','--pickle-file', type=str, default='', help='Pickle file containing the parameters of a population of individuals')
parser.add_argument('-e','--evaluator-file', type=str, default='evaluator.pkl', help='Pickle file containing the evaluator')
parser.add_argument('-o','--output', type=str, default='fI_curve.pkl', help='Output file name (default: fI_curve.pkl)')
parser.add_argument('-R','--replace-axon', type=str, default=None, help='whether to replace the axon (accepted values: "yes" or "no")')
parser.add_argument('-A', '--add-axon-if-missing', type=str, default=None, \
help='whether add an axon if the cell does not have one (accepted values: "yes" or "no")')
parser.add_argument('--delay', default=1000., type=float, help='delay before stimulation onset (default: 1000 ms)')
parser.add_argument('--dur', default=500., type=float, help='stimulation duration (default: 500 ms)')
parser.add_argument('--tran', default=50., type=float, help='transient to be discard after stimulation onset (default: 50 ms)')
args = parser.parse_args(args=sys.argv[1:])
try:
I = np.array([float(args.I)])
except:
if ',' in args.I:
I = np.sort([float(x) for x in args.I.split(',')])
elif ':' in args.I:
tmp = [float(x) for x in args.I.split(':')]
I = np.arange(tmp[0],tmp[1]+tmp[2]/2,tmp[2])
else:
print('Unknown current definition: %s.' % args.I)
sys.exit(1)
new_config_style = False
if args.config_file != '':
new_config_style = True
if new_config_style and args.cell_name == '':
print('You must provide the --cell-name option along with the --config-file option.')
sys.exit(1)
if '*' in args.params_files:
import glob
params_files = glob.glob(args.params_files)
else:
params_files = args.params_files.split(',')
if args.mech_file == '':
if not new_config_style:
print('You must provide the --mech-file option if no configuration file is specified.')
sys.exit(1)
cell_name = args.cell_name
mechanisms = dlu.extract_mechanisms(args.config_file, cell_name)
else:
cell_name = None
mechanisms = json.load(open(args.mech_file,'r'))
if args.pickle_file == '':
population = [json.load(open(params_file,'r')) for params_file in params_files]
working_dir = os.path.split(params_files[0])[0]
else:
if len(params_files) > 1:
print('You cannot specify multiple parameter files and one pickle file.')
sys.exit(1)
population = dlu.individuals_from_pickle(args.pickle_file, args.config_file, cell_name, args.evaluator_file)
working_dir = os.path.split(args.pickle_file)[0]
if working_dir == '':
working_dir = '.'
try:
sim_pars = pickle.load(open(working_dir + '/simulation_parameters.pkl','rb'))
if working_dir == '.':
print('Found pickle file with simulation parameters in current directory.')
else:
print('Found pickle file with simulation parameters in {}.'.format(working_dir))
except:
sim_pars = None
if args.replace_axon is None and args.add_axon_if_missing is None:
print('No pickle file with simulation parameters in {}.'.format(working_dir))
if args.replace_axon is None:
if sim_pars is None:
replace_axon = False
else:
replace_axon = sim_pars['replace_axon']
print('Setting replace_axon = {} as per original optimization.'.format(replace_axon))
else:
if args.replace_axon.lower() in ('y','yes'):
replace_axon = True
elif args.replace_axon.lower() in ('n','no'):
replace_axon = False
else:
print('Unknown value for --replace-axon: "{}".'.format(args.replace_axon))
sys.exit(3)
if args.add_axon_if_missing is None:
if sim_pars is None:
add_axon_if_missing = True
else:
add_axon_if_missing = not sim_pars['no_add_axon']
print('Setting add_axon_if_missing = {} as per original optimization.'.format(add_axon_if_missing))
else:
if args.add_axon_if_missing.lower() in ('y','yes'):
add_axon_if_missing = True
elif args.add_axon_if_missing.lower() in ('n','no'):
add_axon_if_missing = False
else:
print('Unknown value for --add-axon-if-missing: "{}".'.format(args.add_axon_if_missing))
sys.exit(4)
dur = args.dur
delay = args.delay
tran = args.tran
N = len(population)
f = np.zeros((N,len(I)))
no_spikes = np.zeros((N,len(I)))
inverse_first_isi = np.zeros((N,len(I)))
inverse_last_isi = np.zeros((N,len(I)))
spike_times = []
inj_loc = 'soma'
inj_dist = 0
for i,individual in enumerate(population):
worker = lambda Idc: inject_current_step(Idc, delay, dur, args.swc_file, inj_loc, inj_dist, individual, mechanisms,
replace_axon=replace_axon, add_axon_if_missing=add_axon_if_missing,
cell_name=None, neuron=neuron, do_plot=False, verbose=False)
curve = list(map_fun(worker, I))
neuron.h('forall delete_section()')
spks = [np.array(point['spike_times']) for point in curve]
no_spikes[i,:] = [len(x)/dur*1e3 for x in spks]
f[i,:] = [len(x[(x>delay+tran) & (x<delay+dur)])/(dur-tran)*1e3 for x in spks]
inverse_first_isi[i,:] = [1e3/np.diff(t[:2]) if len(t) > 1 else 0 for t in spks]
inverse_last_isi[i,:] = [1e3/np.diff(t[-2:]) if len(t) > 1 else 0 for t in spks]
spike_times.append(spks)
data = {'delay': delay, 'dur': dur, 'tran': tran,
'I': I, 'spike_times': spike_times,
'f': f, 'no_spikes': no_spikes,
'inverse_first_isi': inverse_first_isi,
'inverse_last_isi': inverse_last_isi,
'population': population}
pickle.dump(data, open(args.output,'wb'))
fig,ax = plt.subplots(1,1)
dlg.plot_means_with_errorbars(I*1e-3,no_spikes,mode='sem',ax=ax,color='r',label='All spikes')
dlg.plot_means_with_errorbars(I*1e-3,f,mode='sem',ax=ax,color='k',label='With transient removed')
dlg.plot_means_with_errorbars(I*1e-3,inverse_first_isi,mode='sem',ax=ax,color='b',label='Inverse first ISI')
#dlg.plot_means_with_errorbars(I*1e-3,inverse_last_isi,mode='sem',ax=ax,color='m',label='Inverse last ISI')
plt.xlabel('Current (nA)')
plt.ylabel(r'$f$ (spikes/s)')
plt.legend(loc='best')
folder = os.path.dirname(args.output)
if folder == '':
folder = '.'
plt.savefig(folder + '/' + os.path.basename(args.output).split('.')[0] + '.pdf')
plt.show()