"""
Functions for STN model
"""
__author__ = "Nils A. Koch"
__copyright__ = 'Copyright 2022, Nils A. Koch'
__license__ = "MIT"
import scipy
import json
import copy
import os
import gc
import h5py
import numpy as np
from numba import types, njit
from numba.typed import Dict
from scipy import signal
from Code.Functions.Utility_fxns import stimulus_init, NumpyEncoder
@njit
def STN(V, g, E, I_in, dt, currents_included, stim_time, stim_num, C, shift, scale, b_param, slope_shift,
gating, current, ind_dict, vol):
""" Simulate current and Vm for STN model from Otsuka et al 2004 (https://dx.doi.org/10.1152/jn.00508.2003)
for I_in
Parameters
----------
V : float
Initial membrane potentia
g : dict
maximal conductance for currents
E : dict
Reversal potentials for currents
I_in : array
Current stimulus input to model
dt : float
fixed time step
currents_included : dict
Boolean to include current in simulation or not
stim_time :
length of stimulus ()
stim_num : int
number of different stimulus currents
C : float
membrane capacitance
shift : dict
for each gating parameter
scale : dict
for each gating parameter
b_param : dict
Boltzmann parameter values for each gating parameter
slope_shift : dict
for each gating parameter
gating : array
array of gating values over time
current : array
array of current values over time
ind_dict: dict
Dictionary of indices in arrays for gating variables
vol : float
cell volume
Returns
-------
V_m : array
Simulated membrane potential in response to I_in
"""
# initialize V_m, Ca_conc and E_Ca
V_m = np.zeros((stim_num, stim_time))
Ca_conc = np.ones((stim_num)) * 5e-3
E_Ca = 0.013319502 * np.log(2000 / Ca_conc) * 1000 * np.ones((stim_num)) # mV
# initialize gating
# Na activation
m_inf = (1 / (1 + np.exp((V - shift["m"] - slope_shift["m"] - b_param["m"][0]) / b_param["m"][1]))) ** b_param["m"][2]
gating[ind_dict['m'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * m_inf
# Na inactivation
h_inf = ((1 - b_param["h"][3]) / (1 + np.exp((V - shift["h"] - slope_shift["h"] - b_param["h"][0]) / b_param["h"][1])) + b_param["h"][3]) \
** b_param["h"][2]
gating[ind_dict['h'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * h_inf
# K activation
n_inf = (1 / (1 + np.exp((V - shift["n"] - slope_shift["n"] - b_param["n"][0]) / b_param["n"][1]))) ** \
b_param["n"][2]
gating[ind_dict['n'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * n_inf
# T-type Ca activation
p_inf = (1 / (1 + np.exp((V - shift["p"] - slope_shift["p"] - b_param["p"][0]) / b_param["p"][1]))) ** \
b_param["p"][2]
gating[ind_dict['p'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * p_inf
# T-type Ca inactivation
q_inf = (1 / (1 + np.exp((V - shift["q"] - slope_shift["q"] - b_param["q"][0]) / b_param["q"][1]))) ** \
b_param["q"][2]
gating[ind_dict['q'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * q_inf
# Ca dependent K activation
r_inf = (1 / (1 + np.exp((Ca_conc - shift["r"] - slope_shift["r"] - b_param["r"][0]) / b_param["r"][1]))) ** \
b_param["r"][2]
gating[ind_dict['r'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * r_inf
# A-type activation
a_inf = (1 / (1 + np.exp((V - shift["a"] - slope_shift["a"] - b_param["a"][0]) / b_param["a"][1]))) ** \
b_param["a"][2]
gating[ind_dict['a'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * a_inf
# A-type inactivation
b_inf = (1 / (1 + np.exp((V - shift["b"] - slope_shift["b"] - b_param["b"][0]) / b_param["b"][1]))) ** \
b_param["b"][2]
gating[ind_dict['b'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * b_inf
# mutant A-type activation
a_mut_inf = (1 / (1 + np.exp((V - shift["a_mut"] - slope_shift["a_mut"] - b_param["a_mut"][0]) / b_param["a_mut"][1]))) ** \
b_param["a_mut"][2]
gating[ind_dict['a_mut'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * a_mut_inf
# mutant A-type inactivation
b_mut_inf = (1 / (1 + np.exp((V - shift["b_mut"] - slope_shift["b_mut"] - b_param["b_mut"][0]) / b_param["b_mut"][1]))) ** \
b_param["b_mut"][2]
gating[ind_dict['b_mut'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * b_mut_inf
# L-type Ca activation
c_inf = (1 / (1 + np.exp((V - shift["c"] - slope_shift["c"] - b_param["c"][0]) / b_param["c"][1]))) ** \
b_param["c"][2]
gating[ind_dict['c'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * c_inf
# L-type Ca inactivation 1
d1_inf = (1 / (1 + np.exp((V - shift["d1"] - slope_shift["d1"] - b_param["d1"][0]) / b_param["d1"][1]))) ** \
b_param["d1"][2]
gating[ind_dict['d1'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * d1_inf
# L-type Ca inactivation 2
d2_inf = (1 / (1 + np.exp((Ca_conc - shift["d2"] - slope_shift["d2"] - b_param["d2"][0]) / b_param["d2"][1]))) ** \
b_param["d2"][2]
gating[ind_dict['d2'], :] = np.ones((stim_num), dtype=np.dtype('float64')) * d2_inf
# initialize gating and V_m
current[ind_dict['Na'], :] = np.zeros((stim_num), dtype=np.dtype('float64'))
current[ind_dict['Kd'], :] = np.zeros((stim_num), dtype=np.dtype('float64'))
current[ind_dict['A'], :] = np.zeros((stim_num), dtype=np.dtype('float64'))
current[ind_dict['A_mut'], :] = np.zeros((stim_num), dtype=np.dtype('float64'))
current[ind_dict['Leak'], :] = np.zeros((stim_num), dtype=np.dtype('float64'))
current[ind_dict['L'], :] = np.zeros((stim_num), dtype=np.dtype('float64'))
current[ind_dict['T'], :] = np.zeros((stim_num), dtype=np.dtype('float64'))
V_m[:, 0] = V
gating[ind_dict['Ca_conc'], :] = Ca_conc
gating[ind_dict['E_Ca'], :] = E_Ca
t = 1
while t < stim_time:
# update Ca conc and E_Ca
Ca_conc_dt = ((-5.182134834753951e-06 / vol) * (current[ind_dict['T'], :] + current[ind_dict['L'], :])) / 1000 - (2 * Ca_conc)
Ca_conc = Ca_conc + Ca_conc_dt * dt
Ca_conc[Ca_conc <= 0] = 0
E_Ca = 0.013319502 * np.log(2000 / Ca_conc) * 1000 # mV
gating[ind_dict['Ca_conc'], :] = Ca_conc
gating[ind_dict['E_Ca'], :] = E_Ca
# Na activation
m_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["m"] - slope_shift["m"] - b_param["m"][0]) / b_param["m"][1]))) ** \
b_param["m"][2]
tau_m = 0.2 + 3 / (1 + np.exp((V_m[:, t - 1] - shift["m"] + 53) / 0.7))
m_dot = (m_inf - gating[ind_dict['m'], :]) * (1 / (tau_m * scale["m"]))
gating[ind_dict['m'], :] = gating[ind_dict['m'], :] + m_dot * dt
# Na inactivation
h_inf = ((1 - b_param["h"][3]) / (1 + np.exp((V_m[:, t - 1] - shift["h"] - slope_shift["h"] - b_param["h"][0])
/ b_param["h"][1])) + b_param["h"][3]) ** b_param["h"][2]
tau_h = 24.5 / (np.exp((V_m[:, t - 1] - shift["h"] + 50) / 15) + np.exp(-(V_m[:, t - 1] - shift["h"] + 50) / 16))
h_dot = (h_inf - gating[ind_dict['h'], :]) * (1 / (tau_h * scale["h"]))
gating[ind_dict['h'], :] = gating[ind_dict['h'], :] + h_dot * dt
# K activation
n_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["n"] - slope_shift["n"] - b_param["n"][0]) / b_param["n"][1]))) ** \
b_param["n"][2]
tau_n = 11 / (np.exp((V_m[:, t - 1] - shift["n"] + 40) / 40) + np.exp(-(V_m[:, t - 1] - shift["n"] + 40) / 50))
n_dot = (n_inf - gating[ind_dict['n'], :]) * (1 / (tau_n * scale["n"]))
gating[ind_dict['n'], :] = gating[ind_dict['n'], :] + n_dot * dt
# T-type Ca activation
p_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["p"] - slope_shift["p"] - b_param["p"][0]) / b_param["p"][1]))) ** \
b_param["p"][2]
tau_p = 5 + 0.33 / (np.exp((V_m[:, t - 1] - shift["p"] + 27) / 10) + np.exp(-(V_m[:, t - 1] - shift["p"] + 102) / 15))
p_dot = (p_inf - gating[ind_dict['p'], :]) * (1 / (tau_p * scale["p"]))
gating[ind_dict['p'], :] = gating[ind_dict['p'], :] + p_dot * dt
# T-type Ca inactivation
q_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["q"] - slope_shift["q"] - b_param["q"][0]) / b_param["q"][1]))) ** \
b_param["q"][2]
tau_q = 400 / (np.exp((V_m[:, t - 1] - shift["q"] + 50) / 15) + np.exp(-(V_m[:, t - 1] - shift["q"] + 50) / 16))
q_dot = (q_inf - gating[ind_dict['q'], :]) * (1 / (tau_q * scale["q"]))
gating[ind_dict['q'], :] = gating[ind_dict['q'], :] + q_dot * dt
# Ca dependent K activation
r_inf = (1 / (1 + np.exp((Ca_conc * 1e-6 - shift["r"] - slope_shift["r"] - b_param["r"][0]) / b_param["r"][1]))) ** \
b_param["r"][2]
tau_r = 2
r_dot = (r_inf - gating[ind_dict['r'], :]) * (1 / (tau_r * scale["r"]))
gating[ind_dict['r'], :] = gating[ind_dict['r'], :] + r_dot * dt
# A-type activation
a_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["a"] - slope_shift["a"] - b_param["a"][0]) / b_param["a"][1]))) ** \
b_param["a"][2]
tau_a = 1 + 1 / (1 + np.exp((V_m[:, t - 1] - shift["a"] + 40) / 0.5))
a_dot = (a_inf - gating[ind_dict['a'], :]) * (1 / (tau_a * scale["a"]))
gating[ind_dict['a'], :] = gating[ind_dict['a'], :] + a_dot * dt
# A-type inactivation
b_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["b"] - slope_shift["b"] - b_param["b"][0]) / b_param["b"][1]))) ** \
b_param["b"][2]
tau_b = 200 / (np.exp((V_m[:, t - 1] - shift["b"] + 60) / 30) + np.exp(-(V_m[:, t - 1] - shift["b"] + 40) / 10))
b_dot = (b_inf - gating[ind_dict['b'], :]) * (1 / (tau_b * scale["b"]))
gating[ind_dict['b'], :] = gating[ind_dict['b'], :] + b_dot * dt
# mutant A-type activation
a_mut_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["a_mut"] - slope_shift["a_mut"] - b_param["a_mut"][0]) / b_param["a_mut"][1]))) ** \
b_param["a_mut"][2]
tau_a_mut = 1 + 1 / (1 + np.exp((V_m[:, t - 1] - shift["a_mut"] + 40) / 0.5))
a_mut_dot = (a_mut_inf - gating[ind_dict['a_mut'], :]) * (1 / (tau_a_mut * scale["a_mut"]))
gating[ind_dict['a_mut'], :] = gating[ind_dict['a_mut'], :] + a_mut_dot * dt
# mutant A-type inactivation
b_mut_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["b_mut"] - slope_shift["b_mut"] - b_param["b_mut"][0]) / b_param["b_mut"][1]))) ** \
b_param["b_mut"][2]
tau_mut_b = 200 / (np.exp((V_m[:, t - 1] - shift["b_mut"] + 60) / 30) + np.exp(-(V_m[:, t - 1] - shift["b_mut"] + 40) / 10))
b_mut_dot = (b_mut_inf - gating[ind_dict['b_mut'], :]) * (1 / (tau_mut_b * scale["b_mut"]))
gating[ind_dict['b_mut'], :] = gating[ind_dict['b_mut'], :] + b_mut_dot * dt
# L-type Ca activation
c_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["c"] - slope_shift["c"] - b_param["c"][0]) / b_param["c"][1]))) ** \
b_param["c"][2]
tau_c = 45 + 10 / (np.exp((V_m[:, t - 1] - shift["c"] + 27) / 20) + np.exp(-(V_m[:, t - 1] - shift["c"] + 50) / 15))
c_dot = (c_inf - gating[ind_dict['c'], :]) * (1 / (tau_c * scale["c"]))
gating[ind_dict['c'], :] = gating[ind_dict['c'], :] + c_dot * dt
# L-type Ca inactivation 1
d1_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["d1"] - slope_shift["d1"] - b_param["d1"][0]) / b_param["d1"][1]))) ** \
b_param["d1"][2]
tau_d1 = 400 + 500 / (np.exp((V_m[:, t - 1] - shift["d1"] + 40) / 15) + np.exp(-(V_m[:, t - 1] - shift["d1"] + 20) / 20))
d1_dot = (d1_inf - gating[ind_dict['d1'], :]) * (1 / (tau_d1 * scale["d1"]))
gating[ind_dict['d1'], :] = gating[ind_dict['d1'], :] + d1_dot * dt
# L-type Ca inactivation 2
d2_inf = (1 / (1 + np.exp((Ca_conc * 1e-6 - shift["d2"] - slope_shift["d2"] - b_param["d2"][0]) / b_param["d2"][1]))) ** \
b_param["d2"][2]
tau_d2 = 130 # dt
d2_dot = (d2_inf - gating[ind_dict['d2'], :]) * (1 / (tau_d2 * scale["d2"]))
gating[ind_dict['d2'], :] = gating[ind_dict['d2'], :] + d2_dot * dt
# update currents
current[ind_dict['Leak'], :] = g["Leak"] * (V_m[:, t - 1] - E["Leak"]) * currents_included["Leak"]
current[ind_dict['Na'], :] = g["Na"] * (gating[ind_dict['m'], :] ** 3) * gating[ind_dict['h'], :] \
* (V_m[:, t - 1] - E["Na"]) * currents_included["Na"]
current[ind_dict['Kd'], :] = g["Kd"] * (gating[ind_dict['n'], :] ** 4) * (V_m[:, t - 1] - E["K"]) * \
currents_included["Kd"]
current[ind_dict['L'], :] = g["L"] * ((gating[ind_dict['c'], :]) ** 2) * gating[ind_dict['d1'],:] \
* gating[ind_dict['d2'],:] * (V_m[:, t - 1] - E_Ca) * currents_included["L"]
current[ind_dict['T'], :] = g["T"] * (gating[ind_dict['p'], :] ** 2) * gating[ind_dict['q'], :] * \
(V_m[:, t - 1] - E_Ca) * currents_included["T"]
current[ind_dict['A'], :] = g["A"] * (gating[ind_dict['a'], :] ** 2) * gating[ind_dict['b'], :] * \
(V_m[:, t - 1] - E["K"]) * currents_included["A"]
current[ind_dict['A_mut'], :] = g["A_mut"] * (gating[ind_dict['a_mut'], :] ** 2) * gating[ind_dict['b_mut'], :] * \
(V_m[:, t - 1] - E["K"]) * currents_included["A_mut"]
current[ind_dict['Ca_K'], :] = g["Ca_K"] * (gating[ind_dict['r'], :]) ** 2 * (V_m[:, t - 1] - E["K"]) * \
currents_included["Ca_K"]
# update dV/dt
V_dot = (1 / C) * (I_in[t, :] - (current[ind_dict['Leak'], :] + current[ind_dict['Na'], :] +
current[ind_dict['Kd'], :] + current[ind_dict['L'], :] +
current[ind_dict['T'], :] + current[ind_dict['A'], :] +
current[ind_dict['A_mut'], :] +current[ind_dict['Ca_K'], :]))
# update V_m
V_m[:, t] = V_m[:, t - 1] + V_dot * dt
t += 1
return V_m
def STN_mut(V_init, g, E, I_in, dt, currents_included, stim_time, stim_num, C, shift, scale, b_param, slope_shift,
gating, current, prominence, desired_AUC_width, mutations, mut, folder, high, low, number_steps,
initial_period, sec, vol):
""" Simulate mutations for STN model
Parameters
----------
V_init : float
Initial membrane potential
g : dict
maximal conductance for currents
E : dict
Reversal potentials for currents
I_in : array
Current stimulus input to model
dt : float
fixed time step
currents_included : dict
Boolean to include current in simulation or not
stim_time :
length of stimulus ()
stim_num : int
number of different stimulus currents
C : float
membrane capacitance
shift : dict
for each gating parameter
scale : dict
for each gating parameter
b_param : dict
Boltzmann parameter values for each gating parameter
slope_shift : dict
for each gating parameter
gating : array
array of gating values over time
current : array
array of current values over time
prominence : float
required spike prominence
desired_AUC_width : float
width of AUC from rheobase to rheobase+desired_AUC_width
mutations : dict
dictionary of mutation effects
mut : str
mutation
folder : str
folder to save hdf5 file in
high : float
upper current step magnitude
low : float
lowest current step magnitude
number_steps : int
number of current steps between low and high
initial_period : float
length of I = 0 before current step
sec:
length of current step in seconds
vol : float
cell volume
Returns
-------
saves mut.hdf5 file to folder
"""
# setup for simulation
g_multi = copy.deepcopy(g)
E_multi = copy.deepcopy(E)
b_param_multi = copy.deepcopy(b_param)
shift_multi = copy.deepcopy(shift)
scale_multi = copy.deepcopy(scale)
slope_shift_multi = copy.deepcopy(slope_shift)
currents_included_multi = copy.deepcopy(currents_included)
gating2 = copy.deepcopy(gating)
current2 = copy.deepcopy(current)
g_effect = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in g_multi.items():
g_effect[k] = v
b_param_effect = Dict.empty(key_type=types.unicode_type, value_type=types.float64[:], )
for k, v in b_param_multi.items():
b_param_effect[k] = np.array(v)
shift2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in shift_multi.items():
shift2[k] = v
scale2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in scale_multi.items():
scale2[k] = v
slope_shift2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in slope_shift_multi.items():
slope_shift2[k] = v
E2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in E_multi.items():
E2[k] = v
currents_included2 = Dict.empty(key_type=types.unicode_type, value_type=types.boolean, )
for k, v in currents_included_multi.items():
currents_included2[k] = v
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(
['m', 'h', 'n', 'a', 'b','a_mut', 'b_mut', 'c', 'd1', 'd2', 'p', 'q', 'r', 'Ca_conc', 'E_Ca']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'A_mut', 'L', 'T', 'Leak', 'Ca_K']):
ind_dict[var] = i
i += 1
#mutation effect
g_effect['A_mut'] = g['A_mut'] * mutations[mut]['g_ratio']
b_param_effect['a_mut'][0] = b_param_effect['a_mut'][0] + mutations[mut]['activation_Vhalf_diff']
b_param_effect['a_mut'][1] = b_param_effect['a_mut'][1] * mutations[mut]['activation_k_ratio']
# initial simulation for current range from low to high
V_m = STN(V_init * np.ones((stim_num)), g_effect, E2, I_in, dt, currents_included2, stim_time, stim_num, C, shift2,
scale2, b_param_effect, slope_shift2, gating2, current2, ind_dict, vol)
stim_start = np.int(initial_period * 1 / dt)
min_spike_height = V_m[0, stim_start] + prominence
# create hdf5 file for the mutation and save data and metadata
fname = os.path.join(folder, "{}.hdf5".format(mut.replace(" ", "_")))
with h5py.File(fname, "a") as f:
data = f.create_group("data")
data.create_dataset('V_m', data=V_m, dtype='float64', compression="gzip", compression_opts=9)
metadata = {'I_high': high, 'I_low': low, 'stim_num': number_steps, 'stim_time': stim_time, 'sec': sec, 'dt': dt,
'initial_period': initial_period, 'g': json.dumps(dict(g_effect)), 'E': json.dumps(dict(E2)), 'C': C,
'V_init': V_init, 'vol': vol, 'ind_dict': json.dumps(dict(ind_dict)),
'b_param': json.dumps(dict(b_param_effect), cls=NumpyEncoder),
'currents': json.dumps(dict(currents_included2)), 'shift': json.dumps(dict(shift2)),
'scale': json.dumps(dict(scale2)), 'slope_shift': json.dumps(dict(slope_shift2)),
'prominence': prominence, 'desired_AUC_width': desired_AUC_width, 'mut': str(mut),
'min_spike_height': min_spike_height,}
data.attrs.update(metadata)
# firing frequency analysis and rheobase
with h5py.File(fname, "r+") as f:
I_mag = np.arange(low, high, (high - low) / stim_num)
analysis = f.create_group("analysis")
dtyp = h5py.special_dtype(vlen=np.dtype('float64'))
for i in np.array(['spike_times', 'ISI', 'Freq','amplitude',]):
analysis.create_dataset(i, (I_mag.shape[0],), dtype=dtyp, compression="gzip", compression_opts=9)
analysis.create_dataset('F_inf', (I_mag.shape[0],), dtype='float64')
analysis.create_dataset('F_null', (I_mag.shape[0],), dtype='float64')
analysis.create_dataset('AP', (I_mag.shape[0],), dtype='bool', compression="gzip", compression_opts=9)
analysis.create_dataset('rheobase', (1,), dtype='float64', compression="gzip", compression_opts=9)
for stim in range(I_mag.shape[0]):
# spike detection
ap_all, ap_prop = scipy.signal.find_peaks(V_m[stim, stim_start:], min_spike_height, prominence=prominence, distance=1 * 1 / dt)
analysis['spike_times'][stim] = ap_all * dt
if analysis['spike_times'][stim].size == 0:
analysis['AP'][stim] = False
else:
analysis['AP'][stim] = True
analysis['amplitude'][stim] = np.array(ap_prop["peak_heights"])
analysis['ISI'][stim] = np.array([x - analysis['spike_times'][stim][i - 1] if i else None for i, x in
enumerate(analysis['spike_times'][stim])][1:]) # msec
analysis['Freq'][stim] = 1 / (analysis['ISI'][stim] / 1000)
if analysis['Freq'][stim].size != 0:
analysis['F_null'][stim] = analysis['Freq'][stim][0]
last_second_spikes = analysis['spike_times'][stim] > np.int(
f['data'].attrs['stim_time'] * dt - (500 + initial_period))
if np.any(last_second_spikes == True):
last_second_spikes_1 = np.delete(last_second_spikes, 0) # remove first spike because it has no ISI/freq
analysis['F_inf'][stim] = np.mean(analysis['Freq'][stim][last_second_spikes_1])
else: # if no ISI detected
analysis['F_null'][stim] = 0
analysis['F_inf'][stim] = 0
# rheobase
# find where the first AP occurs
if np.any(analysis['AP'][:] == True):
# setup current range to be between the current step before and with first AP
I_first_spike = I_mag[np.argwhere(analysis['AP'][:] == True)[0][0]]
I_before_spike = I_mag[np.argwhere(analysis['AP'][:] == True)[0][0] - 1]
analysis.create_dataset('I_first_spike', data=I_first_spike, dtype='float64')
analysis.create_dataset('I_before_spike', data=I_before_spike, dtype='float64')
number_steps_rheo = 100
stim_time, I_in_rheo, stim_num_rheo, V_m_rheo = stimulus_init(I_before_spike, I_first_spike, number_steps_rheo,
initial_period, dt, sec)
# simulate response to this current range and save
V_m_rheo = STN(V_init * np.ones((stim_num_rheo)), g_effect, E2, I_in_rheo, dt, currents_included2,
stim_time, stim_num_rheo, f['data'].attrs['C'], shift2, scale2, b_param_effect, slope_shift2,
np.zeros((len(b_param), stim_num_rheo)), np.zeros((len(currents_included), stim_num_rheo)),
ind_dict, f['data'].attrs['vol'])
f['data'].create_dataset('V_m_rheo', data=V_m_rheo, dtype='float64', compression="gzip", compression_opts=9)
# run AP detection
analysis.create_dataset('spike_times_rheo', (I_in_rheo.shape[1],), dtype=dtyp, compression="gzip",
compression_opts=9)
analysis.create_dataset('AP_rheo', (I_in_rheo.shape[1],), dtype='bool', compression="gzip", compression_opts=9)
for stim in range(I_in_rheo.shape[1]): # range(stim_num):
ap_all, ap_prop = scipy.signal.find_peaks(V_m_rheo[stim, stim_start:], min_spike_height,
prominence=prominence, distance=1 * 1 / dt)
analysis['spike_times_rheo'][stim] = ap_all * dt
if analysis['spike_times_rheo'][stim].size == 0:
analysis['AP_rheo'][stim] = False
else:
analysis['AP_rheo'][stim] = True
# find rheobase as smallest current step to elicit AP and save
analysis['rheobase'][:] = I_in_rheo[-1, np.argwhere(analysis['AP_rheo'][:] == True)[0][0]] * 1000
# AUC
with h5py.File(fname, "a") as f:
F_inf = f['analysis']['F_inf'][:]
try:
# find rheobase of steady state firing (F_inf)
# setup current range to be in current step between no and start of steady state firing
F_inf_start_ind = np.argwhere(F_inf != 0)[0][0]
dI = (high - low) / stim_num
I_start = I_mag[F_inf_start_ind - 1]
I_end = I_mag[F_inf_start_ind]
I_rel = np.arange(I_start, I_end + (I_end - I_start) / 100, (I_end - I_start) / 100)
I_rel = np.reshape(I_rel, (1, I_rel.shape[0]))
I_in = np.zeros((stim_time + np.int(initial_period * 1 / dt), 1)) @ I_rel
I_in[np.int(initial_period * 1 / dt):, :] = np.ones((stim_time, 1)) @ I_rel
stim_num = I_in.shape[1]
# simulate response to this current step range
V_m_inf_rheo = STN(V_init * np.ones(stim_num), g_effect, E2, I_in, dt, currents_included2, stim_time,
stim_num, C, shift2, scale2, b_param_effect, slope_shift2,
np.zeros((len(b_param), stim_num)), np.zeros((len(currents_included), stim_num)),
ind_dict, vol)
# save response, analyse and save
for i in np.array(['spike_times_Finf', 'ISI_Finf', 'Freq_Finf']): # analysis_names:
f['analysis'].require_dataset(i, shape=(I_rel.shape[1],), dtype=dtyp, compression="gzip",
compression_opts=9)
f['analysis'].require_dataset('F_inf_Finf', shape=(I_rel.shape[1],), dtype='float64')
f['analysis'].require_dataset('F_null_Finf', shape=(I_rel.shape[1],), dtype='float64')
for stim in range(I_rel.shape[1]):
# spike peak detection
ap_all, ap_prop = scipy.signal.find_peaks(V_m_inf_rheo[stim, stim_start:], min_spike_height,
prominence=prominence, distance=1 * 1 / dt)
f['analysis']['spike_times_Finf'][stim] = ap_all * dt
f['analysis']['ISI_Finf'][stim] = np.array(
[x - f['analysis']['spike_times_Finf'][stim][i - 1] if i else None for i, x in
enumerate(f['analysis']['spike_times_Finf'][stim])][1:]) # msec
f['analysis']['Freq_Finf'][stim] = 1 / (f['analysis']['ISI_Finf'][stim] / 1000)
if f['analysis']['Freq_Finf'][stim].size != 0:
last_second_spikes = f['analysis']['spike_times_Finf'][stim] > np.int(
f['data'].attrs['stim_time'] * dt - (500 + f['data'].attrs['initial_period']))
if np.any(last_second_spikes == True):
last_second_spikes_1 = np.delete(last_second_spikes,
0) # remove first spike because it has no ISI/freq
f['analysis']['F_inf_Finf'][stim] = np.mean(
f['analysis']['Freq_Finf'][stim][last_second_spikes_1])
else:
f['analysis']['F_null_Finf'][stim] = 0
f['analysis']['F_inf_Finf'][stim] = 0
# find AUC relative to F_inf rheobase
# setup current inputs with current range around F_inf rheobase
I_first_inf = I_rel[0, np.argwhere(f['analysis']['F_inf_Finf'][:] != 0)[0][0]]
I_start = I_first_inf - 0.00005
I_end = I_first_inf + desired_AUC_width * high
I_rel2 = np.arange(I_start, I_end + (I_end - I_start) / 200, (I_end - I_start) / 200)
I_rel2 = np.reshape(I_rel2, (1, I_rel2.shape[0]))
I_in2 = np.zeros((stim_time + np.int(initial_period * 1 / dt), 1)) @ I_rel2
I_in2[np.int(initial_period * 1 / dt):, :] = np.ones((stim_time, 1)) @ I_rel2
stim_num = I_in2.shape[1]
# simulate response to this current step range
V_m_AUC = STN(V_init * np.ones(stim_num), g_effect, E2, I_in2, dt, currents_included2, stim_time, stim_num,
C, shift2, scale2, b_param_effect, slope_shift2, np.zeros((len(b_param), stim_num)),
np.zeros((len(currents_included), stim_num)), ind_dict, vol)
# save data and analysis for this current step range including AUC
f['data'].require_dataset('V_m_AUC', shape=V_m_AUC.shape, dtype='float64')
f['data']['V_m_AUC'][:] = V_m_AUC
for i in np.array(['spike_times_Finf_AUC', 'ISI_Finf_AUC', 'Freq_Finf_AUC']): # analysis_names:
f['analysis'].require_dataset(i, shape=(I_rel2.shape[1],), dtype=dtyp, compression="gzip",
compression_opts=9)
f['analysis'].require_dataset('F_inf_Finf_AUC', shape=(I_rel2.shape[1],), dtype='float64')
f['analysis'].require_dataset('F_null_Finf_AUC', shape=(I_rel2.shape[1],), dtype='float64')
for stim in range(I_rel2.shape[1]):
# spike peak detection
ap_all, ap_prop = scipy.signal.find_peaks(V_m_AUC[stim, stim_start:], min_spike_height,
prominence=prominence, distance=1 * 1 / dt)
f['analysis']['spike_times_Finf_AUC'][stim] = ap_all * dt
f['analysis']['ISI_Finf_AUC'][stim] = np.array(
[x - f['analysis']['spike_times_Finf_AUC'][stim][i - 1] if i else None for i, x in
enumerate(f['analysis']['spike_times_Finf_AUC'][stim])][1:]) # msec
f['analysis']['Freq_Finf_AUC'][stim] = 1 / (f['analysis']['ISI_Finf_AUC'][stim] / 1000)
if f['analysis']['Freq_Finf_AUC'][stim].size != 0:
last_second_spikes = f['analysis']['spike_times_Finf_AUC'][stim] > np.int(
f['data'].attrs['stim_time'] * dt - (500 + f['data'].attrs['initial_period']))
if np.any(last_second_spikes == True):
last_second_spikes_1 = np.delete(last_second_spikes,
0) # remove first spike because it has no ISI/freq
f['analysis']['F_inf_Finf_AUC'][stim] = np.mean(
f['analysis']['Freq_Finf_AUC'][stim][last_second_spikes_1])
else:
f['analysis']['F_null_Finf_AUC'][stim] = 0
f['analysis']['F_inf_Finf_AUC'][stim] = 0
f['analysis'].require_dataset('AUC', shape=(1,), dtype='float64', compression="gzip",
compression_opts=9)
f['analysis']['AUC'][:] = np.trapz(f['analysis']['F_inf_Finf_AUC'][:], I_rel2[:] * 1000, dI * 1000)
except:
f['analysis'].require_dataset('AUC', shape=(1,), dtype='float64', compression="gzip",
compression_opts=9)
f['analysis']['AUC'][:] = 0
# ramp protocol #####################
# setup ramp current
sec = 4
ramp_len = np.int(sec *1000 *1/dt)
stim_time = ramp_len *2
I_amp = np.array([0, -0.000001])
I_amp = np.reshape(I_amp, (1,I_amp.shape[0]))
I_step = np.zeros((stim_time + np.int(initial_period * 1 / dt), 1))@ I_amp
I_step[np.int(initial_period * 1 / dt):, :] = np.ones((stim_time, 1)) @ I_amp
stim_num_step = I_step.shape[1]
start = np.int(initial_period * 1 / dt)
I_step[start:np.int(start+ramp_len),0] = np.linspace(0,high, ramp_len)
I_step[np.int(start + ramp_len):np.int(start + ramp_len*2),0] = np.linspace(high, 0, ramp_len)
stim_time_step = I_step.shape[0]
# simulate response to ramp
V_m_step = STN(-70 * np.ones(stim_num_step), g_effect, E2, I_step, dt, currents_included2, stim_time_step,
stim_num_step, C, shift2, scale2, b_param_effect, slope_shift2,
np.zeros((len(b_param), stim_num_step)), np.zeros((len(currents_included), stim_num_step)), ind_dict,
vol)
# spike peak detection
ap_all, ap_prop = scipy.signal.find_peaks(V_m_step[0, stim_start:], min_spike_height,
prominence=prominence, distance=1 * 1 / dt)
# save ramp firing properties
with h5py.File(fname, "a") as f:
f['data'].create_dataset('V_m_ramp', data=V_m_step, dtype='float64', compression="gzip", compression_opts=9)
if ap_all.shape[0] != 0:
if ap_all[0] < stim_start + ramp_len:
f['analysis'].create_dataset('ramp_I_up', data =I_step[stim_start+ap_all[0], 0] *1000, dtype='float64')
else:
f['analysis'].create_dataset('ramp_I_up', data =np.NaN, dtype='float64')
if ap_all[-1] > stim_start + ramp_len:
f['analysis'].create_dataset('ramp_I_down', data =I_step[stim_start+ap_all[-1], 0]*1000, dtype='float64')
else:
f['analysis'].create_dataset('ramp_I_down', data =np.NaN, dtype='float64')
f['analysis'].create_dataset('hysteresis', data =f['analysis']['ramp_I_down'][()] - f['analysis']['ramp_I_up'][()], dtype='float64')
else: # if no spikes in response to ramp
f['analysis'].create_dataset('ramp_I_up', data =np.NaN, dtype='float64')
f['analysis'].create_dataset('ramp_I_down', data =np.NaN, dtype='float64')
f['analysis'].create_dataset('hysteresis', data =np.NaN, dtype='float64')
if f.__bool__():
f.close()
gc.collect()
def SA_STN(V_init, g, E, I_in, dt, currents_included, stim_time, stim_num, C, shift, scale, b_param, slope_shift,
gating, current, prominence, desired_AUC_width, folder, high, low, number_steps, initial_period, sec, vol,
lin_array, log_array, alt_types, alt_ind, alt):
""" Sensitivity analysis for STN model
Parameters
----------
V_init : float
Initial membrane potential
g : dict
maximal conductance for currents
E : dict
Reversal potentials for currents
I_in : array
Current stimulus input to model
dt : float
fixed time step
currents_included : dict
Boolean to include current in simulation or not
stim_time :
length of stimulus ()
stim_num : int
number of different stimulus currents
C : float
membrane capacitance
shift : dict
for each gating parameter
scale : dict
for each gating parameter
b_param : dict
Boltzmann parameter values for each gating parameter
slope_shift : dict
for each gating parameter
gating : array
array of gating values over time
current : array
array of current values over time
prominence: float
required spike prominence
desired_AUC_width: float
width of AUC from rheobase to rheobase + desired_AUC_width
folder : str
folder to save hdf5 in
high : float
upper current step magnitude
low : float
lowest current step magnitude
number_steps : int
number of current steps between low and high
initial_period: float
length of I = 0 before current step
sec : float
length of current step in seconds
vol : float
cell volume
lin_array : array
array of linear shifts from -10 tp 10
log_array : array
array of log2 values from -1 to 1
alt_types : array
array of arrays of strings of [variable, alteration type]
alt_ind : int
index for alt_type
alt :
index for lin_array of log_array
Returns
-------
saves hdf5 file to folder
"""
# setup for simulation
g_multi = copy.deepcopy(g)
E_multi = copy.deepcopy(E)
b_param_multi = copy.deepcopy(b_param)
shift_multi = copy.deepcopy(shift)
scale_multi = copy.deepcopy(scale)
slope_shift_multi = copy.deepcopy(slope_shift)
currents_included_multi = copy.deepcopy(currents_included)
gating2 = copy.deepcopy(gating)
current2 = copy.deepcopy(current)
g2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in g_multi.items():
g2[k] = v
b_param2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64[:], )
for k, v in b_param_multi.items():
b_param2[k] = np.array(v)
shift2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in shift_multi.items():
shift2[k] = v
scale2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in scale_multi.items():
scale2[k] = v
slope_shift2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in slope_shift_multi.items():
slope_shift2[k] = v
E2 = Dict.empty(key_type=types.unicode_type, value_type=types.float64, )
for k, v in E_multi.items():
E2[k] = v
currents_included2 = Dict.empty(key_type=types.unicode_type, value_type=types.boolean, )
for k, v in currents_included_multi.items():
currents_included2[k] = v
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(
['m', 'h', 'n', 'a', 'b', 'c', 'a_mut', 'b_mut', 'd1', 'd2', 'p', 'q', 'r', 'Ca_conc', 'E_Ca']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'A_mut', 'L', 'T', 'Leak', 'Ca_K']):
ind_dict[var] = i
i += 1
#### alteration #######################################################################
if alt_types[alt_ind, 1] == 'shift':
shift2[alt_types[alt_ind, 0]] = lin_array[alt]
fname = os.path.join(folder,
"{}_{}_{}.hdf5".format(alt_types[alt_ind, 1], alt_types[alt_ind, 0], lin_array[alt]))
elif alt_types[alt_ind, 1] == 'slope':
b_param2[alt_types[alt_ind, 0]][1] = b_param2[alt_types[alt_ind, 0]][1] * log_array[alt]
fname = os.path.join(folder,
"{}_{}_{}.hdf5".format(alt_types[alt_ind, 1], alt_types[alt_ind, 0],
np.int(np.around(log_array[alt], 2) * 100)))
if b_param2[alt_types[alt_ind, 0]][2] != 1:
# adjust slope_shift to account for slope not changing (rotating) around 0.5
V_test = np.arange(-150, 100, 0.001)
if b_param2[alt_types[alt_ind, 0]].shape[0] == 4:
steady_state = ((1 - b_param2[alt_types[alt_ind, 0]][3]) / (
1 + np.exp(
(V_test - b_param2[alt_types[alt_ind, 0]][0]) / b_param2[alt_types[alt_ind, 0]][1])) +
b_param2[alt_types[alt_ind, 0]][3]) ** b_param2[alt_types[alt_ind, 0]][2]
orig = ((1 - b_param[alt_types[alt_ind, 0]][3]) / (
1 + np.exp((V_test - b_param[alt_types[alt_ind, 0]][0]) / b_param[alt_types[alt_ind, 0]][1])) +
b_param[alt_types[alt_ind, 0]][3]) ** b_param[alt_types[alt_ind, 0]][2]
else:
steady_state = (1 / (1 + np.exp(
(V_test - b_param2[alt_types[alt_ind, 0]][0]) / b_param2[alt_types[alt_ind, 0]][1]))) ** \
b_param2[alt_types[alt_ind, 0]][2]
orig = (1 / (1 + np.exp(
(V_test - b_param[alt_types[alt_ind, 0]][0]) / b_param[alt_types[alt_ind, 0]][1]))) ** \
b_param[alt_types[alt_ind, 0]][2]
orig_V_half = V_test[(np.abs(orig - 0.5)).argmin()]
V_half_new = V_test[(np.abs(steady_state - 0.5)).argmin()]
slope_shift2[alt_types[alt_ind, 0]] = orig_V_half - V_half_new
elif alt_types[alt_ind, 1] == 'g':
g2[alt_types[alt_ind, 0]] = g2[alt_types[alt_ind, 0]] * log_array[alt]
fname = os.path.join(folder, "{}_{}_{}.hdf5".format(alt_types[alt_ind, 1], alt_types[alt_ind, 0],
np.int(np.around(log_array[alt], 2) * 100)))
###########################################################################
# initial simulation for current range from low to high
V_m = STN(V_init * np.ones((stim_num)), g2, E2, I_in, dt, currents_included2, stim_time, stim_num, C, shift2,
scale2, b_param2, slope_shift2, gating2, current2, ind_dict, vol)
stim_start = np.int(initial_period * 1 / dt)
min_spike_height = V_m[0, stim_start] + prominence
# create hdf5 file and save data and metadata
with h5py.File(fname, "a") as f:
data = f.create_group("data")
data.create_dataset('V_m', data=V_m, dtype='float64', compression="gzip", compression_opts=9)
metadata = {'I_high': high, 'I_low': low, 'stim_num': number_steps, 'stim_time': stim_time, 'sec': sec, 'dt': dt,
'initial_period': initial_period, 'g': json.dumps(dict(g2)), 'E': json.dumps(dict(E2)), 'C': C,
'V_init': V_init, 'vol': vol, 'ind_dict': json.dumps(dict(ind_dict)),
'b_param': json.dumps(dict(b_param2), cls=NumpyEncoder),
'currents': json.dumps(dict(currents_included2)), 'shift': json.dumps(dict(shift2)),
'var_change': str(alt_types[alt_ind, 1]), 'var': str(alt_types[alt_ind, 0]),
'scale': json.dumps(dict(scale2)), 'slope_shift': json.dumps(dict(slope_shift2)),
'prominence': prominence, 'desired_AUC_width': desired_AUC_width,
'alteration_info': str(alt_types[alt_ind, :]), 'alt_index': str(alt),
'min_spike_height': min_spike_height,}
if alt_types[alt_ind, 1] == 'shift':
meta2 = {'alteration': lin_array[alt]}
else:
meta2 = {'alteration': log_array[alt]}
data.attrs.update(meta2)
# firing frequency analysis and rheobase
with h5py.File(fname, "r+") as f:
I_mag = np.arange(low, high, (high - low) / stim_num)
analysis = f.create_group("analysis")
dtyp = h5py.special_dtype(vlen=np.dtype('float64'))
for i in np.array(['spike_times', 'spike_times_rheo', 'ISI', 'Freq']):
analysis.create_dataset(i, (I_mag.shape[0],), dtype=dtyp, compression="gzip", compression_opts=9)
analysis.create_dataset('F_inf', (I_mag.shape[0],), dtype='float64')
analysis.create_dataset('F_null', (I_mag.shape[0],), dtype='float64')
analysis.create_dataset('AP', (I_mag.shape[0],), dtype='bool', compression="gzip", compression_opts=9)
analysis.create_dataset('AP_rheo', (I_mag.shape[0],), dtype='bool', compression="gzip", compression_opts=9)
for i in np.array(['rheobase']):
analysis.create_dataset(i, (1,), dtype='float64', compression="gzip", compression_opts=9)
try:
for stim in range(I_mag.shape[0]):
# spike detection
ap_all, ap_prop = scipy.signal.find_peaks(V_m[stim, stim_start:], min_spike_height,
prominence=prominence, distance=1 * 1 / dt)
analysis['spike_times'][stim] = ap_all * dt
if analysis['spike_times'][stim].size == 0:
analysis['AP'][stim] = False
else:
analysis['AP'][stim] = True
analysis['ISI'][stim] = np.array([x - analysis['spike_times'][stim][i - 1] if i else None for i, x in
enumerate(analysis['spike_times'][stim])][1:]) # msec
analysis['Freq'][stim] = 1 / (analysis['ISI'][stim] / 1000)
if analysis['Freq'][stim].size != 0:
analysis['F_null'][stim] = analysis['Freq'][stim][0]
if np.argwhere(ap_all >= 1000 * 1 / dt).shape[0] != 0: # if first spike is after 1 ms
half_second_spikes = np.logical_and(analysis['spike_times'][stim] < np.int(
ap_all[np.argwhere(ap_all >= 1000 * 1 / dt)[0][0]] * dt + 500),
analysis['spike_times'][stim] > np.int(
ap_all[np.argwhere(ap_all >= 1000 * 1 / dt)[0][0]] * dt)) # spikes in half second window
if np.sum(half_second_spikes == True) > 1:
half_second_spikes_1 = np.delete(half_second_spikes,
np.argwhere(half_second_spikes == True)[0][0])
analysis['F_inf'][stim] = np.mean(analysis['Freq'][stim][half_second_spikes_1])
else:
analysis['F_inf'][stim] = 0
else:
analysis['F_inf'][stim] = 0
else: # if no ISI detected
analysis['F_null'][stim] = 0
analysis['F_inf'][stim] = 0
except:
print(alt_types[alt_ind, 0], alt_types[alt_ind, 1], alt, 'Freq analysis unsuccessful')
# rheobase
# find where the first AP occurs
if np.any(analysis['AP'][:] == True):
# setup current range to be between the current step before and with first AP
I_first_spike = I_mag[np.argwhere(analysis['AP'][:] == True)[0][0]]
I_before_spike = I_mag[np.argwhere(analysis['AP'][:] == True)[0][0] - 1]
analysis.create_dataset('I_first_spike', data=I_first_spike, dtype='float64')
analysis.create_dataset('I_before_spike', data=I_before_spike, dtype='float64')
number_steps_rheo = 100
stim_time, I_in_rheo, stim_num_rheo, V_m_rheo = stimulus_init(I_before_spike, I_first_spike,
number_steps_rheo,
initial_period, dt, sec)
# simulate response to this current range and save
V_m_rheo = STN(V_init * np.ones((stim_num_rheo)), g2, E2, I_in_rheo, dt, currents_included2, stim_time,
stim_num_rheo, C, shift2, scale2, b_param2, slope_shift2,
np.zeros((len(b_param), stim_num_rheo)), np.zeros((len(currents_included), stim_num_rheo)),
ind_dict, vol)
try:
# run AP detection
for stim in range(I_in_rheo.shape[1]): # range(stim_num):
ap_all, ap_prop = scipy.signal.find_peaks(V_m_rheo[stim, stim_start:], min_spike_height,
prominence=prominence, distance=1 * 1 / dt)
analysis['spike_times_rheo'][stim] = ap_all * dt
if analysis['spike_times_rheo'][stim].size == 0:
analysis['AP_rheo'][stim] = False
else:
analysis['AP_rheo'][stim] = True
# find rheobase as smallest current step to elicit AP and save
analysis['rheobase'][:] = I_in_rheo[-1, np.argwhere(analysis['AP_rheo'][:] == True)[0][0]] * 1000
except:
print(alt_types[alt_ind, 0], alt_types[alt_ind, 1], alt, 'RHEOBASE NO AP')
# AUC
with h5py.File(fname, "a") as f:
# find rheobase of steady state firing (F_inf)
F_inf = f['analysis']['F_inf'][:]
try:
# setup current range to be in current step between no and start of steady state firing
F_inf_start_ind = np.argwhere(F_inf != 0)[0][0]
dI = (high- low) / stim_num
I_start = I_mag[F_inf_start_ind - 1]
I_end = I_mag[F_inf_start_ind]
I_rel = np.arange(I_start, I_end + (I_end - I_start) / 100, (I_end - I_start) / 100)
I_rel = np.reshape(I_rel, (1, I_rel.shape[0]))
I_in = np.zeros((stim_time + np.int(initial_period * 1 / dt), 1)) @ I_rel
I_in[np.int(initial_period * 1 / dt):, :] = np.ones((stim_time, 1)) @ I_rel
stim_num = I_in.shape[1]
# simulate response to this current step range
V_m_inf_rheo = STN(V_init * np.ones(stim_num), g2, E2, I_in, dt, currents_included2, stim_time, stim_num, C,
shift2, scale2, b_param2, slope_shift2, np.zeros((len(b_param), stim_num)),
np.zeros((len(currents_included), stim_num)), ind_dict, vol)
# save response, analyse and save
for i in np.array(['spike_times_Finf', 'ISI_Finf', 'Freq_Finf']):
f['analysis'].require_dataset(i, shape=(I_rel.shape[1],), dtype=dtyp, compression="gzip",
compression_opts=9)
f['analysis'].require_dataset('F_inf_Finf', shape=(I_rel.shape[1],), dtype='float64')
f['analysis'].require_dataset('F_null_Finf', shape=(I_rel.shape[1],), dtype='float64')
for stim in range(I_rel.shape[1]):
# spike peak detection
ap_all, ap_prop = scipy.signal.find_peaks(V_m_inf_rheo[stim, stim_start:], min_spike_height,
prominence=prominence, distance=1 * 1 / dt)
f['analysis']['spike_times_Finf'][stim] = ap_all * dt
f['analysis']['ISI_Finf'][stim] = np.array(
[x - f['analysis']['spike_times_Finf'][stim][i - 1] if i else None for i, x in
enumerate(f['analysis']['spike_times_Finf'][stim])][1:]) # msec
f['analysis']['Freq_Finf'][stim] = 1 / (f['analysis']['ISI_Finf'][stim] / 1000)
if f['analysis']['Freq_Finf'][stim].size != 0:
if np.argwhere(ap_all >= 1000 * 1 / dt).shape[0] != 0: # if first spike is after 1 ms
half_second_spikes = np.logical_and(f['analysis']['spike_times_Finf'][stim] < np.int(
ap_all[np.argwhere(ap_all >= 1000 * 1 / dt)[0][0]] * dt + 500),
f['analysis']['spike_times_Finf'][stim] > np.int(
ap_all[np.argwhere(ap_all >= 1000 * 1 / dt)[0][0]] * dt)) # spikes in half second window
if np.sum(half_second_spikes == True) > 1:
half_second_spikes_1 = np.delete(half_second_spikes,
np.argwhere(half_second_spikes == True)[0][0]) # remove first spike because it has no ISI/freq
f['analysis']['F_inf_Finf'][stim] = np.mean(f['analysis']['Freq_Finf'][stim][half_second_spikes_1])
else:
f['analysis']['F_inf_Finf'][stim] = 0
else:
f['analysis']['F_null_Finf'][stim] = 0
f['analysis']['F_inf_Finf'][stim] = 0
else:
f['analysis']['F_null_Finf'][stim] = 0
f['analysis']['F_inf_Finf'][stim] = 0
# find AUC relative to F_inf rheobase
# setup current inputs with current range around F_inf rheobase
I_first_inf = I_rel[0, np.argwhere(f['analysis']['F_inf_Finf'][:] != 0)[0][0]]
I_start = I_first_inf - 0.00005
I_end = I_first_inf + 0.2 * high
I_rel2 = np.arange(I_start, I_end + (I_end - I_start) / 100, (I_end - I_start) / 100)
I_rel2 = np.reshape(I_rel2, (1, I_rel2.shape[0]))
I_in2 = np.zeros((stim_time + np.int(initial_period * 1 / dt), 1)) @ I_rel2
I_in2[np.int(initial_period * 1 / dt):, :] = np.ones((stim_time, 1)) @ I_rel2
stim_num = I_in2.shape[1]
# simulate response to this current step range
V_m_AUC = STN(V_init * np.ones(stim_num), g2, E2, I_in2, dt, currents_included2, stim_time, stim_num, C,
shift2, scale2, b_param2, slope_shift2, np.zeros((len(b_param), stim_num)),
np.zeros((len(currents_included), stim_num)), ind_dict, vol)
# save data and analysis for this current step range including AUC
f['data'].require_dataset('V_m_AUC', shape=V_m_AUC.shape, dtype='float64')
f['data']['V_m_AUC'][:] = V_m_AUC
dtyp = h5py.special_dtype(vlen=np.dtype('float64'))
for i in np.array(['spike_times_Finf_AUC', 'ISI_Finf_AUC', 'Freq_Finf_AUC']): # analysis_names:
f['analysis'].require_dataset(i, shape=(I_rel2.shape[1],), dtype=dtyp, compression="gzip",
compression_opts=9)
f['analysis'].require_dataset('F_inf_Finf_AUC', shape=(I_rel2.shape[1],), dtype='float64')
f['analysis'].require_dataset('F_null_Finf_AUC', shape=(I_rel2.shape[1],), dtype='float64')
for stim in range(I_rel2.shape[1]):
# spike peak detection
ap_all, ap_prop = scipy.signal.find_peaks(V_m_AUC[stim, stim_start:], min_spike_height,
prominence=prominence, distance=1 * 1 / dt)
f['analysis']['spike_times_Finf_AUC'][stim] = ap_all * dt
f['analysis']['ISI_Finf_AUC'][stim] = np.array(
[x - f['analysis']['spike_times_Finf_AUC'][stim][i - 1] if i else None for i, x in
enumerate(f['analysis']['spike_times_Finf_AUC'][stim])][1:]) # msec
f['analysis']['Freq_Finf_AUC'][stim] = 1 / (f['analysis']['ISI_Finf_AUC'][stim] / 1000)
if f['analysis']['Freq_Finf_AUC'][stim].size != 0:
if np.argwhere(ap_all >= 1000 * 1 / dt).shape[0] != 0: # if first spike is after 1 ms
half_second_spikes = np.logical_and(f['analysis']['spike_times_Finf_AUC'][stim] < np.int(
ap_all[np.argwhere(ap_all >= 1000 * 1 / dt)[0][0]] * dt + 500),
f['analysis']['spike_times_Finf_AUC'][stim] > np.int(
ap_all[np.argwhere(ap_all >= 1000 * 1 / dt)[0][0]] * dt)) # spikes in half second window
if np.sum(half_second_spikes == True) > 1:
half_second_spikes_1 = np.delete(half_second_spikes,
np.argwhere(half_second_spikes == True)[0][
0]) # remove first spike because it has no ISI/freq
f['analysis']['F_inf_Finf_AUC'][stim] = np.mean(
f['analysis']['Freq_Finf_AUC'][stim][half_second_spikes_1])
else:
f['analysis']['F_inf_Finf_AUC'][stim] = 0
else:
f['analysis']['F_null_Finf_AUC'][stim] = 0
f['analysis']['F_inf_Finf_AUC'][stim] = 0
else:
f['analysis']['F_null_Finf_AUC'][stim] = 0
f['analysis']['F_inf_Finf_AUC'][stim] = 0
f['analysis'].require_dataset('AUC', shape=(1,), dtype='float64', compression="gzip",
compression_opts=9)
f['analysis']['AUC'][:] = np.trapz(f['analysis']['F_inf_Finf_AUC'][:],
I_rel2[:] * 1000, dI * 1000)
except:
f['analysis'].require_dataset('AUC', shape=(1,), dtype='float64', compression="gzip",
compression_opts=9)
f['analysis']['AUC'][:] = 0
if f.__bool__():
f.close()
gc.collect()