"""
Tests of features described in Druckmann et. al. 2013 (https://academic.oup.com/cercor/article/23/12/2994/470476)
AP analysis details (from suplementary info): https://github.com/scidash/neuronunit/files/2295064/bhs290supp.pdf
Numbers in class names refer to the numbers in the publication table
"""
from elephant.spike_train_generation import threshold_detection
from neo import AnalogSignal
from numba import jit
from .base import np, pq, ncap, VmTest, scores
# How this file is different to the original.
# Some big functions broken into smaller ones, for greater modularity
# use of numba jit where possible.
# Different inheritance designed to work with optimizer.
per_ms = pq.UnitQuantity('per_ms',1.0/pq.ms,symbol='per_ms')
none_score = {
'mean': None,
'std': None,
'n': 0
}
debug = False #True
[docs]@jit
def get_diff_spikes(vm):
differentiated = np.diff(vm)
spikes = len([np.any(differentiated) > 0.000143667327364])
return spikes
[docs]@jit
def get_diff(vm,axis=None):
if axis is not None:
differentiated = np.diff(vm,axis=axis)
else:
differentiated = np.diff(vm)
return differentiated
[docs]class Druckmann2013AP:
"""
This is a helper class that computes/finds aspects of APs as defined in Druckmann 2013
"""
def __init__(self, waveform, begin_time):
self.waveform = waveform
self.begin_time = begin_time
self.begin_time.units = pq.ms
[docs] def get_beginning(self):
"""
The beginning of a spike was then determined by a crossing of a threshold on the derivative of the voltage (12mV/msec).
:return: the voltage and time of the AP beginning
"""
begining_time = self.begin_time
beginning_voltage = self.waveform[0]
return beginning_voltage, begining_time
[docs] def get_amplitude(self):
"""
The amplitude of a spike is given by the difference between the voltage at the beginning and peak of the spike.
:return: the amplitude value
"""
v_begin, _ = self.get_beginning()
v_peak, _ = self.get_peak()
return (v_peak - v_begin)
[docs] def get_halfwidth(self):
"""
Amount of time in between the first crossing (in the upwards direction) of the
half-height voltage value and the second crossing (in the downwards direction) of
this value, for the first AP. Half-height voltage is the voltage at the beginning
of the AP plus half the AP amplitude.
:return:
"""
v_begin, _ = self.get_beginning()
amp = self.get_amplitude()
half_v = v_begin + amp / 2.0
above_half_v = np.where(self.waveform.magnitude > half_v)[0]
half_start = self.waveform.times[above_half_v[0]]
half_end = self.waveform.times[above_half_v[-1]]
half_width = half_end - half_start
half_width.units = pq.ms
return half_width
[docs] def get_peak(self):
"""
The peak point of the spike is the maximum in between the beginning and the end.
:return: the voltage and time of the peak
"""
if not hasattr(self, 'peak'):
value = self.waveform.max()
time = self.begin_time + self.waveform.times[np.where(self.waveform.magnitude == value)[0]]
time.units = pq.ms
self.peak = { 'value': value, 'time': time }
return self.peak['value'], self.peak['time']
[docs] def get_trough(self):
peak_v, peak_t = self.get_peak()
post_peak_waveform = self.waveform.magnitude[np.where(self.waveform.times > (peak_t - self.begin_time))]
post_peak_waveform = AnalogSignal(post_peak_waveform, units=self.waveform.units, sampling_period=self.waveform.sampling_period)
value = post_peak_waveform.min()
time = peak_t + post_peak_waveform.times[np.where(post_peak_waveform.magnitude == value)[0]]
time = time[0]
time.units = pq.ms
return value, time
[docs]def isolate_code_block(threshold_crosses,start_time,dvdt_threshold_crosses,dvdt_zero_crosses,vm):
'''
The introduction of this function is was not syntactically necissated. The reason for this functions existence
is to support code modularity.
'''
threshold_crosses = threshold_crosses[np.where(threshold_crosses > start_time)]
dvdt_threshold_crosses = dvdt_threshold_crosses[np.where(dvdt_threshold_crosses > start_time)]
dvdt_zero_crosses = dvdt_zero_crosses[np.where(dvdt_zero_crosses > start_time)]
# Normally, there should be at least as many dvdt threshold crosses as there are v threshold crosses
if len(dvdt_threshold_crosses) < len(threshold_crosses):
dvdt_threshold_crosses = threshold_crosses # for slowly rising APs (e.g. muscle) use the vm threshold as the beginning
ap_beginnings = []
prev_beginning = start_time
prev_threshold = start_time
vm_chopped = 0
for ti, curr_thresh in enumerate(threshold_crosses):
prev_dvdt_zero = dvdt_zero_crosses[np.where(dvdt_zero_crosses < curr_thresh)]
if len(prev_dvdt_zero) == 0:
prev_dvdt_zero = start_time
else:
prev_dvdt_zero = prev_dvdt_zero[-1]
earliest_dvdt_thresh_since_prev_ap = dvdt_threshold_crosses[
np.where((dvdt_threshold_crosses > prev_beginning) & (dvdt_threshold_crosses > prev_threshold) & (dvdt_threshold_crosses > prev_dvdt_zero))
]
if len(earliest_dvdt_thresh_since_prev_ap) != 0:
earliest_dvdt_thresh_since_prev_ap = earliest_dvdt_thresh_since_prev_ap[0]
else:
if ti == 0:
earliest_dvdt_thresh_since_prev_ap = prev_beginning
else:
raise Exception("Did not find a dvdt threshold crossing since previous AP")
ap_beginnings.append(earliest_dvdt_thresh_since_prev_ap)
prev_beginning = earliest_dvdt_thresh_since_prev_ap
prev_threshold = curr_thresh
# The number of ap beginnings should match the number aps detected
assert len(np.unique(ap_beginnings)) == len(threshold_crosses)
vm_mag = vm.magnitude
vm_times = vm.times
vm_chopped = np.split(vm_mag, np.isin(vm_times, ap_beginnings).nonzero()[0])
# The waveform should be cut into APs+1 pieces (1st waveform is steady state)
assert len(vm_chopped) == len(threshold_crosses)+1
return vm_chopped, threshold_crosses, ap_beginnings, vm_mag, vm_times
[docs]class Druckmann2013Test(VmTest):
"""
All tests inheriting from this class assume that the subject model:
1. Is at steady state at time 0 (i.e. resume from SS)
2. Starting at t=0, will have a 2s step current injected into soma, at least once
"""
required_capabilities = (ncap.ProducesActionPotentials,)
score_type = scores.ZScore
def __init__(self, current_amplitude, **params):
super(Druckmann2013Test, self).__init__(**params)
self.params = {
'injected_square_current': {
'delay': 1000 * pq.ms,
'duration': 2000 * pq.ms,
'amplitude': current_amplitude
},
'threshold': -20 * pq.mV,
'beginning_threshold': 12.0 * pq.mV/pq.ms,
'ap_window': 10 * pq.ms,
'repetitions': 1,
}
# This will be an array that stores DruckmannAPs
self.APs = None
[docs] def generate_prediction(self, model):
results = []
reps = self.params['repetitions']
for rep in range(reps):
pred = self.generate_repetition_prediction(model)
results.append(pred)
if reps > 1:
return self.aggregate_repetitions(results)
else:
return results[0]
[docs] def generate_repetition_prediction(self, model):
raise NotImplementedError()
[docs] def aggregate_repetitions(self, results):
values = [rep['mean'] for rep in results if rep['mean'] is not None]
units = values[0].units if len(values) > 0 else self.units
if len(values) > 0:
return {
'mean': np.mean(values) * units,
'std': np.std(values) * units,
'n': len(results)
}
return none_score
[docs] def current_length(self):
return self.params['injected_square_current']['duration']
[docs] def get_APs(self, model):
"""
Spikes were detected by a crossing of a voltage threshold (-20 mV).
:param model: model which provides the waveform to analyse
:return: a list of Druckman2013APs
"""
vm = model.get_membrane_potential()
vm_times = vm.times
start_time = self.params['injected_square_current']['delay'].rescale('sec')
end_time = start_time + self.params['injected_square_current']['duration'].rescale('sec')
vm = AnalogSignal(vm.magnitude[np.where(vm_times <= end_time)], sampling_period=vm.sampling_period, units=vm.units)
try:
dvdt = np.array(np.append([0], get_diff(vm, axis=0))) * pq.mV / vm.sampling_period
except:
dvdt = np.array(np.append([0], get_diff(vm))) * pq.mV / vm.sampling_period
dvdt = AnalogSignal(dvdt, sampling_period=vm.sampling_period)
threshold_crosses = threshold_detection(vm,threshold=self.params['threshold'])
dvdt_threshold_crosses = threshold_detection(dvdt,threshold=self.params['beginning_threshold'])
dvdt_zero_crosses = threshold_detection(dvdt, threshold=0 * pq.mV/pq.ms)
vm_chopped, threshold_crosses, ap_beginnings, vm_mag, vm_times = isolate_code_block(
threshold_crosses, \
start_time,dvdt_threshold_crosses,dvdt_zero_crosses,vm \
)
ap_waveforms = []
for i, b in enumerate(ap_beginnings):
if i != len(ap_beginnings)-1:
waveform = vm_chopped[i+1]
else:
# Keep up to 100ms of the last AP
waveform = vm_mag[np.where((vm_times >= b) & (vm_times < b + 100.0*pq.ms))]
waveform = AnalogSignal(waveform, units=vm.units, sampling_rate=vm.sampling_rate)
ap_waveforms.append(waveform)
# Pass in the AP waveforms and the times when they occured
self.APs = []
for i, b in enumerate(ap_beginnings):
self.APs.append(Druckmann2013AP(ap_waveforms[i], ap_beginnings[i]))
return self.APs
[docs] def get_ISIs(self, model=None):
aps = self.get_APs(model)
ap_times = np.array([ap.get_beginning()[1] for ap in aps])
isis = get_diff(ap_times)# np.diff(ap_times)
return isis
[docs]class AP12AmplitudeDropTest(Druckmann2013Test):
"""
1. Drop in AP amplitude (amp.) from first to second spike (mV)
Difference in the voltage value between the amplitude of the first and second AP.
Negative values indicate 2nd AP amplitude > 1st
"""
name = "Drop in AP amplitude from 1st to 2nd AP"
description = "Difference in the voltage value between the amplitude of the first and second AP"
units = pq.mV
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) >= 2:
if debug:
from matplotlib import pyplot as plt
plt.plot(aps[0].waveform)
plt.plot(aps[1].waveform)
plt.show()
return {
'mean': aps[0].get_amplitude() - aps[1].get_amplitude(),
'std': 0,
'n': 1
}
else:
return none_score
[docs]class AP1SSAmplitudeChangeTest(Druckmann2013Test):
"""
2. AP amplitude change from first spike to steady-state (mV)
Steady state AP amplitude is calculated as the mean amplitude of the set of APs
that occurred during the latter third of the current step.
"""
name = "AP amplitude change from 1st AP to steady-state"
description = """Steady state AP amplitude is calculated as the mean amplitude of the set of APs
that occurred during the latter third of the current step."""
units = pq.mV
[docs] def generate_prediction(self, model):
current_start = self.params['injected_square_current']['delay']
start_latter_3rd = current_start + self.current_length() * 2.0 / 3.0
end_latter_3rd = current_start + self.current_length()
aps = self.get_APs(model)
amps = np.array([ap.get_amplitude() for ap in aps]) * pq.mV
ap_times = np.array([ap.get_beginning()[1] for ap in aps]) * pq.ms
ss_aps = np.where(
(ap_times >= start_latter_3rd) &
(ap_times <= end_latter_3rd))
ss_amps = amps[ss_aps]
if len(aps) > 0 and len(ss_amps) > 0:
if debug:
from matplotlib import pyplot as plt
plt.plot(aps[0].waveform)
for i in ss_aps[0]:
plt.plot(aps[i].waveform)
plt.show()
return {
'mean': amps[0] - ss_amps.mean(),
'std': ss_amps.std(),
'n': len(ss_amps)
}
return none_score
[docs]class AP1AmplitudeTest(Druckmann2013Test):
"""
3. AP 1 amplitude (mV)
Amplitude of the first AP.
"""
name = "First AP amplitude"
description = "Amplitude of the first AP"
units = pq.mV
[docs] def generate_prediction(self, model, ap_index=0):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) > ap_index:
amp = aps[ap_index].get_amplitude()
assert 0 * self.units < amp < 200 * self.units
return {
'mean': amp,
'std': 0,
'n': 1
}
else:
return none_score
[docs]class AP1WidthHalfHeightTest(Druckmann2013Test):
"""
4. AP 1 width at half height (ms)
"""
name = "First AP width at its half height"
description = """Amount of time in between the first crossing (in the upwards direction) of the
half-height voltage value and the second crossing (in the downwards direction) of
this value, for the first AP. Half-height voltage is the voltage at the beginning of
the AP plus half the AP amplitude."""
units = pq.ms
[docs] def generate_prediction(self, model, ap_index=0):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) > ap_index:
hw = aps[ap_index].get_halfwidth()
assert 0 * self.units < hw < 100 * self.units
return {
'mean': hw,
'std': 0,
'n': 1
}
return none_score
[docs]class AP1WidthPeakToTroughTest(Druckmann2013Test):
"""
5. AP 1 peak to trough time (ms)
Amount of time between the peak of the first AP and the trough, i.e., the
minimum of the AHP.
"""
name = "AP 1 peak to trough time"
description = """Amount of time between the peak of the first AP and the trough, i.e., the minimum of the AHP"""
units = pq.ms
[docs] def generate_prediction(self, model, ap_index=0):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) > ap_index:
ap = aps[ap_index]
_, peak_t = ap.get_peak()
_, trough_t = ap.get_trough()
width = trough_t - peak_t
if debug:
from matplotlib import pyplot as plt
plt.plot(aps[0].waveform)
plt.xlim(0, 1000)
plt.show()
assert 0 * self.units <= width < 100 * self.units
return {
'mean': width,
'std': 0,
'n': 1
}
return none_score
[docs]class AP1RateOfChangePeakToTroughTest(Druckmann2013Test):
"""
6. AP 1 peak to trough rate of change (mV/ms)
Difference in voltage value between peak and trough divided by the amount of time in
between the peak and trough.
"""
name = "AP 1 peak to trough rate of change"
description = """Difference in voltage value between peak and trough over the amount of time in between the peak and trough."""
units = pq.mV/pq.ms
[docs] def generate_prediction(self, model, ap_index=0):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) > ap_index:
ap = aps[ap_index]
peak_v, peak_t = ap.get_peak()
trough_v, trough_t = ap.get_trough()
width = trough_t - peak_t
if width == 0 * pq.ms:
width = ap.waveform.sampling_period
change = (trough_v - peak_v) / width
assert change < 0 * self.units
return {
'mean': change,
'std': 0,
'n': 1
}
return none_score
[docs]class AP1AHPDepthTest(Druckmann2013Test):
"""
7. AP 1 Fast AHP depth (mV)
Difference between the minimum of voltage at the trough and the voltage value at
the beginning of the AP.
"""
name = "AP 1 Fast AHP depth"
description = """Difference between the minimum of voltage at the trough and the voltage value at
the beginning of the AP."""
units = pq.mV
[docs] def generate_prediction(self, model, ap_index=0):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) > ap_index:
ap = aps[ap_index]
begin_v, _ = ap.get_beginning()
trough_v, _ = ap.get_trough()
change = begin_v - trough_v
if debug:
from matplotlib import pyplot as plt
plt.plot(aps[0].waveform)
plt.xlim(0, 1000)
plt.show()
return {
'mean': change,
'std': 0,
'n': 1
}
else:
return none_score
[docs]class AP2AmplitudeTest(AP1AmplitudeTest):
"""
8. AP 2 amplitude (mV)
Same as :any:`AP1AmplitudeTest` but for second AP
"""
name = "AP 2 amplitude"
description = """Same as :any:`AP1AmplitudeTest` but for second AP"""
[docs] def generate_prediction(self, model):
return super(AP2AmplitudeTest, self).generate_prediction(model, ap_index=1)
[docs]class AP2WidthHalfHeightTest(AP1WidthHalfHeightTest):
"""
9. AP 2 width at half height (ms)
Same as :any:`AP1WidthHalfHeightTest` but for second AP
"""
name = "AP 2 width at half height"
description = """Same as :any:`AP1WidthHalfHeightTest` but for second AP"""
[docs] def generate_prediction(self, model):
return super(AP2WidthHalfHeightTest, self).generate_prediction(model, ap_index=1)
[docs]class AP2WidthPeakToTroughTest(AP1WidthPeakToTroughTest):
"""
10. AP 2 peak to trough time (ms)
Same as :any:`AP1WidthPeakToTroughTest` but for second AP
"""
name = "AP 2 peak to trough time"
description = """Same as :any:`AP1WidthPeakToTroughTest` but for second AP"""
[docs] def generate_prediction(self, model):
return super(AP2WidthPeakToTroughTest, self).generate_prediction(model, ap_index=1)
[docs]class AP2RateOfChangePeakToTroughTest(AP1RateOfChangePeakToTroughTest):
"""
11. AP 2 peak to trough rate of change (mV/ms)
Same as :any:`AP1RateOfChangePeakToTroughTest` but for second AP
"""
name = "AP 2 peak to trough rate of change"
description = """Same as :any:`AP1RateOfChangePeakToTroughTest` but for second AP"""
[docs] def generate_prediction(self, model):
return super(AP2RateOfChangePeakToTroughTest, self).generate_prediction(model, ap_index=1)
[docs]class AP2AHPDepthTest(AP1AHPDepthTest):
"""
12. AP 2 Fast AHP depth (mV)
Same as :any:`AP1AHPDepthTest` but for second AP
"""
name = "AP 2 Fast AHP depth"
description = """Same as :any:`AP1AHPDepthTest` but for second AP"""
[docs] def generate_prediction(self, model):
return super(AP2AHPDepthTest, self).generate_prediction(model, ap_index=1)
[docs]class AP12AmplitudeChangePercentTest(Druckmann2013Test):
"""
13. Percent change in AP amplitude, first to second spike (%)
Difference in AP amplitude between first and second AP divided by the first AP
amplitude.
"""
name = "Percent change in AP amplitude, first to second spike "
description = """Difference in AP amplitude between first and second AP divided by the first AP
amplitude."""
units = pq.dimensionless
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) >= 2:
amp = self.params['injected_square_current']['amplitude']
amp1 = AP1AmplitudeTest(amp).generate_prediction(model)["mean"]
amp2 = AP2AmplitudeTest(amp).generate_prediction(model)["mean"]
change = (amp2-amp1)/amp1 * 100.0;
return {
'mean': change,
'std': 0,
'n': 1
}
else:
return none_score
[docs]class AP12HalfWidthChangePercentTest(Druckmann2013Test):
"""
14. Percent change in AP width at half height, first to second spike (%)
Difference in AP width at half-height between first and second AP divided by the
first AP width at half-height.
"""
name = "Percent change in AP width at half height, first to second spike"
description = """Difference in AP width at half-height between first and second AP divided by the
first AP width at half-height."""
units = pq.dimensionless
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) >= 2:
amp = self.params['injected_square_current']['amplitude']
width1 = AP1WidthHalfHeightTest(amp).generate_prediction(model)["mean"]
width2 = AP2WidthHalfHeightTest(amp).generate_prediction(model)["mean"]
change = (width2-width1)/width1 * 100.0;
return {
'mean': change,
'std': 0,
'n': 1
}
else:
return none_score
[docs]class AP12RateOfChangePeakToTroughPercentChangeTest(Druckmann2013Test):
"""
15. Percent change in AP peak to trough rate of change, first to second spike (%)
Difference in peak to trough rate of change between first and second AP divided
by the first AP peak to trough rate of change.
"""
name = "Percent change in AP peak to trough rate of change, first to second spike"
description = """Difference in peak to trough rate of change between first and second AP divided
by the first AP peak to trough rate of change."""
units = pq.dimensionless
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) >= 2:
amp = self.params['injected_square_current']['amplitude']
roc1 = AP1RateOfChangePeakToTroughTest(amp).generate_prediction(model)["mean"]
roc2 = AP2RateOfChangePeakToTroughTest(amp).generate_prediction(model)["mean"]
change = (roc2-roc1)/roc1 * 100.0;
return {
'mean': change,
'std': 0,
'n': 1
}
else:
return none_score
[docs]class AP12AHPDepthPercentChangeTest(Druckmann2013Test):
"""
16 Percent change in AP fast AHP depth, first to second spike (%)
Difference in depth of fast AHP between first and second AP divided by the first
AP depth of fast AHP.
"""
name = "Percent change in AP fast AHP depth, first to second spike"
description = """Difference in depth of fast AHP between first and second AP divided by the first
AP depth of fast AHP."""
units = pq.dimensionless
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) >= 2:
amp = self.params['injected_square_current']['amplitude']
ap1 = AP1AHPDepthTest(amp).generate_prediction(model)["mean"]
ap2 = AP2AHPDepthTest(amp).generate_prediction(model)["mean"]
change = (ap2-ap1)/ap1 * 100.0;
return {
'mean': change,
'std': 0,
'n': 1
}
else:
return none_score
[docs]class AP1DelayMeanTest(Druckmann2013Test):
"""
18 Average delay to AP 1 (ms)
Mean of the delay to beginning of first AP over experimental repetitions of step
currents.
"""
name = "First AP delay mean"
description = "Mean delay to the first AP"
units = pq.ms
def __init__(self, current_amplitude, repetitions=7, **params):
super(AP1DelayMeanTest, self).__init__(current_amplitude, **params)
self.params['repetitions'] = repetitions
[docs] def generate_repetition_prediction(self, model, ap_index=0):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) > ap_index:
delay = self.params['injected_square_current']['delay']
if debug:
from matplotlib import pyplot as plt
vm = model.get_membrane_potential()
plt.plot(vm.times.magnitude, vm.magnitude)
plt.xlim(1, aps[ap_index].get_beginning()[1].rescale('sec').magnitude + 0.1)
plt.show()
ap_delay = aps[ap_index].get_beginning()[1] - delay
assert ap_delay > -1 * self.units
return {
'mean': ap_delay,
'std': 0,
'n': 1
}
return none_score
[docs]class AP1DelaySDTest(AP1DelayMeanTest):
"""
19 SD of delay to AP 1 (ms)
Standard deviation of the delay to beginning of first AP over experimental
repetitions of step currents.
"""
name = "First AP delay standard deviation"
description = "Standard deviation of delay to the first AP"
units = pq.ms
[docs] def aggregate_repetitions(self, results):
aggregate = super(AP1DelaySDTest, self).aggregate_repetitions(results)
if aggregate['mean'] is not None:
aggregate['mean'] = aggregate['std']
aggregate['std'] = 0 * self.units
assert aggregate['mean'] >= 0 * self.units
return aggregate
return none_score
[docs]class AP2DelayMeanTest(AP1DelayMeanTest):
"""
20 Average delay to AP 2 (ms)
Same as :any:`AP1DelayMeanTest` but for 2nd AP
"""
name = "Second AP delay mean"
description = "Mean of delay to the second AP"
[docs] def generate_repetition_prediction(self, model, ap_index=0):
return super(AP2DelayMeanTest, self).generate_repetition_prediction(model, ap_index=1)
[docs]class AP2DelaySDTest(AP1DelaySDTest):
"""
21 SD of delay to AP 2 (ms)
Same as :any:`AP1DelaySDTest` but for 2nd AP
Only stochastic models will have a non-zero value for this test
"""
name = "Second AP delay standard deviation"
description = "Standard deviation of delay to the second AP"
[docs] def generate_repetition_prediction(self, model, ap_index=0):
return super(AP2DelaySDTest, self).generate_repetition_prediction(model, ap_index=1)
[docs]class Burst1ISIMeanTest(Druckmann2013Test):
"""
22 Average initial burst interval (ms)
Initial burst interval is defined as the average of the first two ISIs, i.e., the average
of the time differences between the first and second AP and the second and third
AP. This feature is the average the initial burst interval across experimental
repetitions.
"""
name = "Initial burst interval mean"
description = "Mean of the initial burst interval"
units = pq.ms
def __init__(self, current_amplitude, repetitions=7, **params):
super(Burst1ISIMeanTest, self).__init__(current_amplitude, **params)
self.params['repetitions'] = repetitions
[docs] def generate_repetition_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
if len(aps) >= 3:
t1 = aps[0].get_beginning()[1]
t2 = aps[1].get_beginning()[1]
t3 = aps[2].get_beginning()[1]
isi1 = t2 - t1
isi2 = t3 - t2
if debug:
print("first 3 aps: %s, %s, %s"%(t1,t2,t3))
print("2 isis: %s, %s" % (isi1, isi2))
isi_mean = (isi1 + isi2) / 2.0
assert isi_mean > 0 * self.units
return {
'mean': isi_mean,
'std': 0,
'n': 1
}
return none_score
[docs]class Burst1ISISDTest(Burst1ISIMeanTest):
"""
23 SD of average initial burst interval (ms)
The standard deviation of the initial burst interval across experimental repetitions.
"""
name = "Initial burst interval std"
description = "StDev of the initial burst interval"
[docs] def aggregate_repetitions(self, results):
aggregate = super(Burst1ISISDTest, self).aggregate_repetitions(results)
if aggregate['mean'] is not None:
aggregate['mean'] = aggregate['std']
aggregate['std'] = 0 * self.units
assert aggregate['mean'] >= 0 * self.units
return aggregate
return none_score
[docs]class InitialAccommodationMeanTest(Druckmann2013Test):
"""
24 Average initial accommodation (%)
Initial accommodation is defined as the percent difference between the spiking rate of the
*first* fifth of the step current and the *third* fifth of the step current.
"""
name = "Initial accomodation mean"
description = "Mean of the initial accomodation"
units = pq.dimensionless
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
current_start = self.params['injected_square_current']['delay']
start_1st_5th = current_start
end_1st_5th = current_start + self.current_length() * 1/5.0
start_3rd_5th = current_start + self.current_length() * 2/5.0
end_3rd_5th = current_start + self.current_length() * 3/5.0
aps = self.get_APs(model)
ap_times = np.array([ap.get_beginning()[1] for ap in aps])
ap_count15 = np.where((ap_times >= start_1st_5th) & (ap_times <= end_1st_5th))[0]
ap_count35 = np.where((ap_times >= start_3rd_5th) & (ap_times <= end_3rd_5th))[0]
if debug:
print("aps in 1st 5th: %s" % (len(ap_count15)))
print("aps in 3rd 5th: %s" % (len(ap_count35)))
if len(ap_count15) > 0:
percent_diff = (len(ap_count35) - len(ap_count15)) / float(len(ap_count15)) * 100.0
return {
'mean': percent_diff,
'std': 0,
'n': 1
}
return none_score
[docs]class SSAccommodationMeanTest(Druckmann2013Test):
"""
25 Average steady-state accommodation (%)
Steady-state accommodation is defined as the percent difference between the spiking rate
of the *first* fifth of the step current and the *last* fifth of the step current.
"""
name = "Steady state accomodation mean"
description = "Mean of the steady state accomodation"
units = pq.dimensionless
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
current_start = self.params['injected_square_current']['delay']
start_1st_5th = current_start
end_1st_5th = current_start + self.current_length() * 1/5.0
start_last_5th = current_start + self.current_length() * 4/5.0
end_last_5th = current_start + self.current_length()
aps = self.get_APs(model)
ap_times = np.array([ap.get_beginning()[1] for ap in aps])
ap_count15 = np.where((ap_times >= start_1st_5th) & (ap_times <= end_1st_5th))[0]
ap_count55 = np.where((ap_times >= start_last_5th) & (ap_times <= end_last_5th))[0]
if len(ap_count15) > 0:
percent_diff = (len(ap_count55) - len(ap_count15)) / float(len(ap_count15)) * 100.0
if debug:
print("aps in 1st 5th: %s" % (len(ap_count15)))
print("aps in last 5th: %s" % (len(ap_count55)))
return {
'mean': percent_diff,
'std': 0,
'n': 1
}
return none_score
[docs]class AccommodationRateToSSTest(Druckmann2013Test):
"""
26 Rate of accommodation to steady-state (percent/ms)
The percent difference between the spiking rate of the *first* fifth of the step current and
*final* fifth of the step current divided by the time taken to first reach the rate of
steady state accommodation.
Note: It's not clear what is meant by "time taken to first reach the rate of steady state
accommodation". Here, it's computed as smallest t of an ISI which is longer than 0.95 of the
mean ISIs in the final fifth of the current step.
"""
name = "ISI Accomodation Rate"
description = "Rate of ISI Accomodation"
units = per_ms
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
current_start = self.params['injected_square_current']['delay']
start_1st_5th = current_start
end_1st_5th = current_start + self.current_length() * 1/5.0
start_last_5th = current_start + self.current_length() * 4/5.0
end_last_5th = current_start + self.current_length()
aps = self.get_APs(model)
ap_times = np.array([ap.get_beginning()[1] for ap in aps])
aps_15 = np.where((ap_times >= start_1st_5th) & (ap_times <= end_1st_5th))[0]
aps_55 = np.where((ap_times >= start_last_5th) & (ap_times <= end_last_5th))[0]
if len(aps_15) > 0 and len(aps_55) >= 2:
percent_diff = (len(aps_55) - len(aps_15)) / float(len(aps_15)) * 100.0
if debug:
print("aps in 1st 5th vs last 5th, percent change: %s" % (percent_diff))
isis = get_diff(ap_times)
isi_times = ap_times[1:]
isis_55 = isis[np.where((isi_times >= start_last_5th) & (isi_times <= end_last_5th))]
ss_isi = np.mean(isis_55)
if debug:
print("mean ISI at SS: %s" % (ss_isi))
nearly_ss_isis = np.where(isis >= 0.95*ss_isi)[0]
if len(nearly_ss_isis) > 0:
first_nearly_ss_isi_time = isi_times[nearly_ss_isis[0]] * pq.ms
if debug:
print("time of first nearly mean SS ISI: %s" % (first_nearly_ss_isi_time))
if first_nearly_ss_isi_time > 0:
return {
'mean': percent_diff / (first_nearly_ss_isi_time - self.params['injected_square_current']['delay']),
'std': 0,
'n': 1
}
return none_score
[docs]class AccommodationAtSSMeanTest(Druckmann2013Test):
"""
27 Average accommodation at steady-state (%)
Accommodation analysis based on a fit of the ISIs to an exponential function:
ISI = A+B*exp(-t/tau). This feature gives the relative size of the constant term (A) to
the term before the exponent (B).
"""
name = "ISI Steady state accomodation mean"
description = "Mean of the ISI steady state accomodation"
units = pq.dimensionless
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
ap_times = np.array([ap.get_beginning()[1] for ap in aps])
if len(aps) >= 4:
isis = get_diff(ap_times)
isi_delays = ap_times[1:] - self.params['injected_square_current']['delay'].rescale('ms').magnitude
isi_delays = isi_delays - isi_delays[0]
def isi_func(t, A, B, tau):
return A + B * np.exp(-t/(1.0*tau))
from lmfit import Model
model = Model(isi_func)
params = model.make_params(A=isis[-1], B=-1.0, tau=10.0)
params['A'].min = 0
params['B'].max = 0
params['tau'].min = 0
result = model.fit(isis, t=isi_delays, params=params)
A = result.best_values["A"]
B = result.best_values["B"]
tau = result.best_values["tau"]
if debug:
from matplotlib import pyplot as plt
print(result.fit_report())
plt.plot(isi_delays, isis, 'bo')
plt.plot(isi_delays, result.best_fit, 'r-')
plt.show()
return {
'mean': self.get_final_result(A, B, tau),
'std': 0,
'n': 1
}
return none_score
[docs] @jit
def get_final_result(self, A, B, tau):
return B / float(A) * 100.0
[docs]class AccommodationRateMeanAtSSTest(AccommodationAtSSMeanTest):
"""
28 Average rate of accommodation during steady-state
Accommodation analysis based on a fit of the ISIs to an exponential function.
This feature is the time constant (tau) of the exponent.
"""
name = "ISI accomodation time constant"
description = "Time constant of the ISI accomodation"
units = pq.ms
[docs] def get_final_result(self, A, B, tau):
return tau * pq.ms
[docs]class ISICVTest(Druckmann2013Test):
"""
29 Average inter-spike interval (ISI) coefficient of variation (CV) (unitless)
Coefficient of variation (mean divided by standard deviation) of the distribution
of ISIs.
"""
name = "ISI CV"
description = "ISI Coefficient of Variation"
units = pq.dimensionless
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
isis = self.get_ISIs(model)
if len(isis) >= 2:
mean = np.mean(isis)
std = np.std(isis)
if debug:
print("isi mean: %s std: %s"%(mean, std))
if std > 1e-5:
assert mean > 0 * self.units
return {
'mean': mean / std,
'std': 0,
'n': 1
}
return none_score
[docs]class ISIBurstMeanChangeTest(Druckmann2013Test):
"""
31 Average change in ISIs during a burst (%)
Difference between the first and second ISI divided by the value of the first ISI.
"""
name = "ISI Burst Mean Change"
description = "ISI Burst Mean Change"
units = pq.dimensionless
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
isis = self.get_ISIs(model)
if len(isis) >= 2:
if debug:
print("ISI1: %s ISI2: %s Change: %s"%(isis[0],isis[1],isis[1] - isis[0]))
if isis[0] > 1e-5:
return {
'mean': (isis[1] - isis[0])/isis[0] * 100.0,
'std': 0,
'n': 1
}
return none_score
[docs]class SpikeRateStrongStimTest(Druckmann2013Test):
"""
32 Average rate, strong stimulus (Hz)
Firing rate of strong stimulus.
"""
name = "Strong Stimulus Firing Rate"
description = "Strong Stimulus Firing Rate"
units = pq.Hz
[docs] def generate_prediction(self, model):
model.inject_square_current(self.params['injected_square_current'])
aps = self.get_APs(model)
duration = self.current_length()
spike_rate = len(aps) / duration
spike_rate.units = pq.Hz
if debug:
print("APs: %s Duration: %s"%(len(aps), duration))
assert 400 * self.units > spike_rate > 0 * self.units
return {
'mean': spike_rate,
'std': 0,
'n': 1
}
[docs]class AP1DelayMeanStrongStimTest(AP1DelayMeanTest):
"""
33 Average delay to AP 1, strong stimulus (ms)
Same as :any:`AP1DelayMeanTest` but for strong stimulus
"""
[docs]class AP1DelaySDStrongStimTest(AP1DelaySDTest):
"""
34 SD of delay to AP 1, strong stimulus (ms)
Same as :any:`AP1DelaySDTest` but for strong stimulus
"""
[docs]class AP2DelayMeanStrongStimTest(AP2DelayMeanTest):
"""
35 Average delay to AP 2, strong stimulus (ms)
Same as :any:`AP2DelayMeanTest` but for strong stimulus
"""
[docs]class AP2DelaySDStrongStimTest(AP2DelaySDTest):
"""
36 SD of delay to AP 2, strong stimulus (ms)
Same as :any:`AP2DelaySDTest` but for strong stimulus
"""
[docs]class Burst1ISIMeanStrongStimTest(Burst1ISIMeanTest):
"""
37 Average initial burst ISI, strong stimulus (ms)
Same as :any:`Burst1ISIMeanTest` but for strong stimulus
"""
[docs]class Burst1ISISDStrongStimTest(Burst1ISISDTest):
"""
38 SD of average initial burst ISI, strong stimulus (ms)
Same as :any:`Burst1ISISDTest` but for strong stimulus
"""