# -*- 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'))