From aa84b11202cc863553bf330b53f20afb9d2f1048 Mon Sep 17 00:00:00 2001
From: "Doorn, Nina (UT-TNW)" <n.doorn-1@utwente.nl>
Date: Wed, 15 Jan 2025 11:36:52 +0100
Subject: [PATCH] Fixed MEA feature names
---
Simulator.py | 299 ++++++++++++++++++++++++++-------------------------
1 file changed, 150 insertions(+), 149 deletions(-)
diff --git a/Simulator.py b/Simulator.py
index 67fb501..e1cb554 100644
--- a/Simulator.py
+++ b/Simulator.py
@@ -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]
--
GitLab