""" 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()