Functions and plotting scripts for each model, sensitivity analysis in each model and KCNA1 mutations in each model

This commit is contained in:
nkoch1 2022-08-28 15:57:33 -04:00
parent 64e663d857
commit 47f5c1db02
42 changed files with 8963 additions and 11 deletions

View File

@ -0,0 +1,960 @@
import h5py
import gc
import scipy
import json
import copy
import os
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 Cb_stellate(V, g, E, I_in, dt, currents_included, stim_time, stim_num, C, shift, scale, b_param, slope_shift,
gating, current, ind_dict):
""" Simulate current and Vm for Cb stellate cell model from Alexander et al 2019 (https://dx.doi.org/10.1523%2FENEURO.0126-19.2019)
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
Returns
-------
V_m : array
Simulated membrane potential in response to I_in
"""
V_m = np.zeros((stim_num, stim_time))
# 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)) * 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)) * 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)) * n_inf
# A-type K activation
n_A_inf = (1 / (1 + np.exp((V - shift["n_A"] - slope_shift["n_A"] - b_param["n_A"][0]) / b_param["n_A"][1]))) ** \
b_param["n_A"][2]
gating[ind_dict['n_A'], :] = np.ones((stim_num)) * n_A_inf
# A-type K inactivation
h_A_inf = (1 / (1 + np.exp((V - shift["h_A"] - slope_shift["h_A"] - b_param["h_A"][0]) / b_param["h_A"][1]))) ** \
b_param["h_A"][2]
gating[ind_dict['h_A'], :] = np.ones((stim_num)) * h_A_inf
# mutant A-type K activation
n_A_mut_inf = (1 / (1 + np.exp((V - shift["n_A_mut"] - slope_shift["n_A_mut"] - b_param["n_A_mut"][0]) / b_param["n_A_mut"][1]))) ** \
b_param["n_A_mut"][2]
gating[ind_dict['n_A_mut'], :] = np.ones((stim_num)) * n_A_mut_inf
# mutant A-type K inactivation
h_A_mut_inf = (1 / (1 + np.exp((V - shift["h_A_mut"] - slope_shift["h_A_mut"] - b_param["h_A_mut"][0]) / b_param["h_A_mut"][1]))) ** \
b_param["h_A_mut"][2]
gating[ind_dict['h_A_mut'], :] = np.ones((stim_num)) * h_A_mut_inf
# T-type Ca activation
m_T_inf = (1 / (1 + np.exp((V - shift["m_T"] - slope_shift["m_T"] - b_param["m_T"][0]) / b_param["m_T"][1]))) ** \
b_param["m_T"][2]
gating[ind_dict['m_T'], :] = np.ones((stim_num)) * m_T_inf
# T-type Ca inactivation
h_T_inf = (1 / (1 + np.exp((V - shift["h_T"] - slope_shift["h_T"] - b_param["h_T"][0]) / b_param["h_T"][1]))) ** \
b_param["h_T"][2]
gating[ind_dict['h_T'], :] = np.ones((stim_num)) * h_T_inf
# initialize currents and Vm
current[ind_dict['Na'], :] = np.zeros((stim_num))
current[ind_dict['Kd'], :] = np.zeros((stim_num))
current[ind_dict['Leak'], :] = np.zeros((stim_num))
current[ind_dict['A'], :] = np.zeros((stim_num))
current[ind_dict['A_mut'], :] = np.zeros((stim_num))
current[ind_dict['T'], :] = np.zeros((stim_num))
V_m[:, 0] = V
t = 1
while t < stim_time:
# 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]
gating[ind_dict['m'], :] = m_inf
gating[ind_dict['m'], :][gating[ind_dict['m'], :] > 1.] = 1.
gating[ind_dict['m'], :][gating[ind_dict['m'], :] < 0.] = 0.
# 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 = 0.1 + (2 * 322 * 46) / (4 * np.pi * (V_m[:, t - 1] - shift["h"] - -74) ** 2 + (46 ** 2))
h_dot = (h_inf - gating[ind_dict['h'], :]) * (1 / (tau_h * scale["h"]))
gating[ind_dict['h'], :] = gating[ind_dict['h'], :] + h_dot * dt
gating[ind_dict['h'], :][gating[ind_dict['h'], :] > 1.] = 1.
gating[ind_dict['h'], :][gating[ind_dict['h'], :] < 0.] = 0.
# 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 = 6 / (1 + np.exp((V_m[:, t - 1] - shift["n"] + 23) / 15))
n_dot = (n_inf - gating[ind_dict['n'], :]) * (1 / (tau_n * scale["n"]))
gating[ind_dict['n'], :] = gating[ind_dict['n'], :] + n_dot * dt
gating[ind_dict['n'], :][gating[ind_dict['n'], :] > 1.] = 1.
gating[ind_dict['n'], :][gating[ind_dict['n'], :] < 0.] = 0.
# A-type K activation
n_A_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["n_A"] - slope_shift["n_A"] - b_param["n_A"][0]) / b_param["n_A"][1]))) ** \
b_param["n_A"][2]
tau_n_A = 5
n_A_dot = (n_A_inf - gating[ind_dict['n_A'], :]) * (1 / (tau_n_A * scale["n_A"]))
gating[ind_dict['n_A'], :] = gating[ind_dict['n_A'], :] + n_A_dot * dt
gating[ind_dict['n_A'], :][gating[ind_dict['n_A'], :] > 1.] = 1.
gating[ind_dict['n_A'], :][gating[ind_dict['n_A'], :] < 0.] = 0.
# A-type K inactivation
h_A_inf = (1 / (1 + np.exp(
(V_m[:, t - 1] - shift["h_A"] - slope_shift["h_A"] - b_param["h_A"][0]) / b_param["h_A"][1]))) ** \
b_param["h_A"][2]
tau_h_A = 10
h_A_dot = (h_A_inf - gating[ind_dict['h_A'], :]) * (1 / (tau_h_A * scale["h_A"]))
gating[ind_dict['h_A'], :] = gating[ind_dict['h_A'], :] + h_A_dot * dt
gating[ind_dict['h_A'], :][gating[ind_dict['h_A'], :] > 1.] = 1.
gating[ind_dict['h_A'], :][gating[ind_dict['h_A'], :] < 0.] = 0.
# mutant A-type K activation
n_A_mut_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["n_A_mut"] - slope_shift["n_A_mut"] - b_param["n_A_mut"][0]) / b_param["n_A_mut"][1]))) ** \
b_param["n_A_mut"][2]
tau_n_A_mut = 5
n_A_mut_dot = (n_A_mut_inf - gating[ind_dict['n_A_mut'], :]) * (1 / (tau_n_A_mut * scale["n_A_mut"]))
gating[ind_dict['n_A_mut'], :] = gating[ind_dict['n_A_mut'], :] + n_A_mut_dot * dt
gating[ind_dict['n_A_mut'], :][gating[ind_dict['n_A_mut'], :] > 1.] = 1.
gating[ind_dict['n_A_mut'], :][gating[ind_dict['n_A_mut'], :] < 0.] = 0.
# mutant A-type K inactivation
h_A_mut_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["h_A_mut"] - slope_shift["h_A_mut"] - b_param["h_A_mut"][0]) / b_param["h_A_mut"][1]))) ** \
b_param["h_A_mut"][2]
tau_h_A_mut = 10
h_A_mut_dot = (h_A_mut_inf - gating[ind_dict['h_A_mut'], :]) * (1 / (tau_h_A_mut * scale["h_A_mut"]))
gating[ind_dict['h_A_mut'], :] = gating[ind_dict['h_A_mut'], :] + h_A_mut_dot * dt
gating[ind_dict['h_A_mut'], :][gating[ind_dict['h_A_mut'], :] > 1.] = 1.
gating[ind_dict['h_A_mut'], :][gating[ind_dict['h_A_mut'], :] < 0.] = 0.
# T-type Ca activation
m_T_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["m_T"] - slope_shift["m_T"] - b_param["m_T"][0]) / b_param["m_T"][1]))) ** \
b_param["m_T"][2]
gating[ind_dict['m_T'], :] = m_T_inf
gating[ind_dict['m_T'], :][gating[ind_dict['m_T'], :] > 1.] = 1.
gating[ind_dict['m_T'], :][gating[ind_dict['m_T'], :] < 0.] = 0.
# T-type Ca inactivation
h_T_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["h_T"] - slope_shift["h_T"] - b_param["h_T"][0]) / b_param["h_T"][1]))) ** \
b_param["h_T"][2]
tau_h_T = 15
h_T_dot = (h_T_inf - gating[ind_dict['h_T'], :]) * (1 / (tau_h_T * scale["h_T"]))
gating[ind_dict['h_T'], :] = gating[ind_dict['h_T'], :] + h_T_dot * dt
gating[ind_dict['h_T'], :][gating[ind_dict['h_T'], :] > 1.] = 1.
gating[ind_dict['h_T'], :][gating[ind_dict['h_T'], :] < 0.] = 0.
# 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['A'], :] = g["A"] * gating[ind_dict['n_A'], :] * gating[ind_dict['h_A'], :] * (
V_m[:, t - 1] - E["K"]) * currents_included["A"]
current[ind_dict['A_mut'], :] = g["A_mut"] * gating[ind_dict['n_A_mut'], :] * gating[ind_dict['h_A_mut'], :] * (
V_m[:, t - 1] - E["K"]) * currents_included["A_mut"]
current[ind_dict['T'], :] = g["T"] * gating[ind_dict['m_T'], :] * gating[ind_dict['h_T'], :] * (
V_m[:, t - 1] - E["Ca"]) * currents_included["T"]
# 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['A'], :] +
current[ind_dict['A_mut'], :] + current[ind_dict['T'], :]))
# update V_m
V_m[:, t] = V_m[:, t - 1] + V_dot * dt
t += 1
return V_m
def Cb_stellate_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):
""" Simulate mutations for Cb Stellate 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
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', 'n_A', 'h_A', 'n_A_mut', 'h_A_mut', 'm_T', 'h_T']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'A_mut', 'T', 'Leak']):
ind_dict[var] = i
i += 1
# effect of mutation
g_effect['A_mut'] = g['A_mut'] * mutations[mut]['g_ratio']
b_param_effect['n_A_mut'][0] = b_param_effect['n_A_mut'][0] + mutations[mut]['activation_Vhalf_diff']
b_param_effect['n_A_mut'][1] = b_param_effect['n_A_mut'][1] * mutations[mut]['activation_k_ratio']
# initial simulation for current range from low to high
V_m = Cb_stellate(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)
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, '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', 'spike_times_rheo', '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('AP_rheo', (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 = Cb_stellate(V_init * np.ones((stim_num_rheo)), g_effect, E2, I_in_rheo, dt, currents_included2,
stim_time, stim_num_rheo, 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'].create_dataset('V_m_rheo', data=V_m_rheo, dtype='float64', compression="gzip", compression_opts=9)
# run AP detection
stim_start = np.int(initial_period * 1 / dt)
for stim in range(I_in_rheo.shape[1]):
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
# find rheobase of steady state firing (F_inf)
with h5py.File(fname, "a") as f:
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 = Cb_stellate(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)
# 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:
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 = Cb_stellate(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)
# 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:
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_rel_acc', shape=(1,), dtype='float64', compression="gzip",
compression_opts=9)
f['analysis']['AUC_rel_acc'][:] = np.trapz(f['analysis']['F_inf_Finf_AUC'][:],
I_rel2[:] * 1000, dI * 1000) # AUC
except:
f['analysis'].require_dataset('AUC_rel_acc', shape=(1,), dtype='float64', compression="gzip",
compression_opts=9)
f['analysis']['AUC_rel_acc'][:] = 0 # AUC
print('exception', f['analysis']['AUC_rel_acc'][()])
# 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 = Cb_stellate(-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)
# 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_Cb_stellate(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, lin_array, log_array, alt_types, alt_ind, alt):
""" Sensitivity analysis for Cb Stellate 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
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', 'n_A', 'h_A', 'n_A_mut', 'h_A_mut', 'm_T', 'h_T']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'A_mut', 'T', 'Leak']):
ind_dict[var] = i
i += 1
# alteration to gating ##########################################################################
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 = Cb_stellate(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)
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, '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, }
data.attrs.update(metadata)
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 = Cb_stellate(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)
# run AP detection
try:
# stim_start = np.int(initial_period * 1 / dt)
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
except:
print(alt_types[alt_ind, 0], alt_types[alt_ind, 1], alt, 'RHEOBASE NO AP')
# 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
# find rheobase of steady state firing (F_inf)
with h5py.File(fname, "a") as f:
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 = Cb_stellate(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)
# 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 = Cb_stellate(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)
# 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']['AUCc'][:] = 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()

View File

@ -0,0 +1,993 @@
import h5py
import scipy
import json
import copy
import os
import gc
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 Cb_stellate_Kv(V, g, E, I_in, dt, currents_included, stim_time, stim_num, C, shift, scale, b_param,
slope_shift, gating, current, ind_dict):
""" Simulate current and Vm for Cb stellate cell model from Alexander et al 2019 (https://dx.doi.org/10.1523%2FENEURO.0126-19.2019)
with Kv1.1 from Ranjan et al. 2019 (https://dx.doi.org/10.3389/fncel.2019.00358) added 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
Returns
-------
V_m : array
Simulated membrane potential in response to I_in
"""
V_m = np.zeros((stim_num, stim_time))
# 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)) * 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)) * 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)) * n_inf
# A-type K activation
n_A_inf = (1 / (1 + np.exp((V - shift["n_A"] - slope_shift["n_A"] - b_param["n_A"][0]) / b_param["n_A"][1]))) ** \
b_param["n_A"][2]
gating[ind_dict['n_A'], :] = np.ones((stim_num)) * n_A_inf
# A-type K inactivation
h_A_inf = (1 / (1 + np.exp((V - shift["h_A"] - slope_shift["h_A"] - b_param["h_A"][0]) / b_param["h_A"][1]))) ** \
b_param["h_A"][2]
gating[ind_dict['h_A'], :] = np.ones((stim_num)) * h_A_inf
# T-type Ca activation
m_T_inf = (1 / (1 + np.exp((V - shift["m_T"] - slope_shift["m_T"] - b_param["m_T"][0]) / b_param["m_T"][1]))) ** \
b_param["m_T"][2]
gating[ind_dict['m_T'], :] = np.ones((stim_num)) * m_T_inf
# T-type Ca inactivation
h_T_inf = (1 / (1 + np.exp((V - shift["h_T"] - slope_shift["h_T"] - b_param["h_T"][0]) / b_param["h_T"][1]))) ** \
b_param["h_T"][2]
gating[ind_dict['h_T'], :] = np.ones((stim_num)) * h_T_inf
# Kv1.1 activation
s_inf = (1 / (1 + np.exp((V - shift["s"] - slope_shift["s"] - b_param["s"][0]) / b_param["s"][1]))) ** \
b_param["s"][2]
gating[ind_dict['s'], :] = np.ones((stim_num)) * s_inf
# Kv1.1 inactivation
u_inf = ((1 - b_param["u"][3]) / (1 + np.exp((V - shift["u"] - slope_shift["u"] - b_param["u"][0])
/ b_param["u"][1])) + b_param["u"][3]) ** b_param["u"][2]
gating[ind_dict['u'], :] = np.ones((stim_num)) * u_inf
# mutant Kv1.1 activation
s_mut_inf = (1 / (1 + np.exp((V - shift["s_mut"] - slope_shift["s_mut"] - b_param["s_mut"][0]) / b_param["s_mut"][1]))) ** \
b_param["s_mut"][2]
gating[ind_dict['s_mut'], :] = np.ones((stim_num)) * s_mut_inf
# mutant Kv1.1 inactivation
u_mut_inf = ((1 - b_param["u_mut"][3]) / (1 + np.exp((V - shift["u"] - slope_shift["u_mut"] - b_param["u_mut"][0])
/ b_param["u_mut"][1])) + b_param["u_mut"][3]) ** b_param["u_mut"][2]
gating[ind_dict['u_mut'], :] = np.ones((stim_num)) * u_mut_inf
# initialize currents and Vm
current[ind_dict['Na'], :] = np.zeros((stim_num))
current[ind_dict['Kd'], :] = np.zeros((stim_num))
current[ind_dict['Kv'], :] = np.zeros((stim_num))
current[ind_dict['Kv_mut'], :] = np.zeros((stim_num))
current[ind_dict['Leak'], :] = np.zeros((stim_num))
current[ind_dict['A'], :] = np.zeros((stim_num))
current[ind_dict['T'], :] = np.zeros((stim_num))
V_m[:, 0] = V
t = 1
while t < stim_time:
# 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]
gating[ind_dict['m'], :] = m_inf
gating[ind_dict['m'], :][gating[ind_dict['m'], :] > 1.] = 1.
gating[ind_dict['m'], :][gating[ind_dict['m'], :] < 0.] = 0.
# 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 = 0.1 + (2 * 322 * 46) / (4 * np.pi * (V_m[:, t - 1] - shift["h"] - -74) ** 2 + (46 ** 2))
h_dot = (h_inf - gating[ind_dict['h'], :]) * (1 / (tau_h * scale["h"]))
gating[ind_dict['h'], :] = gating[ind_dict['h'], :] + h_dot * dt
gating[ind_dict['h'], :][gating[ind_dict['h'], :] > 1.] = 1.
gating[ind_dict['h'], :][gating[ind_dict['h'], :] < 0.] = 0.
# 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 = 6 / (1 + np.exp((V_m[:, t - 1] - shift["n"] + 23) / 15))
n_dot = (n_inf - gating[ind_dict['n'], :]) * (1 / (tau_n * scale["n"]))
gating[ind_dict['n'], :] = gating[ind_dict['n'], :] + n_dot * dt
gating[ind_dict['n'], :][gating[ind_dict['n'], :] > 1.] = 1.
gating[ind_dict['n'], :][gating[ind_dict['n'], :] < 0.] = 0.
# A-type K activation
n_A_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["n_A"] - slope_shift["n_A"] - b_param["n_A"][0]) / b_param["n_A"][1]))) ** \
b_param["n_A"][2]
tau_n_A = 5
n_A_dot = (n_A_inf - gating[ind_dict['n_A'], :]) * (1 / (tau_n_A * scale["n_A"]))
gating[ind_dict['n_A'], :] = gating[ind_dict['n_A'], :] + n_A_dot * dt
gating[ind_dict['n_A'], :][gating[ind_dict['n_A'], :] > 1.] = 1.
gating[ind_dict['n_A'], :][gating[ind_dict['n_A'], :] < 0.] = 0.
# A-type K inactivation
h_A_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["h_A"] - slope_shift["h_A"] - b_param["h_A"][0]) / b_param["h_A"][1]))) ** \
b_param["h_A"][2]
tau_h_A = 10
h_A_dot = (h_A_inf - gating[ind_dict['h_A'], :]) * (1 / (tau_h_A * scale["h_A"]))
gating[ind_dict['h_A'], :] = gating[ind_dict['h_A'], :] + h_A_dot * dt
gating[ind_dict['h_A'], :][gating[ind_dict['h_A'], :] > 1.] = 1.
gating[ind_dict['h_A'], :][gating[ind_dict['h_A'], :] < 0.] = 0.
# T-type Ca activation
m_T_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["m_T"] - slope_shift["m_T"] - b_param["m_T"][0]) / b_param["m_T"][1]))) ** \
b_param["m_T"][2]
gating[ind_dict['m_T'], :] = m_T_inf
gating[ind_dict['m_T'], :][gating[ind_dict['m_T'], :] > 1.] = 1.
gating[ind_dict['m_T'], :][gating[ind_dict['m_T'], :] < 0.] = 0.
# T-type Ca inactivation
h_T_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["h_T"] - slope_shift["h_T"] - b_param["h_T"][0]) / b_param["h_T"][1]))) ** \
b_param["h_T"][2]
tau_h_T = 15
h_T_dot = (h_T_inf - gating[ind_dict['h_T'], :]) * (1 / (tau_h_T * scale["h_T"]))
gating[ind_dict['h_T'], :] = gating[ind_dict['h_T'], :] + h_T_dot * dt
gating[ind_dict['h_T'], :][gating[ind_dict['h_T'], :] > 1.] = 1.
gating[ind_dict['h_T'], :][gating[ind_dict['h_T'], :] < 0.] = 0.
# Kv1.1 activation
s_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["s"] - slope_shift["s"] - b_param["s"][0]) / b_param["s"][1]))) ** \
b_param["s"][2]
temp = 36
vBreak = -46.7
amp1 = 52.7
amp2 = 15.98
vh1 = -49.87
vh2 = -41.64
offset = 0.9
slope1 = 5
slope2 = 24.99
sigswitch = 1 / (1 + np.exp((V_m[:, t - 1] - shift["s"] - vBreak) / 3))
sig1 = sigswitch * amp1 / (1 + np.exp((V_m[:, t - 1] - shift["s"] - vh1) / slope1))
sig2 = (1 - sigswitch) * offset + (amp2 - offset) / (1 + np.exp((V_m[:, t - 1] - shift["s"] - vh2) / slope2))
sTauFunc = sig1 + sig2
s_tau_Q10 = (7.54 * np.exp(-(V_m[:, t - 1] - shift["s"]) / 379.7) * np.exp(-temp / 35.66)) ** ((temp - 25) / 10)
tau_s = sTauFunc / s_tau_Q10
s_dot = (s_inf - gating[ind_dict['s'], :]) * (1 / (tau_s * scale["s"]))
gating[ind_dict['s'], :] = gating[ind_dict['s'], :] + s_dot * dt
gating[ind_dict['s'], :][gating[ind_dict['s'], :] > 1.] = 1.
gating[ind_dict['s'], :][gating[ind_dict['s'], :] < 0.] = 0.
# Kv1.1 inactivation
u_inf = ((1 - b_param["u"][3]) / (1 + np.exp((V_m[:, t - 1] - shift["u"] - slope_shift["u"] - b_param["u"][0])
/ b_param["u"][1])) + b_param["u"][3]) ** b_param["u"][2]
u_tau_Q10 = 2.7 ** ((temp - 25) / 10)
tau_u = (86.86 + (408.78 / (1 + np.exp((V_m[:, t - 1] - shift["u"] + 13.6) / 7.46)))) / u_tau_Q10
u_dot = (u_inf - gating[ind_dict['u'], :]) * (1 / (tau_u * scale["u"]))
gating[ind_dict['u'], :] = gating[ind_dict['u'], :] + u_dot * dt
gating[ind_dict['u'], :][gating[ind_dict['u'], :] > 1.] = 1.
gating[ind_dict['u'], :][gating[ind_dict['u'], :] < 0.] = 0.
# mutant Kv1.1 activation
s_mut_inf = (1 / (1 + np.exp((V_m[:, t - 1] - shift["s_mut"] - slope_shift["s_mut"] - b_param["s_mut"][0]) / b_param["s_mut"][1]))) ** \
b_param["s_mut"][2]
s_mut_dot = (s_mut_inf - gating[ind_dict['s_mut'], :]) * (1 / (tau_s * scale["s_mut"]))
gating[ind_dict['s_mut'], :] = gating[ind_dict['s_mut'], :] + s_mut_dot * dt
gating[ind_dict['s_mut'], :][gating[ind_dict['s_mut'], :] > 1.] = 1.
gating[ind_dict['s_mut'], :][gating[ind_dict['s_mut'], :] < 0.] = 0.
# mutant Kv1.1 inactivation
u_mut_inf = ((1 - b_param["u_mut"][3]) / (1 + np.exp((V_m[:, t - 1] - shift["u_mut"] - slope_shift["u_mut"] - b_param["u_mut"][0])
/ b_param["u_mut"][1])) + b_param["u_mut"][3]) ** b_param["u_mut"][2]
u_mut_dot = (u_mut_inf - gating[ind_dict['u_mut'], :]) * (1 / (tau_u * scale["u_mut"]))
gating[ind_dict['u_mut'], :] = gating[ind_dict['u_mut'], :] + u_mut_dot * dt
gating[ind_dict['u_mut'], :][gating[ind_dict['u_mut'], :] > 1.] = 1.
gating[ind_dict['u_mut'], :][gating[ind_dict['u_mut'], :] < 0.] = 0.
# 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['Kv'], :] = g["Kv"] * gating[ind_dict['s'], :] * gating[ind_dict['u'], :] * (
V_m[:, t - 1] - E["K"]) * currents_included["Kv"]
current[ind_dict['Kv_mut'], :] = g["Kv_mut"] * gating[ind_dict['s_mut'], :] * gating[ind_dict['u_mut'], :] * \
(V_m[:, t - 1] - E["K"]) * currents_included["Kv_mut"]
current[ind_dict['A'], :] = g["A"] * gating[ind_dict['n_A'], :] * gating[ind_dict['h_A'], :] * (
V_m[:, t - 1] - E["K"]) * currents_included["A"]
current[ind_dict['T'], :] = g["T"] * gating[ind_dict['m_T'], :] * gating[ind_dict['h_T'], :] * (
V_m[:, t - 1] - E["Ca"]) * currents_included["T"]
# 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['A'], :] +
current[ind_dict['Kv'], :] + current[ind_dict['Kv_mut'], :] +
current[ind_dict['T'], :]))
# update V_m
V_m[:, t] = V_m[:, t - 1] + V_dot * dt
t += 1
return V_m
def Cb_stellate_Kv_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):
""" Simulate mutations for Cb Stellate Kv 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
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', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'Kv', 'Kv_mut', 'T', 'Leak']):
ind_dict[var] = i
i += 1
# effect of mutation
g_effect['Kv_mut'] = g['Kv_mut'] * mutations[mut]['g_ratio'] # mutation_properties[mut, 2] # g ratio A
b_param_effect['s_mut'][0] = b_param_effect['s_mut'][0] + mutations[mut]['activation_Vhalf_diff']
b_param_effect['s_mut'][1] = b_param_effect['s_mut'][1] * mutations[mut]['activation_k_ratio']
# initial simulation for current range from low to high
V_m = Cb_stellate_Kv(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)
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, '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', 'spike_times_rheo', '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('AP_rheo', (I_mag.shape[0],), dtype='bool', compression="gzip", compression_opts=9)
analysis.create_dataset('rheobase', (1,), dtype='float64', compression="gzip", compression_opts=9)
# Firing analysis
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 = Cb_stellate_Kv(V_init * np.ones((stim_num_rheo)), g_effect, E2, I_in_rheo, f['data'].attrs['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'].create_dataset('V_m_rheo', data=V_m_rheo, dtype='float64', compression="gzip", compression_opts=9)
# run AP detection
stim_start = np.int(initial_period * 1 / dt)
for stim in range(I_in_rheo.shape[1]):
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:
# 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 = Cb_stellate_Kv(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)
# 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 = Cb_stellate_Kv(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)
# 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']):
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_rel_acc', shape=(1,), dtype='float64', compression="gzip",
compression_opts=9)
f['analysis']['AUC_rel_acc'][:] = 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
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 = Cb_stellate_Kv(-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)
# 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_Cb_stellate_Kv(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, lin_array, log_array, alt_types, alt_ind, alt):
""" Sensitivity analysis for Cb Stellate Kv 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
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', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'Kv', 'Kv_mut', 'T', 'Leak']):
ind_dict[var] = i
i += 1
# gating 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 = Cb_stellate_Kv(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)
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, '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,}
data.attrs.update(metadata)
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_names:
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)
analysis.create_dataset('rheobase', (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 = Cb_stellate_Kv(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)
# run AP detection
try:
for stim in range(I_in_rheo.shape[1]):
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
# find rheobase of steady state firing (F_inf)
with h5py.File(fname, "a") as f:
try:
# setup current range to be in current step between no and start of steady state firing
F_inf_start_ind = np.argwhere(f['analysis']['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
# simulate response to this current step range
V_m_inf_rheo = Cb_stellate_Kv(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)
# 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 * f['data'].attrs['I_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 = Cb_stellate_Kv(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)
# 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()

File diff suppressed because it is too large Load Diff

1005
Code/Functions/STN_fxns.py Normal file

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,135 @@
import json
import numpy as np
from numba import njit
@njit
def sqrt_log(x, *p):
""" Function of sum of square root and logarithmic function to fit to fI curves
Parameters
----------
x : array
input x value
*p : array
parameters for function
Returns
-------
p[0] * np.sqrt(x + p[1]) + p[2] * np.log(x + (p[1] + 1)) + p[3]
"""
return p[0] * np.sqrt(x + p[1]) + p[2] * np.log(x + (p[1] + 1)) + p[3]
@njit
def capacitance(diam, C_m):
""" calculate the capacitance of a cell's membrane give the diameter
Parameters
----------
diam: float
diameter of spherical cell body
C_m: float
specific membrane capacitance of cell
Returns
-------
capacitance: float
capacitance of cell
surf_area: flaot
surface area of cell
"""
r = diam * 1E-4 / 2 # cm
l = r
surf_area = 2 * np.pi * r ** 2 + 2 * np.pi * r * l # cm^2
C = C_m * surf_area # uF
return (C, surf_area)
@njit
def stimulus_init(low, high, number_steps, initial_period, dt, sec):
""" initiation of stimulus pulse from stimulus magnitudes I_in_mag
Parameters
----------
low : float
lowest current step magnitude
high : float
upper current step magnitude
number_steps : int
number of current steps between low and high
initial_period : float
length of I=0 before current step
dt : float
time step
sec:
length of current step in seconds
Returns
-------
stim_time: float
length of I_in in samples
I_in: array
I input array
stim_num: int
number of current steps
V_m: array
membrane potential array initialization of shape:(stim_num, stim_time)
"""
stim_time = np.int(1000 * 1 / dt * sec) # 1000 msec/sec * 1/dt
I_in_mag = np.arange(low, high, (high - low) / number_steps)
I_in_mag = np.reshape(I_in_mag, (1, I_in_mag.shape[0]))
I_in = np.zeros((stim_time + np.int(initial_period * 1 / dt), 1)) @ I_in_mag
I_in[np.int(initial_period * 1 / dt):, :] = np.ones((stim_time, 1)) @ I_in_mag
stim_num = I_in.shape[1] # number of different currents injected
stim_time = stim_time + np.int(initial_period * 1 / dt)
V_m = np.zeros((stim_num, stim_time))
return stim_time, I_in, stim_num, V_m
class NumpyEncoder(json.JSONEncoder):
""" Special json encoder for numpy types """
def default(self, obj):
if isinstance(obj, np.integer):
return int(obj)
elif isinstance(obj, np.floating):
return float(obj)
elif isinstance(obj, np.ndarray):
return obj.tolist()
return json.JSONEncoder.default(self, obj)
def init_dict(variable):
""" Initialize dictionaries for simulations
Parameters
----------
variable : array
array of variable name strings to put in dicts
Returns
-------
shift : dict
with every variable name in variable = 0.
scale : dict
with every variable name in variable = 1.
slope_shift : dict
with every variable name in variable = 0.
E : dict
empty dictionary
currents_included : dict
empty dictionary
b_param : dict
with every variable name in variable = np.array of 2 zeros
g : dict
empty dictionary
"""
shift = {}
scale = {}
slope_shift = {}
b_param ={}
g = {}
E ={}
currents_included = {}
for var in variable:
shift[var] = 0.
scale[var] = 1.0
slope_shift[var] = 0.
b_param[var] = np.zeros(3, dtype=np.dtype('float64'))
return shift, scale, slope_shift, E, currents_included, b_param, g

View File

@ -0,0 +1,112 @@
# -*- coding: utf-8 -*-
"""
Created on Sat May 8 16:20:11 2021
@author: nils
"""
import numpy as np
from joblib import Parallel, delayed
import json
import os
from numba import types
from numba.typed import Dict
from Code.Functions.Utility_fxns import capacitance, stimulus_init, init_dict
from Code.Functions.Cb_stellate_fxns import Cb_stellate_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.0015
number_steps = 200
initial_period = 1000
num_gating = 9
num_current = 6
C, surf_area = capacitance(61.4, 1.50148)
variable = np.array(['m', 'h', 'n', 'n_A', 'h_A', 'n_A_mut', 'h_A_mut', 'm_T', 'h_T'])
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'n_A', 'h_A','n_A_mut', 'h_A_mut', 'm_T', 'h_T']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-37., -3, 1])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-40., 4., 1., 0])
b_param['n'][:] = np.array([-23, -5, 1])
b_param['n_A'][:] = np.array([-27, -13.2, 1.])
b_param['h_A'][:] = np.array([-80., 6.5, 1.])
b_param['n_A_mut'][:] = np.array([-27, -13.2, 1.])
b_param['h_A_mut'][:] = np.array([-80., 6.5, 1.])
b_param['m_T'][:] = np.array([-50., -3, 1.])
b_param['h_T'][:] = np.array([-68., 3.75, 1.])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 55.
E["K"] = -80.
E["Ca"] = 22.
E["Leak"] = -70. # as per Molineux et al 2005 and NOT the -38 in Alexander et al 2019
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["A_mut"] = True
currents_included["A"] = True
currents_included["T"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 3.4 * surf_area
g["Kd"] = 9.0556 * surf_area
g["A_mut"] = 15.0159 * surf_area * 0.5
g["A"] = 15.0159 * surf_area* 0.5
g["T"] = 0.45045 * surf_area
g["Leak"] = 0.07407 * surf_area
# save folder
folder = './KCNA1_mutations/Cb_stellate'
if not os.path.isdir(folder):
os.makedirs(folder)
#
# mutation properties
mutations = json.load(open("mutations_effects_dict.json"))
# %%
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(Cb_stellate_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) for mut in list(mutations.keys()))

View File

@ -0,0 +1,121 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 18 10:47:40 2021
@author: nils
"""
import numpy as np
from joblib import Parallel, delayed
import json
import os
from numba import types
from numba.typed import Dict
from Code.Functions.Utility_fxns import capacitance, stimulus_init, init_dict
from Code.Functions.Cb_stellate_fxns_Kv import Cb_stellate_Kv_mut
#%%
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.0015
number_steps = 200
initial_period = 1000
num_gating = 11
num_current = 7
C, surf_area = capacitance(61.4, 1.50148)
variable = np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut'])
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-37., -3, 1])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-40., 4., 1., 0])
b_param['n'][:] = np.array([-23, -5, 1])
b_param['n_A'][:] = np.array([-27, -13.2, 1.])
b_param['h_A'][:] = np.array([-80., 6.5, 1.])
b_param['m_T'][:] = np.array([-50., -3, 1.])
b_param['h_T'][:] = np.array([-68., 3.75, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 55.
E["K"] = -80.
E["Ca"] = 22.
E["Leak"] = -70. # as per Molineux et al 2005 and NOT the -38 in Alexander et al 2019
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["A"] = False
currents_included["T"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 34 * surf_area
g["Kd"] = 9.0556 * surf_area
if currents_included["Kv_mut"] == True:
g["Kv"] = 1.50159 / 2 * surf_area
else:
g["Kv"] = 1.50159 / 2 * surf_area * 2
g["Kv_mut"] = 1.50159 / 2 * surf_area
g["A"] = 0 * surf_area
g["T"] = 0.45045 * surf_area
g["Leak"] = 0.07407 * surf_area
# save folder
folder = './KCNA1_mutations/Cb_stellate_Kv_only'
if not os.path.isdir(folder):
os.makedirs(folder)
#
# mutation properties
mutations = json.load(open("mutations_effects_dict.json"))
# %%
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(Cb_stellate_Kv_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) for mut in list(mutations.keys()))

View File

@ -0,0 +1,119 @@
# -*- coding: utf-8 -*-
"""
Created on Sat May 8 16:20:11 2021
@author: nils
"""
import numpy as np
from joblib import Parallel, delayed
import json
import os
from numba import types
from numba.typed import Dict
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Cb_stellate_fxns_Kv import Cb_stellate_Kv_mut
#%%
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.0015
number_steps = 200
initial_period = 1000
num_gating = 11
num_current = 7
C, surf_area = capacitance(61.4, 1.50148)
variable = np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut'])
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-37., -3, 1])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-40., 4., 1., 0])
b_param['n'][:] = np.array([-23, -5, 1])
b_param['n_A'][:] = np.array([-27, -13.2, 1.])
b_param['h_A'][:] = np.array([-80., 6.5, 1.])
b_param['m_T'][:] = np.array([-50., -3, 1.])
b_param['h_T'][:] = np.array([-68., 3.75, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 55.
E["K"] = -80.
E["Ca"] = 22.
E["Leak"] = -70. # as per Molineux et al 2005 and NOT the -38 in Alexander et al 2019
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["A"] = True
currents_included["T"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 3.4 * surf_area
g["Kd"] = 9.0556 * (1-Kv_ratio) * surf_area
g["Kv"] = 6. * Kv_ratio/2 * surf_area * (2 * int(currents_included["Kv_mut"] == False))
g["Kv_mut"] = 6. * Kv_ratio/2 * surf_area
g["A"] = 15.0159 * surf_area
g["T"] = 0.45045 * surf_area
g["Leak"] = 0.07407 * surf_area
# save folder
folder = './KCNA1_mutations/Cb_stellate_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
#
# mutation properties
mutations = json.load(open("mutations_effects_dict.json"))
# %%
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(Cb_stellate_Kv_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) for mut in list(mutations.keys()))

View File

@ -0,0 +1,119 @@
# -*- coding: utf-8 -*-
"""
Created on Sat May 29 21:10:20 2021
@author: nils
"""
import numpy as np
import os
import json
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
from Code.Functions.Utility import capacitance, stimulus_init, init_dict, NumpyEncoder
from Code.Functions.Pospischil_fxns import Pospischil_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(56.9, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 502
V_init = -70
V_T = -57.9
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param = {}
b_param['m'] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.array([-34.51951036, 4.04059373, 1., 0.])
b_param['n'] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'] = np.array([-57.37, 20.98, 1.])
b_param['p'] = np.array([-45., -9.9998807337, 1.])
b_param['s'] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.array([-31., 5.256, 1., 0.245])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.4
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 58. * surf_area
g["Kd"] = 3.9 * (1 - Kv_ratio) * surf_area
g["M"] = 0.075 * surf_area
if currents_included["Kv_mut"] == True:
g["Kv"] = 3.9 * Kv_ratio / 2 * surf_area
else:
g["Kv"] = 3.9 * Kv_ratio / 2 * surf_area * 2
g["Kv_mut"] = 3.9 * Kv_ratio / 2 * surf_area
g["L"] = 0. * surf_area
g["Leak"] = 0.038 * surf_area
# save folder
folder = './KCNA1_mutations/FS_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
# mutation properties
mutations = json.load(open("mutations_effects_dict.json"))
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(Pospischil_mut)(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift, scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, mutations, mut, folder,
high, low, number_steps, initial_period, sec)
for mut in list(mutations.keys()))

View File

@ -0,0 +1,126 @@
# -*- coding: utf-8 -*-
"""
Created on Sat May 29 21:10:48 2021
@author: nils
"""
import numpy as np
from joblib import Parallel, delayed
import os
import h5py
import json
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Pospischil_fxns import Pospischil_mut
# model parameters
tau_max_p = 608
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.8, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 934
V_init = -70
V_T = -67.9
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param = {}
b_param['m'] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'] = np.array([-57.37, 20.98, 1.0])
b_param['p'] = np.array([-45., -9.9998807337, 1.])
b_param['s'] = np.array([-14.16, -10.15, 1.0])
b_param['u'] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.array([-31., 5.256, 1., 0.245])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -56.2
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 10. * surf_area
g["Kd"] = 2.1 * (1 - Kv_ratio) * surf_area
g["M"] = 0.0098 * surf_area
if currents_included["Kv_mut"] == True:
g["Kv"] = 2.1 * Kv_ratio / 2 * surf_area
else:
g["Kv"] = 2.1 * Kv_ratio / 2 * surf_area * 2
g["Kv_mut"] = 2.1 * Kv_ratio / 2 * surf_area
g["L"] = 0. * surf_area
g["Leak"] = 0.0205 * surf_area
# save folder
folder = './KCNA1_mutations/RS_inhib_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
# mutation properties
mutations = json.load(open("mutations_effects_dict.json"))
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(Pospischil_mut)(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift,
scale, b_param, slope_shift, gating, current, prominence, desired_AUC_width, mutations, mut,
folder, high, low, number_steps, initial_period, sec)
for mut in list(mutations.keys()))

View File

@ -0,0 +1,120 @@
# -*- coding: utf-8 -*-
"""
Created on Sat May 8 16:20:11 2021
@author: nils
"""
import numpy as np
import os
import json
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Pospischil_fxns import Pospischil_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.4, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 608
V_init = -70
V_T = -56.2
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param = {}
b_param['m'] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'] = np.array([-57.37, 20.98, 1.0])
b_param['p'] = np.array([-45., -9.9998807337, 1.])
b_param['s'] = np.array([-14.16, -10.15, 1.0])
b_param['u'] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.array([-31., 5.256, 1., 0.245])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.3
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 56. * surf_area
g["Kd"] = 6. * (1 - Kv_ratio) * surf_area
g["M"] = 0.075 * surf_area
if currents_included["Kv_mut"] == True:
g["Kv"] = 6. * Kv_ratio / 2 * surf_area
else:
g["Kv"] = 6. * Kv_ratio / 2 * surf_area * 2
g["Kv_mut"] = 6. * Kv_ratio / 2 * surf_area
g["L"] = 0. * surf_area
g["Leak"] = 0.0205 * surf_area
# save folder
folder = './KCNA1_mutations/RS_pyramidal_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
# mutation properties
mutations = json.load(open("mutations_effects_dict.json"))
# # %%
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(Pospischil_mut)(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift,
scale, b_param, slope_shift, gating, current, prominence, desired_AUC_width, mutations, mut,
folder, high, low, number_steps, initial_period, sec)
for mut in list(mutations.keys()))

View File

@ -0,0 +1,107 @@
# -*- coding: utf-8 -*-
"""
Created on Sat May 8 16:20:11 2021
@author: nils
"""
import numpy as np
from joblib import Parallel, delayed
import json
import os
from numba import types
from numba.typed import Dict
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.STN_fxns import STN_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.005
number_steps = 100
initial_period = 1000
num_gating = 15
num_current = 8
d = 61.4
r = d/2 * 10**-6 # radius in meters
vol = np.pi * r**3 # volume in liters
C, surf_area = capacitance(d, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'a', 'b', 'a_mut', 'b_mut', 'c', 'd1', 'd2', 'p', 'q', 'r', 'Ca_conc', 'E_Ca']))
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'] = np.array([-40., -1. , 1.])
b_param['h'] = np.zeros(4)
b_param['h'] = np.array([-45.5, 1. , 1., 0.05])
b_param['n'] = np.array([-41., -1., 1.])
b_param['p'] = np.array([-56., -1., 1.])
b_param['q'] = np.array([-85., 1., 1.])
b_param['r'] = np.array([0.17, -0.08, 1.])
b_param['a'] = np.array([-45., -14.7, 1.])
b_param['b'] = np.array([-90., 7.5, 1.])
b_param['a_mut'] = np.array([-45., -14.7, 1.])
b_param['b_mut'] = np.array([-90., 7.5, 1.])
b_param['c'] = np.array([-30.6, -5., 1.])
b_param['d1'] = np.array([-60, 7.5, 1.])
b_param['d2'] = np.array([0.1, 0.02, 1.])
# reversal potentials
E["Na"] = 60.
E["K"] = -90.
E["Leak"] = -60.
# model currents
currents_included["Na"] = True
currents_included["Kd"] =True
currents_included["L"] = True
currents_included["T"] = True
currents_included["Ca_K"] = True
currents_included["A"] =True
currents_included["A_mut"] =True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 49. * surf_area
g["Kd"] = 27. * (1-Kv_ratio)* surf_area
g["A"] = 5. * surf_area *0.5
g["A_mut"] = 5. * surf_area *0.5
g["L"] = 5. * surf_area
g["T"] = 5. * surf_area
g["Ca_K"] = 1 * surf_area
g["Leak"] = 0.035 * surf_area
V_init = -70
# save folder
folder = './KCNA1_mutations/STN'
if not os.path.isdir(folder):
os.makedirs(folder)
#
# mutation properties
mutations = json.load(open("mutations_effects_dict.json"))
prominence = 50
desired_AUC_width = high/5
#%%
Parallel(n_jobs=8, verbose=9)(
delayed(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) for mut in list(mutations.keys()))

View File

@ -0,0 +1,129 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 18 10:53:45 2021
@author: nils
"""
import numpy as np
from joblib import Parallel, delayed
import json
import os
from numba import types
from numba.typed import Dict
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.STN_fxns_Kv import STN_Kv_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.005
number_steps = 100
initial_period = 1000
num_gating = 17
num_current = 9
d = 61.4
r = d/2 * 10**-6 # radius in meters
vol = np.pi * r**3 # volume in liters
C, surf_area = capacitance(d, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'a', 'b', 'c', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut','Ca_conc', 'E_Ca']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'] = np.array([-40., -1. , 1.])
b_param['h'] = np.zeros(4)
b_param['h'] = np.array([-45.5, 1. , 1., 0.05])
b_param['n'] = np.array([-41., -1., 1.])
b_param['p'] = np.array([-56., -1., 1.])
b_param['q'] = np.array([-85., 1., 1.])
b_param['r'] = np.array([0.17, -0.08, 1.])
b_param['a'] = np.array([-45., -14.7, 1.])
b_param['b'] = np.array([-90., 7.5, 1.])
b_param['c'] = np.array([-30.6, -5., 1.])
b_param['d1'] = np.array([-60, 7.5, 1.])
b_param['d2'] = np.array([0.1, 0.02, 1.])
b_param['s'] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'] = np.array([-31., 5.256, 1., 0.245])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 60.
E["K"] = -90.
E["Leak"] = -60.
# model currents
currents_included["Na"] = True
currents_included["Kd"] =True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["L"] = True
currents_included["T"] = True
currents_included["Ca_K"] = True
currents_included["A"] =True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 49. * surf_area
g["Kd"] = 27. * (1-Kv_ratio)* surf_area
if currents_included["Kv_mut"]==True:
g["Kv"] = 27. * Kv_ratio/2 * surf_area
else:
g["Kv"] = 27. * Kv_ratio/2 * surf_area * 2
g["Kv_mut"] = 27. * Kv_ratio/2 * surf_area
g["A"] = 5. * surf_area
g["L"] = 5. * surf_area
g["T"] = 5. * surf_area
g["Ca_K"] = 1 * surf_area
g["Leak"] = 0.035 * surf_area
# save folder
folder = './KCNA1_mutations/STN_Kv_only'
if not os.path.isdir(folder):
os.makedirs(folder)
#
# mutation properties
mutations = json.load(open("mutations_effects_dict.json"))
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(STN_Kv_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) for mut in list(mutations.keys()))

View File

@ -0,0 +1,129 @@
# -*- coding: utf-8 -*-
"""
Created on Sat May 8 16:20:11 2021
@author: nils
"""
import numpy as np
from joblib import Parallel, delayed
import json
import os
from numba import types
from numba.typed import Dict
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.STN_fxns_Kv import STN_Kv_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.005
number_steps = 100
initial_period = 1000
num_gating = 17
num_current = 9
d = 61.4
r = d/2 * 10**-6 # radius in meters
vol = np.pi * r**3 # volume in liters
C, surf_area = capacitance(d, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'a', 'b', 'c', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut','Ca_conc', 'E_Ca']))
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'] = np.array([-40., -1. , 1.])
b_param['h'] = np.zeros(4)
b_param['h'] = np.array([-45.5, 1. , 1., 0.05])
b_param['n'] = np.array([-41., -1., 1.])
b_param['p'] = np.array([-56., -1., 1.])
b_param['q'] = np.array([-85., 1., 1.])
b_param['r'] = np.array([0.17, -0.08, 1.])
b_param['a'] = np.array([-45., -14.7, 1.])
b_param['b'] = np.array([-90., 7.5, 1.])
b_param['c'] = np.array([-30.6, -5., 1.])
b_param['d1'] = np.array([-60, 7.5, 1.])
b_param['d2'] = np.array([0.1, 0.02, 1.])
b_param['s'] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'] = np.array([-31., 5.256, 1., 0.245])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 60.
E["K"] = -90.
E["Leak"] = -60.
# model currents
currents_included["Na"] = True
currents_included["Kd"] =True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["L"] = True
currents_included["T"] = True
currents_included["Ca_K"] = True
currents_included["A"] =True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 49. * surf_area
g["Kd"] = 27. * (1-Kv_ratio)* surf_area
if currents_included["Kv_mut"]==True:
g["Kv"] = 27. * Kv_ratio/2 * surf_area
else:
g["Kv"] = 27. * Kv_ratio/2 * surf_area * 2
g["Kv_mut"] = 27. * Kv_ratio/2 * surf_area
g["A"] = 5. * surf_area
g["L"] = 5. * surf_area
g["T"] = 5. * surf_area
g["Ca_K"] = 1 * surf_area
g["Leak"] = 0.035 * surf_area
V_init = -70
# save folder
folder = './KCNA1_mutations/STN_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
#
# mutation properties
mutations = json.load(open("mutations_effects_dict.json"))
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(STN_Kv_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) for mut in list(mutations.keys()))

View File

@ -0,0 +1,88 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 18 10:47:40 2021
@author: nils
"""
import numpy as np
import os
from Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Cb_stellate_fxns_Kv import Cb_stellate_Kv_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.0015
number_steps = 200
initial_period = 1000
num_gating = 11
num_current = 7
C, surf_area = capacitance(61.4, 1.50148)
variable = np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut'])
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# gating parameters
b_param['m'][:] = np.array([-37., -3, 1])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-40., 4., 1., 0])
b_param['n'][:] = np.array([-23, -5, 1])
b_param['n_A'][:] = np.array([-27, -13.2, 1.])
b_param['h_A'][:] = np.array([-80., 6.5, 1.])
b_param['m_T'][:] = np.array([-50., -3, 1.])
b_param['h_T'][:] = np.array([-68., 3.75, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 55.
E["K"] = -80.
E["Ca"] = 22.
E["Leak"] = -70. # as per Molineux et al 2005 and NOT the -38 in Alexander et al 2019
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["A"] = False
currents_included["T"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 34 * surf_area
g["Kd"] = 9.0556 * surf_area
if currents_included["Kv_mut"] == True:
g["Kv"] = 1.50159 / 2 * surf_area
else:
g["Kv"] = 1.50159 / 2 * surf_area * 2
g["Kv_mut"] = 1.50159 / 2 * surf_area
g["A"] = 0 * surf_area
g["T"] = 0.45045 * surf_area
g["Leak"] = 0.07407 * surf_area
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
mut = 'Cb_stellate_Delta_Kv'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
if not os.path.isdir(folder):
os.makedirs(folder)
Cb_stellate_Kv_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)

View File

@ -0,0 +1,87 @@
# -*- coding: utf-8 -*-
"""
Created on Sat May 8 16:20:11 2021
@author: nils
"""
import numpy as np
import os
from Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Cb_stellate_fxns_Kv import Cb_stellate_Kv_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.0015
number_steps = 200
initial_period = 1000
num_gating = 11
num_current = 7
C, surf_area = capacitance(61.4, 1.50148)
variable = np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut'])
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# gating parameters
b_param['m'][:] = np.array([-37., -3, 1])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-40., 4., 1., 0])
b_param['n'][:] = np.array([-23, -5, 1])
b_param['n_A'][:] = np.array([-27, -13.2, 1.])
b_param['h_A'][:] = np.array([-80., 6.5, 1.])
b_param['m_T'][:] = np.array([-50., -3, 1.])
b_param['h_T'][:] = np.array([-68., 3.75, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 55.
E["K"] = -80.
E["Ca"] = 22.
E["Leak"] = -70. # as per Molineux et al 2005 and NOT the -38 in Alexander et al 2019
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["A"] = True
currents_included["T"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 3.4 * surf_area
g["Kd"] = 9.0556 * (1-Kv_ratio) * surf_area
g["Kv"] = 6. * Kv_ratio/2 * surf_area * (2 * int(currents_included["Kv_mut"] == False))
g["Kv_mut"] = 6. * Kv_ratio/2 * surf_area
g["A"] = 15.0159 * surf_area
g["T"] = 0.45045 * surf_area
g["Leak"] = 0.07407 * surf_area
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
mut = 'Cb_stellate_Kv'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
if not os.path.isdir(folder):
os.makedirs(folder)
Cb_stellate_Kv_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)

View File

@ -0,0 +1,73 @@
import numpy as np
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Cb_stellate_fxns import Cb_stellate_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.0015
number_steps = 200
initial_period = 1000
num_gating = 9
num_current = 6
C, surf_area = capacitance(61.4, 1.50148)
variable = np.array(['m', 'h', 'n', 'n_A', 'h_A', 'n_A_mut', 'h_A_mut', 'm_T', 'h_T'])
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'n_A', 'h_A','n_A_mut', 'h_A_mut', 'm_T', 'h_T']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# gating parameters
b_param['m'][:] = np.array([-37., -3, 1])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-40., 4., 1., 0])
b_param['n'][:] = np.array([-23, -5, 1])
b_param['n_A'][:] = np.array([-27, -13.2, 1.])
b_param['h_A'][:] = np.array([-80., 6.5, 1.])
b_param['n_A_mut'][:] = np.array([-27, -13.2, 1.])
b_param['h_A_mut'][:] = np.array([-80., 6.5, 1.])
b_param['m_T'][:] = np.array([-50., -3, 1.])
b_param['h_T'][:] = np.array([-68., 3.75, 1.])
# reversal potentials
E["Na"] = 55.
E["K"] = -80.
E["Ca"] = 22.
E["Leak"] = -70. # as per Molineux et al 2005 and NOT the -38 in Alexander et al 2019
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["A_mut"] = True
currents_included["A"] = True
currents_included["T"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 3.4 * surf_area
g["Kd"] = 9.0556 * surf_area
g["A_mut"] = 15.0159 * surf_area * 0.5
g["A"] = 15.0159 * surf_area* 0.5
g["T"] = 0.45045 * surf_area
g["Leak"] = 0.07407 * surf_area
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
mut = 'Cb_stellate'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
if not os.path.isdir(folder):
os.makedirs(folder)
Cb_stellate_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)

View File

@ -0,0 +1,100 @@
# -*- coding: utf-8 -*-
"""
Created on Fri May 7 23:37:22 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Pospischil_fxns import Pospischil_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(56.9, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array([-45., -9.9998807337, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.4
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
Kv_ratio = 0.10
g["Na"] = 58. * surf_area
g["Kd"] = 3.9 * (1 - Kv_ratio) * surf_area
g["M"] = 0.075 * surf_area
if currents_included["Kv_mut"] == True:
g["Kv"] = 3.9 * Kv_ratio / 2 * surf_area
else:
g["Kv"] = 3.9 * Kv_ratio / 2 * surf_area * 2
g["Kv_mut"] = 3.9 * Kv_ratio / 2 * surf_area
g["L"] = 0. * surf_area
g["Leak"] = 0.038 * surf_area
tau_max_p = 502
V_init = -70
V_T = -57.9
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
if not os.path.isdir(folder):
os.makedirs(folder)
mut = 'FS_Kv'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
Pospischil_mut(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift, scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, mutations, mut, folder,
high, low, number_steps, initial_period, sec)

View File

@ -0,0 +1,101 @@
# -*- coding: utf-8 -*-
"""
Created on Fri May 7 23:37:22 2021
@author: nils
Original Pospischil FS (no Kv1.1)
"""
# Todo: error when run
import numpy as np
from numba import types
from numba.typed import Dict
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Pospischil_fxns import Pospischil_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(56.9, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array([-45., -9.9998807337, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.4
# currents in model
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = False
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 58. * surf_area
g["Kd"] = 3.9 * surf_area
g["M"] = 0.075 * surf_area
g["Kv"] = 0.
g["Kv_mut"] = 0.
g["L"] = 0.
g["Leak"] = 0.038 * surf_area
tau_max_p = 502
V_init = -70
V_T = -57.9
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
if not os.path.isdir(folder):
os.makedirs(folder)
mut = 'FS'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
Pospischil_mut(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift, scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, mutations, mut, folder,
high, low, number_steps, initial_period, sec)

View File

@ -0,0 +1,114 @@
# -*- coding: utf-8 -*-
"""
Created on Fri May 7 23:37:22 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Pospischil_fxns import Pospischil_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.00035
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.8, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 934
V_init = -70
V_T = -67.9
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array([-45., -9.9998807337, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -56.2
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 10. * surf_area
g["Kd"] = 2.1 * (1 - Kv_ratio) * surf_area
g["M"] = 0.0098 * surf_area
if currents_included["Kv_mut"] == True:
g["Kv"] = 2.1 * Kv_ratio / 2 * surf_area
else:
g["Kv"] = 2.1 * Kv_ratio / 2 * surf_area * 2
g["Kv_mut"] = 2.1 * Kv_ratio / 2 * surf_area
g["L"] = 0. * surf_area
g["Leak"] = 0.0205 * surf_area
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5 #0.0005
folder = '../Neuron_models'
if not os.path.isdir(folder):
os.makedirs(folder)
mut = 'RS_inhib_Kv'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
Pospischil_mut(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift, scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, mutations, mut, folder,
high, low, number_steps, initial_period, sec)

View File

@ -0,0 +1,102 @@
# -*- coding: utf-8 -*-
"""
Created on Fri May 7 23:37:22 2021
@author: nils
Original Pospischil RS inhibitory (no Kv1.1)
"""
import numpy as np
from numba import types
from numba.typed import Dict
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Pospischil_fxns import Pospischil_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.00035
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.8, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 934
V_init = -70
V_T = -67.9
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array([-45., -9.9998807337, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potential
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -56.2
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = False
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
g["Na"] = 10. * surf_area
g["Kd"] = 2.1 * surf_area
g["M"] = 0.0098 * surf_area
g["Kv"] = 0.
g["Kv_mut"] = 0.
g["L"] = 0.
g["Leak"] = 0.0205 * surf_area
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
if not os.path.isdir(folder):
os.makedirs(folder)
mut = 'RS_inhib'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
Pospischil_mut(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift, scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, mutations, mut, folder,
high, low, number_steps, initial_period, sec)

View File

@ -0,0 +1,117 @@
# -*- coding: utf-8 -*-
"""
Created on Fri May 7 23:37:22 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Pospischil_fxns import Pospischil_mut
# model parameters
tau_max_p = 608
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.4, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
V_init = -70
V_T = -56.2
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array([-45., -9.9998807337, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
mut_act_Vhalf_wt = -30.01851851851851
mut_act_k_wt = -7.7333333333333325
s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
s_diff_k = mut_act_k_wt - b_param['s'][1]
b_param['s'][1] = b_param['s'][1] + s_diff_k
b_param['u'][1] = b_param['u'][1] + s_diff_k
b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.3
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = False
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 56. * surf_area
g["Kd"] = 6. * (1 - Kv_ratio) * surf_area
g["M"] = 0.075 * surf_area
if currents_included["Kv_mut"] == True:
g["Kv"] = 6. * Kv_ratio / 2 * surf_area
else:
g["Kv"] = 6. * Kv_ratio / 2 * surf_area * 2
g["Kv_mut"] = 6. * Kv_ratio / 2 * surf_area
g["L"] = 0. * surf_area
g["Leak"] = 0.0205 * surf_area
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
if not os.path.isdir(folder):
os.makedirs(folder)
mut = 'RS_pyramidal_Kv'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
Pospischil_mut(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift, scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, mutations, mut, folder,
high, low, number_steps, initial_period, sec)

View File

@ -0,0 +1,99 @@
# -*- coding: utf-8 -*-
"""
Created on Fri May 7 23:37:22 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Pospischil_fxns import Pospischil_mut
# model parameters
tau_max_p = 608
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.4, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
V_init = -70
V_T = -56.2
# initialize parameters
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array(
[-45., -9.9998807337, 1.]) # -45 with 10 mV shift so that it contributes to resting potential
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.3
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = False
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 56. * surf_area
g["Kd"] = 6. * surf_area
g["M"] = 0.075 * surf_area
g["Kv"] = 0.
g["Kv_mut"] = 0.
g["L"] = 0.
g["Leak"] = 0.0205 * surf_area
prominence = 50
min_spike_height = 0
desired_AUC_width = high / 5
folder = '../Neuron_models'
if not os.path.isdir(folder):
os.makedirs(folder)
mut = 'RS_pyramidal'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
Pospischil_mut(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift, scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, mutations, mut, folder,
high, low, number_steps, initial_period, sec)

View File

@ -0,0 +1,111 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 18 10:53:45 2021
@author: nils
"""
import numpy as np
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.STN_fxns_Kv import STN_Kv_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.005
number_steps = 100
initial_period = 1000
num_gating = 17
num_current = 9
d = 61.4
r = d/2 * 10**-6 # radius in meters
vol = np.pi * r**3 # volume in liters
C, surf_area = capacitance(d, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'a', 'b', 'c', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut','Ca_conc', 'E_Ca']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = {}
i = 0
for var in np.array(['m', 'h', 'n', 'a', 'b', 'c', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut', 'Ca_conc','E_Ca']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'L', 'Kv', 'Kv_mut', 'T', 'Leak', 'Ca_K']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'] = np.array([-40., -1. , 1.])
b_param['h'] = np.zeros(4)
b_param['h'] = np.array([-45.5, 1. , 1., 0.05])
b_param['n'] = np.array([-41., -1., 1.])
b_param['p'] = np.array([-56., -1., 1.])
b_param['q'] = np.array([-85., 1., 1.])
b_param['r'] = np.array([0.17, -0.08, 1.])
b_param['a'] = np.array([-45., -14.7, 1.])
b_param['b'] = np.array([-90., 7.5, 1.])
b_param['c'] = np.array([-30.6, -5., 1.])
b_param['d1'] = np.array([-60, 7.5, 1.])
b_param['d2'] = np.array([0.1, 0.02, 1.])
b_param['s'] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 60.
E["K"] = -90.
E["Leak"] = -60.
# model currents
currents_included["Na"] = True
currents_included["Kd"] =True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["L"] = True
currents_included["T"] = True
currents_included["Ca_K"] = True
currents_included["A"] =True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 49. * surf_area
g["Kd"] = 27. * (1-Kv_ratio)* surf_area
if currents_included["Kv_mut"]==True:
g["Kv"] = 27. * Kv_ratio/2 * surf_area
else:
g["Kv"] = 27. * Kv_ratio/2 * surf_area * 2
g["Kv_mut"] = 27. * Kv_ratio/2 * surf_area
g["A"] = 5. * surf_area
g["L"] = 5. * surf_area
g["T"] = 5. * surf_area
g["Ca_K"] = 1 * surf_area
g["Leak"] = 0.035 * surf_area
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
mut = 'STN_Delta_Kv'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
if not os.path.isdir(folder):
os.makedirs(folder)
STN_Kv_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)

View File

@ -0,0 +1,105 @@
import numpy as np
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.STN_fxns_Kv import STN_Kv_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.005
number_steps = 100
initial_period = 1000
num_gating = 17
num_current = 9
d = 61.4
r = d/2 * 10**-6 # radius in meters
vol = np.pi * r**3 # volume in liters
C, surf_area = capacitance(d, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'a', 'b', 'c', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut','Ca_conc', 'E_Ca']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = {}
i = 0
for var in np.array(['m', 'h', 'n', 'a', 'b', 'c', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut', 'Ca_conc','E_Ca']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'L', 'Kv', 'Kv_mut', 'T', 'Leak', 'Ca_K']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'] = np.array([-40., -1., 1.])
b_param['h'] = np.zeros(4)
b_param['h'] = np.array([-45.5, 1., 1., 0.05])
b_param['n'] = np.array([-41., -1., 1.])
b_param['p'] = np.array([-56., -1., 1.])
b_param['q'] = np.array([-85., 1., 1.])
b_param['r'] = np.array([0.17, -0.08, 1.])
b_param['a'] = np.array([-45., -14.7, 1.])
b_param['b'] = np.array([-90., 7.5, 1.])
b_param['c'] = np.array([-30.6, -5., 1.])
b_param['d1'] = np.array([-60, 7.5, 1.])
b_param['d2'] = np.array([0.1, 0.02, 1.])
b_param['s'] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 60.
E["K"] = -90.
E["Leak"] = -60.
# model currents
currents_included["Na"] = True
currents_included["Kd"] =True
currents_included["Kv"] = True
currents_included["Kv_mut"] = True
currents_included["L"] = True
currents_included["T"] = True
currents_included["Ca_K"] = True
currents_included["A"] =True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 49. * surf_area
g["Kd"] = 27. * (1-Kv_ratio)* surf_area
if currents_included["Kv_mut"]==True:
g["Kv"] = 27. * Kv_ratio/2 * surf_area
else:
g["Kv"] = 27. * Kv_ratio/2 * surf_area * 2
g["Kv_mut"] = 27. * Kv_ratio/2 * surf_area
g["A"] = 5. * surf_area
g["L"] = 5. * surf_area
g["T"] = 5. * surf_area
g["Ca_K"] = 1 * surf_area
g["Leak"] = 0.035 * surf_area
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
mut = 'STN_Kv'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
if not os.path.isdir(folder):
os.makedirs(folder)
STN_Kv_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)

View File

@ -0,0 +1,94 @@
import numpy as np
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.STN_fxns import STN_mut
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.005
number_steps = 100
initial_period = 1000
num_gating = 15
num_current = 8
d = 61.4
r = d/2 * 10**-6 # radius in meters
vol = np.pi * r**3 # volume in liters
C, surf_area = capacitance(d, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'a', 'b', 'a_mut', 'b_mut', 'c', 'd1', 'd2', 'p', 'q', 'r', 'Ca_conc', 'E_Ca']))
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = {}
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', 'L', 'A_mut', 'T', 'Leak', 'Ca_K']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'] = np.array([-40., -1. , 1.])
b_param['h'] = np.zeros(4)
b_param['h'] = np.array([-45.5, 1. , 1., 0.05])
b_param['n'] = np.array([-41., -1., 1.])
b_param['p'] = np.array([-56., -1., 1.])
b_param['q'] = np.array([-85., 1., 1.])
b_param['r'] = np.array([0.17, -0.08, 1.])
b_param['a'] = np.array([-45., -14.7, 1.])
b_param['b'] = np.array([-90., 7.5, 1.])
b_param['a_mut'] = np.array([-45., -14.7, 1.])
b_param['b_mut'] = np.array([-90., 7.5, 1.])
b_param['c'] = np.array([-30.6, -5., 1.])
b_param['d1'] = np.array([-60, 7.5, 1.])
b_param['d2'] = np.array([0.1, 0.02, 1.])
# reversal potentials
E["Na"] = 60.
E["K"] = -90.
E["Leak"] = -60.
# model currents
currents_included["Na"] = True
currents_included["Kd"] =True
currents_included["L"] = True
currents_included["T"] = True
currents_included["Ca_K"] = True
currents_included["A"] =True
currents_included["A_mut"] =True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 49. * surf_area
g["Kd"] = 27. * (1-Kv_ratio)* surf_area
g["A"] = 5. * surf_area *0.5
g["A_mut"] = 5. * surf_area *0.5
g["L"] = 5. * surf_area
g["T"] = 5. * surf_area
g["Ca_K"] = 1 * surf_area
g["Leak"] = 0.035 * surf_area
V_init = -70
prominence = 50
min_spike_height = 0
desired_AUC_width = high/5
folder = '../Neuron_models'
mut = 'STN'
mutations = {mut:{'g_ratio': 1, 'activation_Vhalf_diff': 0, 'activation_k_ratio':0}}
if not os.path.isdir(folder):
os.makedirs(folder)
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)

20
Code/README.md Normal file
View File

@ -0,0 +1,20 @@
### Python code repository for Koch et al. 2022
Modelling on the effects of ion current property changes on firing behaviour in diverse neuronal models and in with Kv1.1 mutations that cause episodic ataxia type 1.
## ./Functions
- Python function files for models
- from Pospischil et al. ()
- from Alexander et al. () with and without added Kv1.1
- from Otsuka et al () with and without added Kv1.1
- Utility_fxns for common useful functions
## ./Neuron_models
- Scripts for simulation of each model to current steps and ramp currents
## ./Sensitivity_Analysis
- Scripts for sensitivity analysis in each model
## ./KCNA1_mutations
- Scripts for simulation of KCNA1 mutations described in Lauxmann et al. (2021)

View File

@ -0,0 +1,110 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 4 08:24:33 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Cb_stellate_fxns import SA_Cb_stellate
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 9
num_current = 6
C, surf_area = capacitance(61.4, 1.50148)
variable = np.array(['m', 'h', 'n', 'n_A', 'h_A', 'n_A_mut', 'h_A_mut', 'm_T', 'h_T'])
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'n_A', 'h_A', 'n_A_mut', 'h_A_mut', 'm_T', 'h_T']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'n_A', 'h_A', 'n_A_mut', 'h_A_mut', 'm_T', 'h_T']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'A_mut', 'T', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-37., -3, 1])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-40., 4., 1., 0])
b_param['n'][:] = np.array([-23, -5, 1])
b_param['n_A'][:] = np.array([-27, -13.2, 1.])
b_param['h_A'][:] = np.array([-80., 6.5, 1.])
b_param['n_A_mut'][:] = np.array([-27, -13.2, 1.])
b_param['h_A_mut'][:] = np.array([-80., 6.5, 1.])
b_param['m_T'][:] = np.array([-50., -3, 1.])
b_param['h_T'][:] = np.array([-68., 3.75, 1.])
# reversal potentials
E["Na"] = 55.
E["K"] = -80.
E["Ca"] = 22.
E["Leak"] = -70. # as per Molineux et al 2005 and NOT the -38 in Alexander et al 2019
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["A_mut"] = False
currents_included["A"] = True
currents_included["T"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 3.4 * surf_area
g["Kd"] = 9.0556 * surf_area
g["A_mut"] = 0.
g["A"] = 15.0159 * surf_area
g["T"] = 0.45045 * surf_area
g["Leak"] = 0.07407 * surf_area
# save folder
folder = './SA/Cb_stellate'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n', 'n_A', 'h_A'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'A', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Cb_stellate)(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, lin_array, log_array, alt_types, alt_ind, alt)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1,116 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 4 08:24:33 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Cb_stellate_fxns_Kv import SA_Cb_stellate_Kv
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 11
num_current = 7
C, surf_area = capacitance(61.4, 1.50148)
variable = np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut'])
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'Kv', 'Kv_mut', 'T', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-37., -3, 1])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-40., 4., 1., 0])
b_param['n'][:] = np.array([-23, -5, 1])
b_param['n_A'][:] = np.array([-27, -13.2, 1.])
b_param['h_A'][:] = np.array([-80., 6.5, 1.])
b_param['m_T'][:] = np.array([-50., -3, 1.])
b_param['h_T'][:] = np.array([-68., 3.75, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 55.
E["K"] = -80.
E["Ca"] = 22.
E["Leak"] = -70. # as per Molineux et al 2005 and NOT the -38 in Alexander et al 2019
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = False
currents_included["A"] = False
currents_included["T"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 34 * surf_area
g["Kd"] = 9.0556 * surf_area
if currents_included["Kv_mut"] == True:
g["Kv"] = 1.50159 / 2 * surf_area
else:
g["Kv"] = 1.50159 / 2 * surf_area * 2
g["Kv_mut"] = 1.50159 / 2 * surf_area
g["A"] = 0 * surf_area
g["T"] = 0.45045 * surf_area
g["Leak"] = 0.07407 * surf_area
# save folder
folder = './SA/Cb_stellate_Kv_only'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n', 's', 'u'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'Kv', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Cb_stellate_Kv)(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, lin_array, log_array, alt_types, alt_ind, alt)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1,139 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 4 08:24:33 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.Cb_stellate_fxns_Kv import SA_Cb_stellate_Kv
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 11
num_current = 7
C, surf_area = capacitance(61.4, 1.50148)
variable = np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut'])
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'n_A', 'h_A', 'm_T', 'h_T', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'Kv', 'Kv_mut', 'T', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-37., -3, 1])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-40., 4., 1., 0])
b_param['n'][:] = np.array([-23, -5, 1])
b_param['n_A'][:] = np.array([-27, -13.2, 1.])
b_param['h_A'][:] = np.array([-80., 6.5, 1.])
b_param['m_T'][:] = np.array([-50., -3, 1.])
b_param['h_T'][:] = np.array([-68., 3.75, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 55.
E["K"] = -80.
E["Ca"] = 22.
E["Leak"] = -70. # as per Molineux et al 2005 and NOT the -38 in Alexander et al 2019
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = False
currents_included["A"] = True
currents_included["T"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 3.4 * surf_area
g["Kd"] = 9.0556 * (1-Kv_ratio) * surf_area
g["Kv"] = 9.0556 * Kv_ratio * surf_area
g["Kv_mut"] = 0.
g["A"] = 15.0159 * surf_area
g["T"] = 0.45045 * surf_area
g["Leak"] = 0.07407 * surf_area
# save folder
folder = './Sensitivity_Analysis/Cb_stellate_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n', 'n_A', 'h_A', 's', 'u'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'A', 'Kv', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Cb_stellate_Kv)(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, lin_array, log_array, alt_types, alt_ind, alt)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))
# #%% Get pd Dataframes for certain variables
# import pandas as pd
# AUC = pd.DataFrame(columns=mutations.keys())
# AUC_rel = pd.DataFrame(columns=mutations.keys())
# rheobase = pd.DataFrame(columns=mutations.keys())
# rheobase_fit = pd.DataFrame(columns=mutations.keys())
# rheobase_null_fit = pd.DataFrame(columns=mutations.keys())
# for mut in list(mutations.keys()):
# fname = os.path.join(folder, "{}.hdf5".format(mut))
# f = h5py.File(fname, "r")
# AUC['{}'.format(mut)] = f['analysis']['AUC']
# AUC_rel['{}'.format(mut)] = f['analysis']['AUC_rel']
# rheobase['{}'.format(mut)] = f['analysis']['rheobase']
# rheobase_fit['{}'.format(mut)] = f['analysis']['rheobase_fit']
# rheobase_null_fit['{}'.format(mut)] = f['analysis']['rheobase_null_fit']
# AUC.to_json(os.path.join(folder, 'AUC.json'))
# AUC_rel.to_json(os.path.join(folder, 'AUC_rel.json'))
# rheobase.to_json(os.path.join(folder, 'rheobase.json'))
# rheobase_fit.to_json(os.path.join(folder, 'rheobase_fit.json'))
# rheobase_null_fit.to_json(os.path.join(folder, 'rheobase_null_fit.json'))

View File

@ -0,0 +1,120 @@
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 5 18:14:31 2021
@author: nils
Original Pospischil FS Sensitivity analysis (no Kv1.1)
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict, NumpyEncoder
from Code.Functions.Pospischil_fxns import SA_Pospischil
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(56.9, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 502
V_init = -70
V_T = -57.9 # Threshold
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# create dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array([-45., -9.9998807337, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.4
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = False
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 58. * surf_area
g["Kd"] = 3.9 * surf_area
g["M"] = 0.075 * surf_area
g["Kv"] = 0.
g["Kv_mut"] = 0.
g["L"] = 0.
g["Leak"] = 0.038 * surf_area
# folder to save to
folder = './Sensitivity_Analysis/FS'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Pospischil)(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift, scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, folder, high, low,
number_steps, initial_period, sec, lin_array, log_array, alt_types, alt_ind, alt)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1,153 @@
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 5 18:14:31 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Utility import capacitance, stimulus_init, init_dict, NumpyEncoder
from Code.Functions.Pospischil_fxns import SA_Pospischil
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(56.9, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 502
V_init = -70
V_T = -57.9
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array([-45., -9.9998807337, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# mut_act_Vhalf_wt = -30.01851851851851
# mut_act_k_wt = -7.7333333333333325
# s_diff_Vhalf = mut_act_Vhalf_wt - b_param['s'][0]
# s_diff_k = mut_act_k_wt - b_param['s'][1]
# b_param['s'][1] = b_param['s'][1] + s_diff_k
# b_param['u'][1] = b_param['u'][1] + s_diff_k
# b_param['s'][0] = b_param['s'][0] + s_diff_Vhalf
# b_param['u'][0] = b_param['u'][0] + s_diff_Vhalf
# b_param['s_mut'][1] = b_param['s_mut'][1] + s_diff_k
# b_param['u_mut'][1] = b_param['u_mut'][1] + s_diff_k
# b_param['s_mut'][0] = b_param['s_mut'][0] + s_diff_Vhalf
# b_param['u_mut'][0] = b_param['u_mut'][0] + s_diff_Vhalf
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.4
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.10
g["Na"] = 58. * surf_area
g["Kd"] = 3.9 * (1 - Kv_ratio) * surf_area
g["M"] = 0.075 * surf_area
g["Kv"] = 3.9 * Kv_ratio * surf_area
g["Kv_mut"] = 0.
g["L"] = 0. * surf_area
g["Leak"] = 0.038 * surf_area
# save folder
folder = './Sensitivity_Analysis/FS_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n', 's', 'u'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'Kv', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,20, base=2)
log_array = np.append(log_array,1)
# %% multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Pospischil)(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift,
scale, b_param, slope_shift, gating, current, prominence, desired_AUC_width, folder, high,
low, number_steps, initial_period, sec, lin_array, log_array, alt_types, alt_ind, alt)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))
# #%% Get pd Dataframes for certain variables
# import pandas as pd
# AUC = pd.DataFrame(columns=mutations.keys())
# AUC_rel = pd.DataFrame(columns=mutations.keys())
# rheobase = pd.DataFrame(columns=mutations.keys())
# rheobase_fit = pd.DataFrame(columns=mutations.keys())
# rheobase_null_fit = pd.DataFrame(columns=mutations.keys())
# for mut in list(mutations.keys()):
# fname = os.path.join(folder, "{}.hdf5".format(mut))
# f = h5py.File(fname, "r")
# AUC['{}'.format(mut)] = f['analysis']['AUC']
# AUC_rel['{}'.format(mut)] = f['analysis']['AUC_rel']
# rheobase['{}'.format(mut)] = f['analysis']['rheobase']
# rheobase_fit['{}'.format(mut)] = f['analysis']['rheobase_fit']
# rheobase_null_fit['{}'.format(mut)] = f['analysis']['rheobase_null_fit']
# AUC.to_json(os.path.join(folder, 'AUC.json'))
# AUC_rel.to_json(os.path.join(folder, 'AUC_rel.json'))
# rheobase.to_json(os.path.join(folder, 'rheobase.json'))
# rheobase_fit.to_json(os.path.join(folder, 'rheobase_fit.json'))
# rheobase_null_fit.to_json(os.path.join(folder, 'rheobase_null_fit.json'))

View File

@ -0,0 +1,118 @@
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 5 18:14:31 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict, NumpyEncoder
from Code.Functions.Pospischil_fxns import SA_Pospischil
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.00035
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.8, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 934
V_init = -70
V_T = -67.9
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array([-45., -9.9998807337, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -56.2
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = False
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 10. * surf_area
g["Kd"] = 2.1 * surf_area
g["M"] = 0.0098 * surf_area
g["Kv"] = 0.
g["Kv_mut"] = 0.
g["L"] = 0.
g["Leak"] = 0.0205 * surf_area
# save folder
folder = './Sensitivity_Analysis/RS_inhib'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Pospischil)(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift,
scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, folder, high, low,
number_steps, initial_period, sec, lin_array, log_array, alt_types, alt_ind, alt)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1,116 @@
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 5 18:14:31 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict, NumpyEncoder
from Code.Functions.Pospischil_fxns import SA_Pospischil
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.00035
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.8, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 934
V_init = -70
V_T = -67.9
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'][:] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'][:] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'][:] = np.array([-57.37, 20.98, 1.])
b_param['p'][:] = np.array([-45., -9.9998807337, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -56.2
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 10. * surf_area
g["Kd"] = 2.1 * (1 - Kv_ratio) * surf_area
g["M"] = 0.0098 * surf_area
g["Kv"] = 2.1 * Kv_ratio * surf_area
g["Kv_mut"] = 0.
g["L"] = 0. * surf_area
g["Leak"] = 0.0205 * surf_area
# save folder
folder = './Sensitivity_Analysis/RS_inhib_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n', 's', 'u'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'Kv', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Pospischil)(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift,
scale, b_param, slope_shift, gating, current, prominence, desired_AUC_width, folder, high,
low, number_steps, initial_period, sec, lin_array, log_array, alt_types, alt_ind, alt)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1,100 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 4 08:24:33 2021
@author: nils
"""
import numpy as np
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict, NumpyEncoder
from Code.Functions.Pospischil_fxns import SA_Pospischil
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.4, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 608
V_init = -70
V_T = -56.2
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# gating parameters
b_param = {}
b_param['m'] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'] = np.array([-57.37, 20.98, 1.0])
b_param['p'] = np.array([-45., -9.9998807337, 1.])
b_param['s'] = np.array([-14.16, -10.15, 1.0])
b_param['u'] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.array([-31., 5.256, 1., 0.245])
# reversal potential
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.3
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = False
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
g["Na"] = 56. * surf_area
g["Kd"] = 6. * surf_area
g["M"] = 0.075 * surf_area
g["Kv"] = 0.
g["Kv_mut"] = 0.
g["L"] = 0.
g["Leak"] = 0.0205 * surf_area
folder = './Sensitivity_Analysis/RS_pyramidal'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Pospischil)(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift,
scale,
b_param, slope_shift, gating, current, prominence, desired_AUC_width, folder, high, low,
number_steps, initial_period, sec, lin_array, log_array, alt_types, alt_ind, alt)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1,115 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 4 08:24:33 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict, NumpyEncoder
from Code.Functions.Pospischil_fxns import SA_Pospischil
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 10
num_current = 7
C, surf_area = capacitance(61.4, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(
np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']))
tau_max_p = 608
V_init = -70
V_T = -56.2
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# create dictionary
ind_dict = Dict.empty(key_type=types.unicode_type, value_type=types.int64, )
i = 0
for var in np.array(['m', 'h', 'n', 'q', 'r', 'p', 's', 'u', 's_mut', 'u_mut']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'M', 'Kv', 'Kv_mut', 'L', 'Leak']):
ind_dict[var] = i
i += 1
# gating parameters
b_param = {}
b_param['m'] = np.array([-34.33054521, -8.21450277, 1.42295686])
b_param['h'] = np.array([-34.51951036, 4.04059373, 1., 0.05])
b_param['n'] = np.array([-63.76096946, -13.83488194, 7.35347425])
b_param['q'] = np.array([-39.03684525, -5.57756176, 2.25190197])
b_param['r'] = np.array([-57.37, 20.98, 1.0])
b_param['p'] = np.array([-45., -9.9998807337, 1.])
b_param['s'] = np.array([-14.16, -10.15, 1.0])
b_param['u'] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 50.
E["K"] = -90.
E["Ca"] = 120.
E["Leak"] = -70.3
# model currents
currents_included["Na"] = True
currents_included["Kd"] = True
currents_included["Kv"] = True
currents_included["Kv_mut"] = False
currents_included["L"] = False
currents_included["M"] = True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.1
g["Na"] = 56. * surf_area
g["Kd"] = 6. * (1 - Kv_ratio) * surf_area
g["M"] = 0.075 * surf_area
g["Kv"] = 6. * Kv_ratio * surf_area
g["Kv_mut"] = 0.
g["L"] = 0. * surf_area
g["Leak"] = 0.0205 * surf_area
# save folder
folder = './Sensitivity_Analysis/RS_pyramidal_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n', 's', 'u'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'Kv', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Pospischil)(V_init, V_T, g, E, I_in, dt, currents_included, stim_time, stim_num, C, tau_max_p, shift,
scale, b_param, slope_shift, gating, current, prominence, desired_AUC_width, folder, high,
low, number_steps, initial_period, sec, lin_array, log_array, alt_types, alt_ind, alt)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1,116 @@
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 5 19:44:04 2021
@author: nils
original STN model SA
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.STN_fxns import SA_STN
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 15
num_current = 8
d = 61.4
r = d/2 * 10**-6 # radius in meters
vol = np.pi * r**3 # volume in liters
C, surf_area = capacitance(d, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'a', 'b', 'a_mut', 'b_mut', 'c', 'd1', 'd2', 'p', 'q', 'r', 'Ca_conc', 'E_Ca']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
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', 'L', 'A_mut', 'T', 'Leak', 'Ca_K']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-40., -1., 1.])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-45.5, 1., 1., 0.05])
b_param['n'][:] = np.array([-41., -1., 1.])
b_param['p'][:] = np.array([-56., -1., 1.])
b_param['q'][:] = np.array([-85., 1., 1.])
b_param['r'][:] = np.array([0.17, -0.08, 1.])
b_param['a'][:] = np.array([-45., -14.7, 1.])
b_param['b'][:] = np.array([-90., 7.5, 1.])
b_param['a_mut'][:] = np.array([-45., -14.7, 1.])
b_param['b_mut'][:] = np.array([-90., 7.5, 1.])
b_param['c'][:] = np.array([-30.6, -5., 1.])
b_param['d1'][:] = np.array([-60, 7.5, 1.])
b_param['d2'][:] = np.array([0.1, 0.02, 1.])
# reversal potentials
E["Na"] = 60.
E["K"] = -90.
E["Leak"] = -60.
# model currents
currents_included["Na"] = True
currents_included["Kd"] =True
currents_included["L"] = True
currents_included["T"] = True
currents_included["Ca_K"] = True
currents_included["A"] =True
currents_included["A_mut"] =False
currents_included["Leak"] = True
# model conductances
g["Na"] = 49. * surf_area
g["Kd"] = 57. * surf_area
g["A"] = 5. * surf_area
g["A_mut"] = 0.
g["L"] = 5. * surf_area
g["T"] = 5. * surf_area
g["Ca_K"] = 1 * surf_area
g["Leak"] = 0.035 * surf_area
# save folder
folder = '../Sensitivity_Analysis/STN'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n', 'a', 'b'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'A', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% run SA with multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(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)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1,129 @@
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 5 19:44:04 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.STN_fxns_Kv import SA_Kv_STN
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 17
num_current = 9
d = 61.4
r = d/2 * 10**-6 # radius in meters
vol = np.pi * r**3 # volume in liters
C, surf_area = capacitance(d, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'a', 'b', 'c', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut','Ca_conc', 'E_Ca']))
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
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', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut', 'Ca_conc','E_Ca']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'L', 'Kv', 'Kv_mut', 'T', 'Leak', 'Ca_K']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-40., -1. , 1.])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-45.5, 1. , 1., 0.05])
b_param['n'][:] = np.array([-41., -1., 1.])
b_param['p'][:] = np.array([-56., -1., 1.])
b_param['q'][:] = np.array([-85., 1., 1.])
b_param['r'][:] = np.array([0.17, -0.08, 1.])
b_param['a'][:] = np.array([-45., -14.7, 1.])
b_param['b'][:] = np.array([-90., 7.5, 1.])
b_param['c'][:] = np.array([-30.6, -5., 1.])
b_param['d1'][:] = np.array([-60, 7.5, 1.])
b_param['d2'][:] = np.array([0.1, 0.02, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 60.
E["K"] = -90.
E["Leak"] = -60.
# model currents
currents_included["Na"] = True
currents_included["Kd"] =True
currents_included["Kv"] = True
currents_included["Kv_mut"] = False
currents_included["L"] = True
currents_included["T"] = True
currents_included["Ca_K"] = True
currents_included["A"] =True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.01
g["Na"] = 49. * surf_area
g["Kd"] = 57. * surf_area
g["Kv"] = 0.5 * surf_area
g["Kv_mut"] = 0
g["A"] = 0. * surf_area
g["L"] = 5. * surf_area
g["T"] = 5. * surf_area
g["Ca_K"] = 1 * surf_area
g["Leak"] = 0.035 * surf_area
V_init = -70
# save folder
folder = './Sensitivity_Analysis/STN_Kv_only'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n', 's','u'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'Kv', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% run SA with multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Kv_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)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1,126 @@
# -*- coding: utf-8 -*-
"""
Created on Sat Jun 5 19:44:04 2021
@author: nils
"""
import numpy as np
from numba import types
from numba.typed import Dict
from joblib import Parallel, delayed
import os
from Code.Functions.Utility import capacitance, stimulus_init, init_dict
from Code.Functions.STN_fxns_Kv import SA_Kv_STN
# model parameters
dt = 0.01
sec = 2
low = 0
high = 0.001
number_steps = 200
initial_period = 1000
num_gating = 17
num_current = 9
d = 61.4
r = d/2 * 10**-6 # radius in meters
vol = np.pi * r**3 # volume in liters
C, surf_area = capacitance(d, 1)
stim_time, I_in, stim_num, V_m = stimulus_init(low, high, number_steps, initial_period, dt, sec)
shift, scale, slope_shift, E, currents_included, b_param, g = init_dict(np.array(['m', 'h', 'n', 'a', 'b', 'c', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut','Ca_conc', 'E_Ca']))
V_init = -70
# initialize arrays
current = np.zeros((num_current, stim_num))
gating = np.zeros((num_gating, stim_num))
# initialize dictionary
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', 'd1', 'd2', 'p', 'q', 'r', 's', 'u', 's_mut', 'u_mut', 'Ca_conc','E_Ca']):
ind_dict[var] = i
i += 1
i = 0
for var in np.array(['Na', 'Kd', 'A', 'L', 'Kv', 'Kv_mut', 'T', 'Leak', 'Ca_K']):
ind_dict[var] = i
i += 1
# gating parameters
b_param['m'][:] = np.array([-40., -1. , 1.])
b_param['h'] = np.zeros(4)
b_param['h'][:] = np.array([-45.5, 1. , 1., 0.05])
b_param['n'][:] = np.array([-41., -1., 1.])
b_param['p'][:] = np.array([-56., -1., 1.])
b_param['q'][:] = np.array([-85., 1., 1.])
b_param['r'][:] = np.array([0.17, -0.08, 1.])
b_param['a'][:] = np.array([-45., -14.7, 1.])
b_param['b'][:] = np.array([-90., 7.5, 1.])
b_param['c'][:] = np.array([-30.6, -5., 1.])
b_param['d1'][:] = np.array([-60, 7.5, 1.])
b_param['d2'][:] = np.array([0.1, 0.02, 1.])
b_param['s'][:] = np.array([-14.16, -10.15, 1.])
b_param['u'] = np.zeros(4)
b_param['u'][:] = np.array([-31., 5.256, 1., 0.245])
b_param['s_mut'][:] = np.array([-14.16, -10.15, 1.])
b_param['u_mut'] = np.zeros(4)
b_param['u_mut'][:] = np.array([-31., 5.256, 1., 0.245])
# reversal potentials
E["Na"] = 60.
E["K"] = -90.
E["Leak"] = -60.
# model currents
currents_included["Na"] = True
currents_included["Kd"] =True
currents_included["Kv"] = True
currents_included["Kv_mut"] = False
currents_included["L"] = True
currents_included["T"] = True
currents_included["Ca_K"] = True
currents_included["A"] =True
currents_included["Leak"] = True
# model conductances
Kv_ratio = 0.01
g["Na"] = 49. * surf_area
g["Kd"] = 57. * (1-Kv_ratio)* surf_area
g["Kv"] = 57. * Kv_ratio * surf_area
g["Kv_mut"] = 0
g["A"] = 5. * surf_area
g["L"] = 5. * surf_area
g["T"] = 5. * surf_area
g["Ca_K"] = 1 * surf_area
g["Leak"] = 0.035 * surf_area
# save folder
folder = '../Sensitivity_Analysis/STN_Kv'
if not os.path.isdir(folder):
os.makedirs(folder)
#%% setup for one-factor-at-a-time SA
var = np.array(['m', 'h', 'n', 'a', 'b', 's','u'])
type_names = np.append(np.array(['shift' for i in range(var.shape[0])]),
np.array(['slope' for i in range(var.shape[0])]))
cur = np.array(['Na', 'Kd', 'A','Kv', 'Leak'])
type_names = np.append(type_names, np.array(['g' for i in range(cur.shape[0])]))
var = np.append(var, var)
var = np.append(var, cur)
alt_types = np.c_[var, type_names]
lin_array = np.arange(-10, 11, 1)
log_array = np.logspace(-1,1,21, base=2)
# %% run SA with multiprocessing
prominence = 50
desired_AUC_width = high/5
Parallel(n_jobs=8, verbose=9)(
delayed(SA_Kv_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)
for alt_ind in range(alt_types.shape[0]) for alt in range(21))

View File

@ -0,0 +1 @@
{"I407M": {"activation_Vhalf_diff": 44.2, "activation_k_ratio": 1.0, "g_ratio": 0.43523316100000004}, "G311D": {"activation_Vhalf_diff": 0.0, "activation_k_ratio": 1.0, "g_ratio": 0.333333333}, "V408A 1": {"activation_Vhalf_diff": 1.0, "activation_k_ratio": 1.0, "g_ratio": 0.658536585}, "V404I": {"activation_Vhalf_diff": 11.3, "activation_k_ratio": 1.25, "g_ratio": 1.0}, "V408A 3": {"activation_Vhalf_diff": 0.9, "activation_k_ratio": 1.228915663, "g_ratio": 1.0}, "F249C": {"activation_Vhalf_diff": 0.0, "activation_k_ratio": 1.0, "g_ratio": 0.001}, "F414S": {"activation_Vhalf_diff": 14.9, "activation_k_ratio": 1.0, "g_ratio": 0.079493088}, "F184C 2": {"activation_Vhalf_diff": 25.9, "activation_k_ratio": 0.87654321, "g_ratio": 0.151}, "E325D 1": {"activation_Vhalf_diff": 60.4, "activation_k_ratio": 2.160493827, "g_ratio": 0.038}, "G311S": {"activation_Vhalf_diff": 31.7, "activation_k_ratio": 1.654320988, "g_ratio": 0.22899999999999998}, "V408A 2": {"activation_Vhalf_diff": 0.4, "activation_k_ratio": 0.9259259259999999, "g_ratio": 0.498}, "E325D 2": {"activation_Vhalf_diff": 21.8, "activation_k_ratio": 1.518072289, "g_ratio": 1.0}, "C185W 1": {"activation_Vhalf_diff": 4.4, "activation_k_ratio": 1.0, "g_ratio": 0.298791019}, "A242P": {"activation_Vhalf_diff": -4.8, "activation_k_ratio": 0.947368421, "g_ratio": 0.1}, "T226R": {"activation_Vhalf_diff": 0.0, "activation_k_ratio": 1.0, "g_ratio": 1.0}, "R307C": {"activation_Vhalf_diff": 11.4, "activation_k_ratio": 1.0, "g_ratio": 0.254608295}, "T226M": {"activation_Vhalf_diff": 14.8, "activation_k_ratio": 1.0, "g_ratio": 0.039}, "V408L": {"activation_Vhalf_diff": 0.0, "activation_k_ratio": 1.0, "g_ratio": 1.0}, "R167M": {"activation_Vhalf_diff": 10.0, "activation_k_ratio": 1.0, "g_ratio": 0.283246978}, "wt": {"activation_Vhalf_diff": 0.0, "activation_k_ratio": 1.0, "g_ratio": 1.0}, "F249I": {"activation_Vhalf_diff": -0.6, "activation_k_ratio": 0.7654320990000001, "g_ratio": 0.01}, "N255D": {"activation_Vhalf_diff": 31.0, "activation_k_ratio": 1.0, "g_ratio": 0.180904523}, "E283K": {"activation_Vhalf_diff": 9.3, "activation_k_ratio": 1.054054054, "g_ratio": 0.666666667}, "F414C": {"activation_Vhalf_diff": 2.2, "activation_k_ratio": 1.216666667, "g_ratio": 1.0}, "R324T": {"activation_Vhalf_diff": 0.0, "activation_k_ratio": 1.0, "g_ratio": 1.0}, "R239S": {"activation_Vhalf_diff": 0.0, "activation_k_ratio": 1.0, "g_ratio": 1.0}, "F303V": {"activation_Vhalf_diff": 36.0, "activation_k_ratio": 1.333333333, "g_ratio": 0.035714286}, "C185W 2": {"activation_Vhalf_diff": 0.8, "activation_k_ratio": 0.739130435, "g_ratio": 1.0}, "F184C 1": {"activation_Vhalf_diff": 28.0, "activation_k_ratio": 1.0, "g_ratio": 0.39024390200000003}, "I262M": {"activation_Vhalf_diff": 0.0, "activation_k_ratio": 1.0, "g_ratio": 0.075268817}, "R417Stop": {"activation_Vhalf_diff": 8.5, "activation_k_ratio": 1.684210526, "g_ratio": 0.02}, "V174F": {"activation_Vhalf_diff": 35.1, "activation_k_ratio": 0.75308642, "g_ratio": 0.076}, "I262T": {"activation_Vhalf_diff": 33.6, "activation_k_ratio": 2.136363636, "g_ratio": 0.3}, "I177N": {"activation_Vhalf_diff": -2.8, "activation_k_ratio": 1.226666667, "g_ratio": 1.0}, "T226A": {"activation_Vhalf_diff": 14.3, "activation_k_ratio": 0.740740741, "g_ratio": 0.024}}

View File

@ -112,17 +112,17 @@ Nils A. Koch\textsuperscript{1,2}, Lukas Sonnenberg\textsuperscript{1,2}, Ulrike
\subsection{Author Contributions - Each author must be identified with at least one of the following: Designed research, Performed research, Contributed unpublished reagents/ analytic tools, Analyzed data, Wrote the paper.} \subsection{Author Contributions - Each author must be identified with at least one of the following: Designed research, Performed research, Contributed unpublished reagents/ analytic tools, Analyzed data, Wrote the paper.}
\notenk{Adjust as you deem appropriate}\\ \notenk{Adjust as you deem appropriate}\\
NK, LS, JB Designed Research; NK, LS, UBSH, SL, JB Designed Research;
NK Performed research; NK Performed research;
NK, LS Analyzed data; NK, LS Analyzed data;
NK, LS, JB Wrote the paper NK, LS,UBSH, SL, JB Wrote the paper
\subsection{Correspondence should be addressed to (include email address)} \subsection{Correspondence should be addressed to (include email address)}
\ \notenk{Nils oder Jan?} \ \notenk{Nils oder Jan?}
\subsection{Number of Figures} \subsection{Number of Figures}
5 5
\subsection{Number of Tables} \subsection{Number of Tables}
2 3
\subsection{Number of Multimedia} \subsection{Number of Multimedia}
0 0
\subsection{Number of words for Abstract} \subsection{Number of words for Abstract}
@ -389,9 +389,9 @@ The effects of altered ion channel properties on firing is generally influenced
\bibliographystyle{eneuro} %dcu \bibliographystyle{eneuro} %dcu
\newpage %\newpage
\section*{Figure/Table/Extended Data Legends} \section*{Figures}
% \textit{Figures must be numbered independently of tables and multimedia and cited in the manuscript. Do not duplicate data by presenting it both in the text and in a figure. % \textit{Figures must be numbered independently of tables and multimedia and cited in the manuscript. Do not duplicate data by presenting it both in the text and in a figure.
% A title should be part of the legend and not lettered onto the figure. A legend must be included in the manuscript document after the reference list, and should include enough detail to be intelligible without reference to the text. Specific individuals contributions to data acquisition, analysis, or other responsibility resulting in a figure may be included at the end of each legend. Please use the heading “Figure Contributions” and state each contribution with the authors full name. % A title should be part of the legend and not lettered onto the figure. A legend must be included in the manuscript document after the reference list, and should include enough detail to be intelligible without reference to the text. Specific individuals contributions to data acquisition, analysis, or other responsibility resulting in a figure may be included at the end of each legend. Please use the heading “Figure Contributions” and state each contribution with the authors full name.
% Figure Contributions: John Smith performed the experiments; Jane Jones analyzed the data. % Figure Contributions: John Smith performed the experiments; Jane Jones analyzed the data.
@ -401,8 +401,53 @@ The effects of altered ion channel properties on firing is generally influenced
% Remove top and right borderlines that to not contain measuring metrics from all graph/histogram figure panels (i.e., do not box the panels in). Do not include any two-bar graphs/histograms; instead state those values in the text. % Remove top and right borderlines that to not contain measuring metrics from all graph/histogram figure panels (i.e., do not box the panels in). Do not include any two-bar graphs/histograms; instead state those values in the text.
% All illustrations documenting results must include a bar to indicate the scale. All labels used in a figure should be explained in the legend. The migration of protein molecular weight size markers or nucleic acid size markers must be indicated and labeled appropriately (e.g., “kD”, “nt”, “bp”) on all figure panels showing gel electrophoresis.} % All illustrations documenting results must include a bar to indicate the scale. All labels used in a figure should be explained in the legend. The migration of protein molecular weight size markers or nucleic acid size markers must be indicated and labeled appropriately (e.g., “kD”, “nt”, “bp”) on all figure panels showing gel electrophoresis.}
add from manuscript text before submission \setcounter{figure}{0}
%\setcounter{figure}{0}
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{Figures/diversity_in_firing.pdf}
\linespread{1.}\selectfont
\caption[]{Diversity in Neuronal Model Firing. Spike trains (left), frequency-current (fI) curves (right) for Cb stellate (A), RS inhibitory (B), FS (C), RS pyramidal (D), RS inhibitory +\Kv (E), Cb stellate +\Kv (F), FS +\Kv (G), RS pyramidal +\Kv (H), STN +\Kv (I), Cb stellate \(\Delta\)\Kv (J), STN \(\Delta\)\Kv (K), and STN (L) neuron models. Black marker on the fI curves indicate the current step at which the spike train occurs. The green marker indicates the current at which firing begins in response to an ascending current ramp, whereas the red marker indicates the current at which firing ceases in response to a descending current ramp (see \Cref{fig:ramp_firing}).}
\label{fig:diversity_in_firing}
\end{figure}
\begin{figure}[tp]
\centering
\includegraphics[width=0.5\linewidth]{Figures/firing_characterization_arrows.pdf}
\linespread{1.}\selectfont
\caption[]{Characterization of firing with AUC and rheobase. (A) The area under the curve (AUC) of the repetitive firing frequency-current (fI) curve. (B)
Changes in firing as characterized by \(\Delta\)AUC and \(\Delta\)rheobase occupy 4 quadrants separated by no changes in AUC and rheobase. Representative schematic fI curves in red with respect to a reference fI curve (blue) depict the general changes associated with each quadrant.}
\label{fig:firing_characterization}
\end{figure}
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{Figures/AUC_correlation.pdf}
\linespread{1.}\selectfont
\caption[]{Effects of altered channel kinetics on AUC in various neuron models. The fI curves corresponding to shifts in FS \(+\)\Kv model delayed rectifier K half activation \(V_{1/2}\) (A), changes \Kv activation slope factor \(k\) in the FS \(+\)\Kv model (B), and changes in maximal conductance of delayed rectifier K current in the STN \(+\)\Kv model (C) are shown. The \ndAUC of fI curves is plotted against delayed rectifier K half activation potential (\(\Delta V_{1/2}\); B), \Kv activation slope factor \(k\) (k/\(\textrm{k}_{WT}\); E) and maximal conductance \(g\) of the delayed rectifier K current (g/\(\textrm{g}_{WT}\); H) for all models (thin lines) with relationships from the fI curve examples (A, D, G respectively) highlighted by thick lines with colors corresponding to the box highlighting each set of fI curves. The Kendall rank correlation (Kendall \(\tau\)) coefficients between shifts in half maximal potential \(V_{1/2}\) and \ndAUC (C), slope factor k and \ndAUC (F) as well as maximal current conductances and \ndAUC (I) for each model and current property is computed. The relationships between \(\Delta V_{1/2}\), k/\(\textrm{k}_{WT}\), and g/\(\textrm{g}_{WT}\) and \ndAUC for the Kendall rank correlations highlighted in the black boxes are depicted in (B), (E) and (H) respectively.}
\label{fig:AUC_correlation}
\end{figure}
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{Figures/rheobase_correlation.pdf}
\linespread{1.}\selectfont
\caption[]{Effects of altered channel kinetics on rheobase. The fI curves corresponding to shifts in FS \(+\)\Kv model \Kv activation \(V_{1/2}\) (A), changes \Kv inactivation slope factor \(k\) in the Cb stellate \(+\)\Kv model (B), and changes in maximal conductance of the leak current in the Cb stellate model (C) are shown. The \drheo of fI curves is plotted against \Kv half activation potential (\(\Delta V_{1/2}\); B), \Kv inactivation slope factor \(k\) (k/\(\textrm{k}_{WT}\); E) and maximal conductance \(g\) of the leak current (g/\(\textrm{g}_{WT}\); H) for all models (thin lines) with relationships from the fI curve examples (A, D, G respectively) highlighted by thick lines with colors corresponding to the box highlighting each set of fI curves. The Kendall rank correlation (Kendall \(\tau\)) coefficients between shifts in half maximal potential \(V_{1/2}\) and \drheo (C), slope factor k and \drheo (F) as well as maximal current conductances and \drheo (I) for each model and current property is computed. The relationships between \(\Delta V_{1/2}\), k/\(\textrm{k}_{WT}\), and g/\(\textrm{g}_{WT}\) and \drheo for the Kendall rank correlations highlighted in the black boxes are depicted in (B), (E) and (H) respectively..}
\label{fig:rheobase_correlation}
\end{figure}
\begin{figure}[tp]
\centering
\includegraphics[width=\linewidth]{Figures/simulation_model_comparison.pdf}
\linespread{1.}\selectfont
\caption[]{Effects of episodic ataxia type~1 associated \Kv mutations on firing. Effects of \Kv mutations on AUC (percent change in normalized \(\Delta\)AUC) and rheobase (\(\Delta\)Rheobase) compared to wild type for RS pyramidal +\Kv (A), RS inhibitory +\Kv (B), FS +\Kv (C), Cb stellate (D), Cb stellate +\Kv (E), Cb stellate \(\Delta\)\Kv (F), STN (G), STN +\Kv (H) and STN \(\Delta\)\Kv (I) models. V174F, F414C, E283K, and V404I mutations are highlighted in color for each model. Pairwise Kendall rank correlation coefficients (Kendall \(\tau\)) between the effects of \Kv mutations on rheobase and on AUC are shown in J and K respectively. Marker shape is indicative of model/firing type, and grey dashed lines denote the quadrants of firing characterization (see \Cref{fig:firing_characterization}).}
\label{fig:simulation_model_comparision}
\end{figure}
%\captionof{figure}{Characterization of firing with AUC and rheobase. (A) The area under the curve (AUC) of the repetitive firing frequency-current (fI) curve. (B) Changes in firing as characterized by \(\Delta\)AUC and \(\Delta\)rheobase occupy 4 quadrants separated by no changes in AUC and rheobase. Representative schematic fI curves in blue with respect to a reference fI curve (black) depict the general changes associated with each quadrant.} %\captionof{figure}{Characterization of firing with AUC and rheobase. (A) The area under the curve (AUC) of the repetitive firing frequency-current (fI) curve. (B) Changes in firing as characterized by \(\Delta\)AUC and \(\Delta\)rheobase occupy 4 quadrants separated by no changes in AUC and rheobase. Representative schematic fI curves in blue with respect to a reference fI curve (black) depict the general changes associated with each quadrant.}
% %
%\captionof{figure}{Diversity in Neuronal Model Firing. Spike trains (left), frequency-current (fI) curves (right) for Cb stellate (A), RS inhibitory (B), FS (C), RS pyramidal (D), RS inhibitory +\Kv (E), Cb stellate +\Kv (F), FS +\Kv (G), RS pyramidal +\Kv (H), STN +\Kv (I), Cb stellate \(\Delta\)\Kv (J), STN \(\Delta\)\Kv (K), and STN (L) neuron models. Black marker on the fI curves indicate the current step at which the spike train occurs. The green marker indicates the current at which firing begins in response to an ascending current ramp, whereas the red marker indicates the current at which firing ceases in response to a descending current ramp.} %\captionof{figure}{Diversity in Neuronal Model Firing. Spike trains (left), frequency-current (fI) curves (right) for Cb stellate (A), RS inhibitory (B), FS (C), RS pyramidal (D), RS inhibitory +\Kv (E), Cb stellate +\Kv (F), FS +\Kv (G), RS pyramidal +\Kv (H), STN +\Kv (I), Cb stellate \(\Delta\)\Kv (J), STN \(\Delta\)\Kv (K), and STN (L) neuron models. Black marker on the fI curves indicate the current step at which the spike train occurs. The green marker indicates the current at which firing begins in response to an ascending current ramp, whereas the red marker indicates the current at which firing ceases in response to a descending current ramp.}
@ -414,17 +459,23 @@ add from manuscript text before submission
%\captionof{figure}{Effects of episodic ataxia type~1 associated \Kv mutations on firing. Effects of \Kv mutations on AUC (\(AUC_{contrast}\)) and rheobase (\(\Delta\)rheobase) compared to wild type for RS pyramidal +\Kv (A), RS inhibitory +\Kv (B), FS +\Kv (C), Cb stellate (D), Cb stellate +\Kv (E), Cb stellate \(\Delta\)\Kv (F), STN (G), STN +\Kv (H) and STN \(\Delta\)\Kv (I) models V174F, F414C, E283K, and V404I mutations are highlighted in color for each model. Pairwise Kendall rank correlation coefficients (Kendall \(\tau\)) between the effects of \Kv mutations on rheobase and on AUC are shown in J and K respectively.} %\captionof{figure}{Effects of episodic ataxia type~1 associated \Kv mutations on firing. Effects of \Kv mutations on AUC (\(AUC_{contrast}\)) and rheobase (\(\Delta\)rheobase) compared to wild type for RS pyramidal +\Kv (A), RS inhibitory +\Kv (B), FS +\Kv (C), Cb stellate (D), Cb stellate +\Kv (E), Cb stellate \(\Delta\)\Kv (F), STN (G), STN +\Kv (H) and STN \(\Delta\)\Kv (I) models V174F, F414C, E283K, and V404I mutations are highlighted in color for each model. Pairwise Kendall rank correlation coefficients (Kendall \(\tau\)) between the effects of \Kv mutations on rheobase and on AUC are shown in J and K respectively.}
%\newpage %\newpage
\subsection*{Tables} \section*{Tables}
% \textit{All tables must be numbered independently of figures, multimedia, and 3D models and cited in the manuscript. Do not duplicate data by presenting it both in the text and in a table. % \textit{All tables must be numbered independently of figures, multimedia, and 3D models and cited in the manuscript. Do not duplicate data by presenting it both in the text and in a table.
% Each table should include a title and legend; legends should be included in the manuscript file after the reference list. Legends should include sufficient detail to be intelligible without reference to the text and define all symbols and include essential information. % Each table should include a title and legend; legends should be included in the manuscript file after the reference list. Legends should include sufficient detail to be intelligible without reference to the text and define all symbols and include essential information.
% Each table should be double-spaced. Multiple-part tables (A and B sections with separate subtitles) should be avoided, especially when there are two [different] sets [or types] of column headings. % Each table should be double-spaced. Multiple-part tables (A and B sections with separate subtitles) should be avoided, especially when there are two [different] sets [or types] of column headings.
% Do not use color or shading, bold or italic fonts, or lines to highlight information. Indention of text and sometimes, additional space between lines is preferred. Tables with color or shading in the table body will need to be processed as a figure.} % Do not use color or shading, bold or italic fonts, or lines to highlight information. Indention of text and sometimes, additional space between lines is preferred. Tables with color or shading in the table body will need to be processed as a figure.}
\setcounter{table}{0} \setcounter{table}{0}
%\input{g_table} %add from manuscript text before submission
add from manuscript text before submission \input{g_table}
\input{gating_table}
\input{statistical_table}
\FloatBarrier
%\newpage %\newpage
\subsection*{Extended Data} \section*{Extended Data}
\beginsupplement \beginsupplement
\begin{figure}[tp]%described \begin{figure}[tp]%described
\centering \centering