Skip to content
Snippets Groups Projects
Commit aa84b112 authored by Doorn, Nina (UT-TNW)'s avatar Doorn, Nina (UT-TNW)
Browse files

Fixed MEA feature names

parent 1c4e4ebb
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 150 additions and 149 deletions
Simulator.py
+ 150
− 149
View file @ aa84b112
......@@ -10,68 +10,69 @@ from scipy.stats import norm
from itertools import combinations
import numpy as np
def MEAnetSimulate(parameter_set, simtime=165*second, transient=5*second):
# network parameters
def MEAnetsimulate(parameter_set, simtime=165*second, transient=5*second):
# Network parameters
Nl = 10
N = Nl * Nl # number of neurons
# neuron parameters
area = 300 * umetre ** 2
Cm = (1 * ufarad * cm ** -2) * area # membrane capacitance
El = -39.2 * mV # Nernst potential of leaky ions
EK = -80 * mV # Nernst potential of potassium
ENa = 70 * mV # Nernst potential of sodium
gam_na = parameter_set[1].item()
g_na = gam_na * (50 * msiemens * cm ** -2) * area # maximal conductance of sodium channels
gam_kd = parameter_set[2].item()
g_kd = gam_kd * (5 * msiemens * cm ** -2) * area # maximal conductance of potassium
gl = (0.3 * msiemens * cm ** -2) * area # maximal leak conductance
VT = -30.4 * mV # alters firing threshold of neurons
sigma = parameter_set[0].item() * mV # standard deviation of the noisy voltage fluctuations
N = Nl * Nl # number of neurons
# Neuron parameters
area = 300 * umetre ** 2
Cm = (1 * ufarad * cm ** -2) * area # membrane capacitance
E_l = -39.2 * mV # nernst potential of leaky ions
E_K = -80 * mV # nernst potential of potassium
E_Na = 70 * mV # nernst potential of sodium
gam_na = parameter_set[1].item() # modifies sodium channel conductance
g_na = gam_na * (50 * msiemens * cm ** -2) * area # maximal conductance of sodium channels
gam_kd = parameter_set[2].item() # modifies potassium channel conductance
g_kd = gam_kd * (5 * msiemens * cm ** -2) * area # maximal conductance of delayed rectifyer potassium channels
g_l = (0.3 * msiemens * cm ** -2) * area # maximal leak conductance
V_T = -30.4 * mV # adapts firing threshold of neurons
sigma = parameter_set[0].item() * mV # standard deviation of the noisy voltage fluctuations
# Adaptation parameters
E_AHP = EK # Nernst potential of afterhyperpolarization current
g_AHP = parameter_set[3].item() * nS # maximal conductance of sAHP channels
tau_Ca = 6000 * ms # recovery time constant AHP channels
alpha_Ca = 0.00035 # strength of the spike-frequency adaptation --DS MODEL VALUE: 0.00005
E_AHP = E_K # nernst potential of afterhyperpolarization potassium current
g_AHP = parameter_set[3].item() * nS # maximal conductance of sAHP potassium channels
tau_Ca = 6000 * ms # recovery time constant sAHP channels
alpha_Ca = 0.00035 # strength of the spike-frequency adaptation
# synapse parameters
g_ampa = parameter_set[4].item() * nS # maximal conductance of AMPA channels
g_nmda = parameter_set[5].item() * nS
E_ampa = 0 * mV # Nernst potentials of synaptic channels
g_ampa = parameter_set[4].item() * nS # maximal conductance of AMPA channels
g_nmda = parameter_set[5].item() * nS # maximal conductance of NMDA channels
E_ampa = 0 * mV # nernst potentials of synaptic channels
E_nmda = 0 * mV
tau_ampa = 2 * ms # recovery time constant of ampa conductance
taus_nmda = 100 * ms # decay time constant of nmda conductance
taux_nmda = 2 * ms # rise time constant of nmda conductance
tau_ampa = 2 * ms # recovery time constant of ampa conductance
taus_nmda = 100 * ms # decay time constant of nmda conductance
taux_nmda = 2 * ms # rise time constant of nmda conductance
alpha_nmda = 0.5 * kHz
prob = parameter_set[6].item() # connection probability
sd = 0.7 # SD of normal distribution of synaptic weights
Maxdelay = 25 * ms # maximum conduction delay between neurons furthest away
prob = parameter_set[6].item() # connection probability
sd = 0.7 # SD of normal distribution of synaptic weights
max_delay = 25 * ms # maximum conduction delay between neurons furthest away
tau_d = parameter_set[7].item() * ms # Recovery time constant of synaptic depression
U = parameter_set[8].item() # Magnitude of STD
tau_d = parameter_set[7].item() * ms # Recovery time constant of synaptic depression (STD)
U = parameter_set[8].item() # strength of STD
tau_ar = 700 * ms # Time constant of asynchronous release
Uar = parameter_set[9].item() # Strength of asynchronous release
Umax = 0.5 / ms # Saturation level of asynchronous release
x0 = 5 # Number of neurotransmitters in a vesicle
tau_ar = 700 * ms # time constant of asynchronous release
U_ar = parameter_set[9].item() # strength of asynchronous release
U_max = 0.5 / ms # saturation level of asynchronous release
x0 = 5 # number of neurotransmitters in a vesicle
# BUILD NETWORK
# neuron model
eqs = Equations('''
dV/dt = noise + (-gl*(V-El)-g_na*(m*m*m)*h*(V-ENa)-g_kd*(n*n*n*n)*(V-EK)+I-I_syn+I_AHP)/Cm : volt
dV/dt = noise + (-g_l*(V-E_l)-g_na*(m*m*m)*h*(V-E_Na)-g_kd*(n*n*n*n)*(V-E_K)+I-I_syn+I_AHP)/Cm : volt
dm/dt = alpha_m*(1-m)-beta_m*m : 1
dh/dt = alpha_h*(1-h)-beta_h*h : 1
dn/dt = (alpha_n*(1-n)-beta_n*n) : 1
dhp/dt = 0.128*exp((17.*mV-V+VT)/(18.*mV))/ms*(1.-hp)-4./(1+exp((30.*mV-V+VT)/(5.*mV)))/ms*h : 1
alpha_m = 0.32*(mV**-1)*4*mV/exprel((13*mV-V+VT)/(4*mV))/ms : Hz
beta_m = 0.28*(mV**-1)*5*mV/exprel((V-VT-40*mV)/(5*mV))/ms : Hz
alpha_h = 0.128*exp((17*mV-V+VT)/(18*mV))/ms : Hz
beta_h = 4./(1+exp((40*mV-V+VT)/(5*mV)))/ms : Hz
alpha_n = 0.032*(mV**-1)*5*mV/exprel((15*mV-V+VT)/(5*mV))/ms : Hz
beta_n = .5*exp((10*mV-V+VT)/(40*mV))/ms : Hz
noise = sigma*(2*gl/Cm)**.5*randn()/sqrt(dt) :volt/second (constant over dt)
dhp/dt = 0.128*exp((17.*mV-V+V_T)/(18.*mV))/ms*(1.-hp)-4./(1+exp((30.*mV-V+V_T)/(5.*mV)))/ms*h : 1
alpha_m = 0.32*(mV**-1)*4*mV/exprel((13*mV-V+V_T)/(4*mV))/ms : Hz
beta_m = 0.28*(mV**-1)*5*mV/exprel((V-V_T-40*mV)/(5*mV))/ms : Hz
alpha_h = 0.128*exp((17*mV-V+V_T)/(18*mV))/ms : Hz
beta_h = 4./(1+exp((40*mV-V+V_T)/(5*mV)))/ms : Hz
alpha_n = 0.032*(mV**-1)*5*mV/exprel((15*mV-V+V_T)/(5*mV))/ms : Hz
beta_n = .5*exp((10*mV-V+V_T)/(40*mV))/ms : Hz
noise = sigma*(2*g_l/Cm)**.5*randn()/sqrt(dt) :volt/second (constant over dt)
I : amp
I_syn = I_ampa + I_nmda: amp
I_ampa = g_ampa*(V-E_ampa)*(s_ampa) : amp
......@@ -90,8 +91,8 @@ def MEAnetSimulate(parameter_set, simtime=165*second, transient=5*second):
method='exponential_euler')
# Initialize neuron parameters
P.V = -39 * mV # approximately resting membrane potential
P.I = '(rand() -0.5)* 19 * pA' # Make neurons heterogeneously excitable
P.V = -39 * mV # approximately resting membrane potential
P.I = '(rand() -0.5)* 19 * pA' # make neurons heterogeneously excitable
# Position neurons on a grid
grid_dist = 45 * umeter
......@@ -112,7 +113,7 @@ def MEAnetSimulate(parameter_set, simtime=165*second, transient=5*second):
eqs_onpre = '''
x_nmda += 1
x_d *= (1-U)
uar += Uar*(Umax-uar)
uar += U_ar*(U_max-uar)
s_ampa += w * x_d
'''
# Make synapses
......@@ -120,57 +121,56 @@ def MEAnetSimulate(parameter_set, simtime=165*second, transient=5*second):
Conn.connect(p=prob)
Conn.w[:] = 'clip(1.+sd*randn(), 0, 2)'
Conn.x_d[:] = 1
Vmax = (sqrt(Nl ** 2 + Nl ** 2) * grid_dist) / Maxdelay
Conn.delay = '(sqrt((x_pre - x_post)**2 + (y_pre - y_post)**2))/Vmax'
V_max = (sqrt(Nl ** 2 + Nl ** 2) * grid_dist) / max_delay
Conn.delay = '(sqrt((x_pre - x_post)**2 + (y_pre - y_post)**2))/V_max'
# SET UP MONITORS AND RUN
recordstring = ['V']
dt2 = defaultclock.dt # Allows for chanching the timestep of recording
record_string = ['V']
trace = StateMonitor(P, recordstring, record=True, dt=dt2)
dt2 = defaultclock.dt # allows for changing the timestep of recording
trace = StateMonitor(P, record_string, record=True, dt=dt2)
spikes = SpikeMonitor(P)
try:
run(simtime, report='text', profile=True) # run
run(simtime, report='text', profile=True) # run
except:
simulation_num = int(float(sys.argv[1]))
torch.save([nan, nan, nan, nan, nan, nan, nan, nan], 'Results' + str(simulation_num) + '.pt')
print("Simulation failed, results are set to nan and python will be terminated")
exit()
# set up a filter to filter the voltage signal
fs = 1 / (dt2 / second)
fc = 100 # Cut-off frequency of the filter
w = fc / (fs / 2) # Normalize the frequency
# Set up a filter to filter the voltage signal
fs = 1 / (dt2 / second) # nyquist frequency
fc = 100 # cut-off frequency of the high-pass filter
w = fc / (fs / 2) # normalize the frequency
b, a = signal.butter(5, w, 'high')
voltagetraces = zeros((12, len(trace.t)))
elec_grid_dist = (grid_dist * (Nl - 1)) / 4 # electrode grid size (there are 12 electrodes)
elec_range = 3 * grid_dist # measurement range of each electrode
elec_grid_dist = (grid_dist * (Nl - 1)) / 4 # electrode grid size (there are 12 electrodes)
elec_range = 3 * grid_dist # measurement range of each electrode
comp_dist = ((Nl - 1) * grid_dist - elec_grid_dist * 3) / 2
elecranges = pickle.load(open("elecranges.dat", "rb"))
elecranges = pickle.load(open("elecranges.dat", "rb"))
k = 0
MaxAPs = len(spikes.t)
max_APs = len(spikes.t) # maximum amount of detected APs
APs = np.array([0, 0])
for i in [1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14]:
templist = elecranges[i, :]
x_electrode = i % 4 * elec_grid_dist + comp_dist
y_electrode = i // 4 * elec_grid_dist + comp_dist
templist = [j for j in templist if j != 0]
Voltage = np.zeros(len(trace.V[0]))
voltage = np.zeros(len(trace.V[0]))
for l in range(len(templist)):
Voltage += trace[templist[l]].V / mV * 1 / (
sqrt((x_electrode - P.x[templist[l]]) ** 2 + (y_electrode - P.y[templist[l]]) ** 2) / (
grid_dist * 0.2))
Voltage = Voltage - mean(Voltage)
Voltagefilt = signal.filtfilt(b, a, Voltage) # high pass filter
threshold = 4 * np.sqrt(np.mean(Voltagefilt ** 2)) # threshold to detect APs
APstemp, _ = signal.find_peaks(abs(Voltagefilt), height=threshold)
for j in range(len(APstemp)):
APs = np.append(APs, [k, APstemp[j]])
voltage += trace[templist[l]].V / mV * 1 / (
sqrt((x_electrode - P.x[templist[l]]) ** 2 + (y_electrode - P.y[templist[l]]) ** 2) / (grid_dist * 0.2))
voltage = voltage - mean(voltage)
voltagefilt = signal.filtfilt(b, a, voltage) # high pass filter
threshold = 4 * np.sqrt(np.mean(voltagefilt ** 2)) # threshold to detect APs
APs_temp, _ = signal.find_peaks(abs(voltagefilt), height=threshold)
for j in range(len(APs_temp)):
APs = np.append(APs, [k, APs_temp[j]])
templist = None
k += 1
......@@ -179,32 +179,32 @@ def MEAnetSimulate(parameter_set, simtime=165*second, transient=5*second):
return APs, simtime, transient, fs
def ComputeFeatures(APs, simtime, transient, fs):
numelectrodes = 12 # number of electrodes
timeBin = 25 * ms # timebin to compute network firing rate
def compute_features(APs, simtime, transient, fs):
numelectrodes = 12 # number of electrodes
timbin = 25 * ms # timebin to compute network firing rate
# initialize
APsinbin = zeros((numelectrodes, int(floor((simtime - transient) / timeBin))))
spikerate = zeros((numelectrodes, int(floor((simtime - transient) / timeBin))))
# Initialize
APs_inbin = zeros((numelectrodes, int(floor((simtime - transient) / timbin))))
spikerate = zeros((numelectrodes, int(floor((simtime - transient) / timbin))))
# delete transient
# Delete transient
APs_wot = APs[APs[:, 1] > transient * fs / second, :]
APs_wot[:, 1] = APs_wot[:, 1] - transient * fs / second
# calculate the network firing rate
# Calculate the network firing rate
for l in range(numelectrodes):
APt = APs_wot[APs_wot[:, 0] == l, :] # go through one electrode first
APt = APs_wot[APs_wot[:, 0] == l, :] # go through one electrode first
APt = APt[:, 1]
binsize = timeBin * fs / second
for k in range(int(floor((simtime - transient) / timeBin))):
APsinbin[l, k] = sum((APt > (k * binsize)) & (APt < ((k + 1) * binsize)))
spikerate[l, k] = APsinbin[l, k] * second / timeBin
binsize = timbin * fs / second
for k in range(int(floor((simtime - transient) / timbin))):
APs_inbin[l, k] = sum((APt > (k * binsize)) & (APt < ((k + 1) * binsize)))
spikerate[l, k] = APs_inbin[l, k] * second / timbin
APsinbintot = sum(APsinbin, axis=0)
APs_inbintot = sum(APs_inbin, axis=0)
spikeratetot = sum(spikerate, axis=0)
# make smoother spikerate
width = 11
# Smoothen the spikerate by convolution with gaussian kernel
width = 11
sigma = 3.0
x = np.arange(0, width, 1, float)
x = x - width // 2
......@@ -212,93 +212,93 @@ def ComputeFeatures(APs, simtime, transient, fs):
kernel /= np.sum(kernel)
spikeratesmooth = np.convolve(spikeratetot, kernel, mode='same')
# detect fragmented bursts on smoothed spikerate
MBth = (1/16) * max(spikeratetot)
peaks, ph = find_peaks(spikeratesmooth, height=MBth, prominence=(1 / 20) * max(spikeratetot))
# Detect fragmented bursts on smoothed spikerate
MB_th = (1/16) * max(spikeratetot)
peaks, ph = find_peaks(spikeratesmooth, height=MB_th, prominence=(1 / 20) * max(spikeratetot))
# set parameters for burst detection
actelec = sum(mean(spikerate, axis=1) > 0.02) # calculate the number of active electrodes
# Set parameters for burst detection
act_elec = sum(mean(spikerate, axis=1) > 0.02) # calculate the number of active electrodes
Startth = 0.25 * max(spikeratetot) # threshold to start a burst
Tth = int((50 * ms) / timeBin) # how long it has to surpass threshold
Eth = 0.5 * actelec # active electrode number threshold
Stopth = (1 / 50) * max(spikeratetot) # threshold to stop a burst
start_th = 0.25 * max(spikeratetot) # spikerate threshold to start a burst
t_th = int((50 * ms) / timbin) # how long it has to surpass threshold for to start burst
e_th = 0.5 * act_elec # how many electrodes need to be active in the burst
stop_th = (1 / 50) * max(spikeratetot) # threshold to end a burst
# initialize
# Initialize burst detection
i = 0
NBcount = 0
maxNB = 1000
NBs = zeros((maxNB, 4))
# detect NBs
while (i + Tth) < len(spikeratetot):
if (all(spikeratetot[i:i + Tth] > Startth)) \
& (sum(sum(APsinbin[:, i:i + Tth], axis=1) > Tth) > Eth):
NBs[NBcount, 2] = NBs[NBcount, 2] + sum(APsinbintot[i:i + Tth])
NBs[NBcount, 0] = i
i = i + Tth
while any(spikeratetot[i:i + 2 * Tth] > Stopth):
NBs[NBcount, 2] = NBs[NBcount, 2] + APsinbintot[i]
NB_count = 0
max_NBs = 1000 # maximum amount of to be detected bursts
NBs = zeros((max_NBs, 4))
# Detect NBs
while (i + t_th) < len(spikeratetot):
if (all(spikeratetot[i:i + t_th] > start_th)) \
& (sum(sum(APs_inbin[:, i:i + t_th], axis=1) > t_th) > e_th):
NBs[NB_count, 2] = NBs[NB_count, 2] + sum(APs_inbintot[i:i + t_th])
NBs[NB_count, 0] = i
i = i + t_th
while any(spikeratetot[i:i + 2 * t_th] > stop_th):
NBs[NB_count, 2] = NBs[NB_count, 2] + APs_inbintot[i]
i = i + 1
NBs[NBcount, 3] = sum((peaks > NBs[NBcount, 0]) & (peaks < i))
NBs[NBcount, 1] = i
NBcount = NBcount + 1
NBs[NB_count, 3] = sum((peaks > NBs[NB_count, 0]) & (peaks < i))
NBs[NB_count, 1] = i
NB_count = NB_count + 1
else:
i = i + 1
NBs = NBs[0:NBcount, :]
NBs = NBs[0:NB_count, :]
MNBR = NBcount * 60 * second / (simtime - transient)
NBdurations = (array(NBs[:, 1]) - array(NBs[:, 0])) * timeBin / second
MNBR = NB_count * 60 * second / (simtime - transient)
NBdurations = (array(NBs[:, 1]) - array(NBs[:, 0])) * timbin / second
MNBD = mean(NBdurations)
PSIB = sum(NBs[:, 2] / len(APs_wot)) * 100
MFR = len(APs_wot) / ((simtime - transient) / second) / numelectrodes
IBI = (array(NBs[1:, 0]) - array(NBs[0:-1, 1])) * timeBin / second
IBI = (array(NBs[1:, 0]) - array(NBs[0:-1, 1])) * timbin / second
CVIBI = np.std(IBI) / np.mean(IBI)
if NBcount == 0:
if NB_count == 0:
MNBD = 0.0
MNMBs = 0.0
NFBs = 0
else:
NFBs = sum(NBs[:, 3]) / NBcount
NFBs = sum(NBs[:, 3]) / NB_count
if NBcount < 2:
if NB_count < 2:
CVIBI = 0.0
# MAC as defined by Maheswaranathan
# Calculate MAC metric as defined by Maheswaranathan
yf = fft(spikeratetot)
xf = fftfreq(len(spikeratetot), timeBin / second)[:len(spikeratetot) // 2]
xf = fftfreq(len(spikeratetot), timbin / second)[:len(spikeratetot) // 2]
MAC = max(np.abs(yf[1:len(spikeratetot)])) / np.abs(yf[0])
# Calculate cross-correlation between binarized spike trains
# binarize
# Binarize the spike trains
all_combinations = list(combinations(list(arange(numelectrodes)), 2))
Transtimebin = 0.2 * second # timebin to transform spiketrains to binary
bintimeseries = list(range(0, int((simtime - transient) / ms), int(Transtimebin / ms)))
binarysignal = zeros((numelectrodes, len(bintimeseries)))
trans_timebin = 0.2 * second # timebin to transform spiketrains to binary
bin_timeseries = list(range(0, int((simtime - transient) / ms), int(trans_timebin / ms)))
binary_signal = zeros((numelectrodes, len(bin_timeseries)))
for i in range(numelectrodes):
signal = APsinbin[i, :]
grosssignal = [sum(signal[x:x + int((Transtimebin / timeBin))]) for x in
range(0, len(signal), int(Transtimebin / timeBin))]
binarysignal[i, :] = [1 if x > 0 else 0 for x in grosssignal]
signal = APs_inbin[i, :]
grosssignal = [sum(signal[x:x + int((trans_timebin / timbin))]) for x in
range(0, len(signal), int(trans_timebin / timbin))]
binary_signal[i, :] = [1 if x > 0 else 0 for x in grosssignal]
# Calculate coefficients between every pair of electrodes
coefficients = []
N = len(binarysignal[0, :])
N = len(binary_signal[0, :])
for i, j in all_combinations:
signal1 = binarysignal[i, :]
signal2 = binarysignal[j, :]
signal1 = binary_signal[i, :]
signal2 = binary_signal[j, :]
if (i != j) & (not list(signal1) == list(signal2)):
coefficients.append((N * sum(signal1 * signal2) - sum(signal1) * (sum(signal2)))
* ((N * sum(signal1 ** 2) - sum(signal1) ** 2) ** (-0.5))
* ((N * sum(signal2 ** 2) - sum(signal2) ** 2) ** (-0.5)))
meancorr = mean(coefficients)
sdcorr = std(coefficients)
mean_corr = mean(coefficients)
sd_corr = std(coefficients)
if not coefficients:
meancorr = 1
sdcorr = 0
mean_corr = 1
sd_corr = 0
# Compute continuous ISI arrays
time_vector = np.arange(0, (simtime - transient) / second, 1/fs)
......@@ -322,16 +322,17 @@ def ComputeFeatures(APs, simtime, transient, fs):
# Compute ISI measures
meanisi_array = np.nanmean(isi_arrays, axis=0)
meanISI = np.nanmean(meanisi_array)
sdmeanISI = np.nanstd(meanisi_array)
sdtimeISI = np.nanmean(np.nanstd(isi_arrays, axis=0))
mean_ISI = np.nanmean(meanisi_array)
sdmean_ISI = np.nanstd(meanisi_array)
sdtime_ISI = np.nanmean(np.nanstd(isi_arrays, axis=0))
# Calculate the ISI-distance and ISI correlations
isi_distances = np.zeros(len(all_combinations))
ISI_distances = np.zeros(len(all_combinations))
isicoefficients = np.zeros(len(all_combinations))
N = len(isi_arrays[0, :])
j = 0
# iterate through the electrode combinations
# Iterate through the electrode combinations
for electrode1_key, electrode2_key in all_combinations:
# Get the ISI arrays for the selected electrodes
isit1_wn = isi_arrays[electrode1_key, :]
......@@ -342,7 +343,7 @@ def ComputeFeatures(APs, simtime, transient, fs):
isit1 = isit1[~isnan(isit2_wn)]
isi_diff = isit1 / isit2
isi_distances[j] = np.mean(np.where(isi_diff <= 1, abs(isi_diff - 1), -1 * (1 / isi_diff - 1)))
ISI_distances[j] = np.mean(np.where(isi_diff <= 1, abs(isi_diff - 1), -1 * (1 / isi_diff - 1)))
if (i != j) & (not list(isit1) == list(isit2)):
isicoefficients[j] = ((N * sum(isit1 * isit2) - sum(isit1) * (sum(isit2)))
......@@ -351,10 +352,10 @@ def ComputeFeatures(APs, simtime, transient, fs):
j += 1
meanisicorr = mean(isicoefficients)
sdisicorr = std(isicoefficients)
isi_distance = np.mean(isi_distances)
mean_ISIcorr = mean(isicoefficients)
sd_ISIcorr = std(isicoefficients)
ISI_distance = np.mean(ISI_distances)
return [MFR, MNBR, MNBD, PSIB, NFBs, CVIBI, meancorr, sdcorr, meanisicorr, sdisicorr, isi_distance, meanISI, sdmeanISI, sdtimeISI, MAC]
return [MFR, MNBR, MNBD, PSIB, NFBs, CVIBI, mean_corr, sd_corr, mean_ISIcorr, sd_ISIcorr, ISI_distance, mean_ISI, sdmean_ISI, sdtime_ISI, MAC]
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment