Subir archivos a 'functions'

functions/f_plots.py

@@ -0,0 +1,133 @@
+
+import matplotlib
+from matplotlib import rcParams
+matplotlib.rcParams['text.usetex'] = True
+import matplotlib.pyplot as plt
+import pylab as pl
+from copy import copy
+
+
+### Format for figures
+import json
+P= json.load(open("cfg/PLOSmpl.json")) # import plot parameter
+
+matplotlib.rcParams.update([(j,k) for (j,k) in P.items()
+                            if j in matplotlib.rcParams.keys()])
+
+matplotlib.rcParams['lines.linewidth']=0.065
+fig_wide=matplotlib.rcParams["figure.figsize"][0]
+fig_height=matplotlib.rcParams["figure.figsize"][1]
+
+
+### includes scale bar in the figures
+import matplotlib.offsetbox
+from matplotlib.lines import Line2D
+
+class AnchoredHScaleBarY(matplotlib.offsetbox.AnchoredOffsetbox):
+    """ size: length of bar in data units
+        extent : height of bar ends in axes units """
+    def __init__(self, size=1, extent = 0.03, label="", loc=2, ax=None,
+                 pad=0.4, borderpad=0.5, ppad = 0, sep=2, prop=None, 
+                 frameon=True, linekw={}, **kwargs):
+        if not ax:
+            ax = plt.gca()
+        trans = ax.get_yaxis_transform()
+        size_bar = matplotlib.offsetbox.AuxTransformBox(trans)
+        line = Line2D([0,0],[0,size], **linekw)
+        size_bar.add_artist(line)
+        txt = matplotlib.offsetbox.TextArea(label, minimumdescent=False,textprops={'fontsize':matplotlib.rcParams['ytick.labelsize']})
+        self.vpac = matplotlib.offsetbox.HPacker(children=[size_bar,txt],  
+                                 align="center", pad=ppad, sep=sep) 
+        matplotlib.offsetbox.AnchoredOffsetbox.__init__(self, loc, pad=pad, 
+                 borderpad=borderpad, child=self.vpac, prop=prop, frameon=frameon,
+                 **kwargs)
+
+        
+### Format for figures
+#######################################################################
+########################## Creating nice figure
+def fancy_fig_4_stimulation(sim,title=[],check_fit=False,add_scale=[],fwidealt=0,popt=[],tad=[]):
+    
+    t_v=[float(i_n) for i_n in range(0,92)]
+    lt=int(len(sim.a_Results.t))
+    if fwidealt==0:
+        f1 = plt.figure(facecolor="1",figsize=(fig_wide,fig_height*0.7))
+    else:
+        f1 = plt.figure(facecolor="1",figsize=(fwidealt,fig_height*0.7))
+    ax1 = plt.subplot2grid((5,1), (0, 0), colspan=2,rowspan=1)
+    axp = plt.subplot2grid((5,1), (1, 0), colspan=2,rowspan=4,sharex=ax1)
+    ##### Plot Input
+    ### Python 3 doesn't work with this...
+    ax1.plot(sim.a_Results.t[::20]/1000.0, sim.c_neuron.noisy_current[::20],color="black")
+    ax1.set_ylim([min(sim.c_neuron.noisy_current)-0.01,max(sim.c_neuron.noisy_current)+0.45])
+    ax1.get_xaxis().set_visible(False)
+    ax1.get_yaxis().set_visible(False)
+    #ax2.set_ylabel("Input",color="white",labelpad=10)
+    ax1.set_ylabel(r'I [uA]',labelpad=10)
+    ax1.set_title('input', loc='left')
+    locatory1 = MaxNLocator(nbins=2) # with 3 bins you will have 4 ticks
+    ax1.yaxis.set_major_locator(locatory1)
+    ax1.spines['bottom'].set_color('white')
+    ax1.spines['top'].set_color('white')
+    ax1.spines['left'].set_color('white')
+    ax1.spines['right'].set_color('white')
+    if len(add_scale)<1:
+        add_scale=[0.2,'0.2nA']
+        
+    box1 =ax1.get_position()
+    xxpos=1.12
+    yypos=0.7
+    ob = AnchoredHScaleBarY(ax=ax1,size=add_scale[0], label=add_scale[1],bbox_to_anchor=(xxpos,yypos), bbox_transform=ax1.transAxes, loc="upper right", frameon=False,
+                            pad=0.6,sep=4, linekw=dict(color="k", linewidth=0.8))
+
+    ax1.add_artist(ob)
+    for t in ax1.xaxis.get_ticklines(): t.set_color('white')
+    #### Plot complete trace
+    axp.plot(sim.a_Results.t[0:lt]/1000.0,sim.a_Results.V[0:lt],color=[0,0,0])
+#     ### show fit
+#     sim_exp=sim
+#     v_t=sim.a_Results.t
+#     sp_ampl,sp_V_max,sp_V_min=extract_spike_ampl(sim_exp.a_Results.V, sim_exp.fr[2][0:-1])
+#     tad, popt,tpeak,perr=doub_expon_fit_2(v_t[sim_exp.fr[2][0:-1]],sp_V_max,std_error=True)
+#     mqe=mean_squared_error(sp_V_max[len(sp_V_max)-len(doub_expon_fun_2(tad,*popt)):], doub_expon_fun_2(tad,*popt))
+#     # f1,axp=fancy_fig_4_stimulation(sim,title='40Hz stim '+s_Cell+' msqe='+str(mqe),check_fit=True)
+#     c=0
+#     while c < 4 and mqe>1.5:
+#         print('Trying improvement!')
+#         poptx=popt
+#         # poptx[0]=popt[0]/10
+#         tad, popt,tpeak,perr=doub_expon_fit_2(v_t[sim_exp.fr[2][0:-1]],sp_V_max,p=poptx,std_error=True)
+#         mqe=mean_squared_error(sp_V_max[len(sp_V_max)-len(doub_expon_fun_2(tad,*popt)):], doub_expon_fun_2(tad,*popt))
+#         c+=1
+    
+    if len(popt)!=0:
+        sim=sim_exp
+        sp_ampl,sp_V_max,sp_V_min=extract_spike_ampl(sim_exp.a_Results.V, sim_exp.fr[2][0:-1])
+        tad, popt,tpeak,perr=doub_expon_fit_2(v_t[sim_exp.fr[2][0:-1]],sp_V_max,std_error=True)
+        mqe=mean_squared_error(sp_V_max[len(sp_V_max)-len(doub_expon_fun_2(tad,*popt)):], doub_expon_fun_2(tad,*popt))
+        # f1,axp=fancy_fig_4_stimulation(sim,title='40Hz stim '+s_Cell+' msqe='+str(mqe),check_fit=True)
+        c=0
+        while c < 4 and mqe>1.5:
+            print('Trying improvement!')
+            poptx=popt
+            # poptx[0]=popt[0]/10
+            tad, popt,tpeak,perr=doub_expon_fit_2(v_t[sim_exp.fr[2][0:-1]],sp_V_max,p=poptx,std_error=True)
+            mqe=mean_squared_error(sp_V_max[len(sp_V_max)-len(doub_expon_fun_2(tad,*popt)):], doub_expon_fun_2(tad,*popt))
+            c+=1
+        axp.plot(tad/1000.0,doub_expon_fun_2(tad,*popt),'--',linewidth=2,color='gray')
+
+    axp.set_ylabel(r'V (mV)',labelpad=5)
+    axp.set_xlim([-10/1000.0,max(sim.a_Results.t[0:lt])/1000.0])
+    axp.set_ylim([-100,80])
+    axp.set_xlabel('time (sec)',labelpad=5)
+    locatory11 = MaxNLocator(nbins=3) # with 3 bins you will have 4 ticks
+    axp.yaxis.set_major_locator(locatory11)
+    #f1.subplots_adjust(bottom=0.17)
+    if title ==[]:
+        pass
+    else:
+        f1.suptitle(title)
+    if check_fit:
+        return f1,axp
+    else:
+        return f1

functions/f_post_simulation_analysis.py

@@ -0,0 +1,140 @@
+import numpy as np
+from scipy.signal import find_peaks_cwt
+from scipy.signal import argrelextrema
+from scipy.optimize import curve_fit
+
+################important define!
+global precision_convergence_ss
+
+precision_convergence_ss=0.00001
+
+def tic():
+    #Homemade version of matlab tic and toc functions
+    import time
+    global startTime_for_tictoc
+    startTime_for_tictoc = time.time()
+
+def toc():
+    import time
+    if 'startTime_for_tictoc' in globals():
+        print( "Elapsed time is " + str(time.time() - startTime_for_tictoc) + " seconds.")
+    else:
+        print( "Toc: start time not set")
+
+def detect_peaks(signal, threshold=0.5, spacing=3):#https://github.com/MonsieurV/py-findpeaks/blob/master/tests/libs/tony_beltramelli_detect_peaks.py
+    limit=None
+    data=signal
+    len = data.size
+    x = np.zeros(len+2*spacing)
+    x[:spacing] = data[0]-1.e-6
+    x[-spacing:] = data[-1]-1.e-6
+    x[spacing:spacing+len] = data
+    peak_candidate = np.zeros(len)
+    peak_candidate[:] = True
+    for s in range(spacing):
+        start = spacing - s - 1
+        h_b = x[start : start + len]  # before
+        start = spacing
+        h_c = x[start : start + len]  # central
+        start = spacing + s + 1
+        h_a = x[start : start + len]  # after
+        peak_candidate = np.logical_and(peak_candidate, np.logical_and(h_c > h_b, h_c > h_a))
+
+    peak_indexes = np.argwhere(peak_candidate)
+    peak_indexes = peak_indexes.reshape(peak_indexes.size)
+    if limit is not None:
+        peak_indexes = ind[data[ind] > limit]
+    return peak_indexes
+
+def firing_rate(t,V):
+    fr=np.zeros(size(t))
+    t_sp=np.zeros(size(t))
+    d=[]
+    sp=[]
+    i_n=[]
+    d=detect_peaks(V)
+    if d==[] or len(d)<2:
+        fr=np.zeros(size(t))
+        t_sp=np.zeros(size(t))
+    else:
+        sp=V[d]>thresh
+        if np.count_nonzero(sp)>2:
+            i=np.transpose(sp)*d
+            i_n=i[np.nonzero(i)]
+            t_sp=t[i_n]
+            t_sp=t_sp.reshape(size(t_sp),1)
+            fr=1/(np.diff(t_sp, n=1, axis=0)/1000)
+            t_sp=t_sp[1:]
+        else:
+            fr=np.zeros(size(t))
+            t_sp=t
+
+    return t_sp,fr,i_n
+
+def extract_spike_ampl(V, i_nn,compress=[]):
+    c=0
+    sp_V_max=[]
+    sp_V_min=[]
+    sp_ampl=[]
+    if compress !=[]:
+        i_n=[int(i/compress) for i in i_nn]
+    else:
+        i_n=i_nn
+    d_i_ni=int(i_n[0]/2)
+
+    for i in i_n:
+        if c<len(i_n)-1:
+            d_i_nd=int((i_n[c+1]-i)/2)
+        else:
+            d_i_nd=int((len(V)-i)/2)
+        sp_V_max.append(max(V[i-d_i_ni:i+d_i_nd]))
+        sp_V_min.append(min(V[i-d_i_ni:i+d_i_nd]))
+        sp_ampl.append(max(V[i-d_i_ni:i+d_i_nd])-min(V[i-d_i_ni:i+d_i_nd]))
+        c+=1
+        d_i_ni=int((i_n[c-1]-i)/2)
+
+    return sp_ampl,sp_V_max,sp_V_min 
+
+
+
+def doub_expon_fun_2(t,tau1,tau2,a,b,c):
+    if b > 0 and a > 0 and c>0 and tau1>0 and tau2>0:## Constraint to have all values positive
+        return a * np.exp(-t/tau1 ) + b * np.exp(-t/tau2 )+c
+    else:
+        return 1e10
+    
+def doub_expon_fit_2(t,x,p=[],std_error=False):
+    x=np.array(x)
+    t=np.array(t)
+    perr=[]
+    im=np.argmax(x)
+    popt=[]
+    tad=[]
+    tpeak=[]
+    adapt=False
+    if x[-1]<x[im]:
+        adapt=True
+    if adapt:
+        tad=t[im:]-t[im]
+        frad=x[im:]
+        frad02=x[im:]
+        frad01=x[im:int(len(x[im:])/10)]
+        half02=(max(frad02)-min(frad02))/2
+        half01=(max(frad01)-min(frad01))/2
+        halfid01=(np.abs(frad01-(max(frad01)-half01))).argmin()
+        halfid02=(np.abs(frad02-(max(frad02)-half02))).argmin()
+        ss_data01=frad[-int(len(frad01)/5):]
+        ss_data02=frad[-int(len(frad02)/5):]
+        ss_data01=ss_data01[abs(ss_data01 - np.mean(ss_data01)) < 2 * np.std(ss_data01)]
+        ss_data02=ss_data02[abs(ss_data02 - np.mean(ss_data02)) < 2 * np.std(ss_data02)]
+        if p!=[]:
+            popt,pcov=curve_fit(doub_expon_fun_2,tad.flatten(),frad.flatten(), p0=p)
+        else:
+            popt,pcov=curve_fit(doub_expon_fun_2,tad.flatten(),frad.flatten(), p0=[tad[halfid01],tad[halfid02]*10,x[im]-mean(ss_data01[-int(len(ss_data01)/5):]),(x[im]-mean(ss_data02[-int(len(ss_data02)/5):])),mean(ss_data02[-int(len(ss_data02)/5):])])
+
+        tpeak=t[im]
+        perr=np.sqrt(np.diag(pcov))
+    if std_error:
+        return tad, popt,tpeak, perr
+    else:
+        return tad, popt,tpeak

functions/f_saving_plots.py

@@ -0,0 +1,21 @@
+import pylab as pl
+import numpy as np
+from matplotlib.backends.backend_pdf import PdfPages
+import matplotlib.pyplot as plt
+import pickle
+
+def saving_pdf_figure(fig,fig_name,md_dir,pickle_f=True):
+    ax=plt.figure()
+    ax=fig
+    ax2=fig
+    pdffig = PdfPages(fig_name+".pdf")
+    ax.savefig(fig_name+'.png',dpi=200)
+    fig.savefig(pdffig, format="pdf")
+    metadata = pdffig.infodict()
+    metadata['Title'] = md_dir['Title']
+    metadata['Author'] = md_dir['Author']
+    metadata['Subject'] = md_dir['Subject']
+    metadata['Keywords'] = md_dir['Keywords']
+    pdffig.close()
+    if pickle_f:
+        pickle.dump(ax2, file(fig_name+'.pickle', 'w'))

functions/f_stimulations_simulations.py

@@ -0,0 +1,685 @@
+from scipy.signal import butter, lfilter, freqz
+
+### Class simulation used to store simulationss..
+class Simulation():
+        def __init__(self):
+            self.c_neuron=[]#Link to neuron or actual instance (Should contain model and parameters used)
+            self.s_Description=[]#Couple of lines describing the purpose of the simulations
+            self.d_Protocol={
+                        's_Type':[],
+                        's_Stimulus_features':[],
+                        'v_Stimulus_features':[]
+                            }
+            self.d_Configuration={
+                        'n_max_step_integration':[],
+                        'v_time_integration':[],
+                        's_ODE_Solver':[],
+                        'n_compressed_storage_factor':[]
+                                }
+            self.d_dinSys=[]
+            self.d_dinSys={
+                        'b_Convergence':[],
+                        'n_tConvergence':[],
+                        'Vars_ss':[],
+                        'Stable_ss':[],
+                        'Unstable_ss':[],
+                        's_Classification':[],
+                        'm_Jacobian4EqPoints':[],
+                        's_Stability4EqPoints':[],
+                        'm_Vect_Field':[],
+                        'n_resol_Vect_Field':[],
+                        'n_tempresol_Vect_Field':[],
+                        'm_Vect_Field_aprox_trajectory':[],
+                        'sv_state_vars_vectfield':[],
+                        'm_Vect_Field_firing_mode':[],
+                        'n_aprox_trajectory_converg_resol':[],
+                        'm_fr_overview':[]
+                            }
+            self.a_Results=[]#Results of simulation
+            self.m_Results_dinSys=[]#Results of simulation
+            self.fr=[]#Results of simulation
+            self.Adaptation=[]#If adaptation present.. Some parameters
+            self.Adaptation={
+                        'popt':[],
+                        't_peak':[],
+                        't_adapt':[],
+                            }
+
+
+class Results():
+    def __init__(self,s_Results,v_Results,compress=[],fr=[]):#s_Results contains the strings naming results, and v_Results the actual reults to be stored
+        a,b=shape(v_Results)
+        c=len(s_Results)
+        if compress==[]:
+            count=0
+            if b==c:
+                while count<len(s_Results):
+                    setattr(self, s_Results[count], v_Results[:,count])
+                    count+=1
+            if a==c:
+                while count<len(s_Results):
+                    setattr(self, s_Results[count], v_Results[count,:])
+                    count+=1
+            if a==b:
+                print("Warning, in Results storage.. columns are taken as the different results")
+                while count<len(s_Results):
+                    setattr(self, s_Results[count], v_Results[:,count])
+                    count+=1
+            if c!=a and c!=b:
+                print("Sorry! couldn't create array of results because s_Results (String of features) doesn't have the same number of entries as v_Results (Matrix with the features)")
+                while count<len(s_Results):
+                    setattr(self, s_Results[count], v_Results[:,count])
+                    count+=1
+        else:
+            count=0
+            i_spl=[]
+            i_splv=[]
+            if fr!=[]:
+                i_sp=fr[2]
+                if i_sp!=[]:
+                    i_spl=i_sp.tolist()
+                i_splv=i_spl
+                for i in range(-compress,compress):
+                    if i!=0 and i_sp!=[]:
+                        i_sptemp=i_sp+i
+                        i_spl=i_spl+i_sptemp.tolist()
+
+            if b==c:
+                v_normal=arange(0,a-1,compress).tolist()
+                v_fin=v_normal+i_spl
+                v_fin1=list(set(v_fin))
+                v_fin2=sort(v_fin1)
+
+                vect_i_sp=[]
+                for i in i_splv:
+                    vect_i_sp.append(v_fin2.tolist().index(i))
+
+                fr1=[]
+                fr1.append(fr[0])
+                fr1.append(fr[1])
+                fr1.append(vect_i_sp)
+                setattr(self, 'fr', fr1)
+
+                while count<len(s_Results):
+                    setattr(self, s_Results[count], v_Results[v_fin2,count])
+                    count+=1
+            if a==c:
+                v_normal=arange(0,b-1,compress).tolist()
+                v_fin=v_normal+i_spl
+                v_fin1=list(set(v_fin))
+                v_fin2=sort(v_fin1)
+
+                vect_i_sp=[]
+                for i in i_splv:
+                    vect_i_sp.append(v_fin2.tolist().index(i))
+                fr1=[]
+                fr1.append(fr[0])
+                fr1.append(fr[1])
+                fr1.append(vect_i_sp)
+                setattr(self, 'fr', fr1)
+
+                while count<len(s_Results):
+                    setattr(self, s_Results[count], v_Results[v_fin2,count])
+                    count+=1
+            if a==b:
+                print("Warning, in Results storage.. columns are taken as the different results")
+                while count<len(s_Results):
+                    setattr(self, s_Results[count], v_Results[:,count])
+                    count+=1
+            if c!=a and c!=b:
+                print("Sorry! couldn't create array of results because s_Results (String of features) doesn't have the same number of entries as v_Results (Matrix with the features)")
+                while count<len(s_Results):
+                    setattr(self, s_Results[count], v_Results[:,count])
+                    count+=1
+
+class Simulation_Light():
+        def __init__(self):
+            self.c_neuron=[]#Link to neuron or actual instance (Should contain model and parameters used)
+            self.fr={
+                        't':[],
+                        'sp_c':[],
+                        'sp_t':[],
+                        't_I_stim':[],
+                        's_I_stim':[],
+                        'Fr':[]
+                            }#Results of simulation
+            self.state_vars={
+                        'Na_i':[],
+                        'K_i':[],
+                        'K_o':[],
+                        'i_p':[]
+                            }#Results of simulation
+            self.d_dinSys={
+                        'saddle_node_Iapp':[],
+                        'saddle_node_Iapp_error':[]
+                            }#Results of simulation
+            self.stim_ft={
+                        's_Stimulus_features':[],
+                        'v_Stimulus_features':[]
+                            }
+            self.d_Configuration={
+                        'n_max_step_integration':[],
+                        'v_time_integration':[],
+                        's_ODE_Solver':[],
+                        'n_compressed_storage_factor':[]
+                                }
+
+def static_input_simulation(neuron,exp_duration,step_size,interval0=0):
+    ####Very fast define the stimulation protocol
+    I_exp1 = lambda t: 0 if t<interval0 else step_size
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+
+    #####Initializa the structure to save results#########
+    sim=Simulation()
+    sim.c_neuron=neuron
+    sim.s_Description="Constant input of magnitude "+str(step_size)+", to "+neuron.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Static_Input"
+    sim.d_Protocol['s_Stimulus_features']=["Input_Magnitude","Exp_Duration"]
+    sim.d_Protocol['v_Stimulus_features']=[step_size,exp_duration]
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if t<"+str(interval0)+" else "+str(step_size)
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = neuron.stimulate_neuron(t,neuron.current_state,I_exp1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+    ######And store them into the structure of results########
+    sim.fr=[t_sp,fr,i_n]
+    sim.a_Results=res
+    return sim
+
+def step_current_simulation(neuron,step_size,exp_duration,interval0=500,compress=[]):
+    ####Very fast define the stimulation protocol
+    I_exp1 = lambda t: 0 if t<interval0 else step_size
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+
+    #####Initializa the structure to save results#########
+    sim=Simulation()
+    sim.c_neuron=neuron
+    sim.s_Description="Step current of magnitude "+str(step_size)+", to "+neuron.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Step_Current"
+    sim.d_Protocol['s_Stimulus_features']=["Step_size","Exp_Duration","Time_Before_Step"]
+    sim.d_Protocol['v_Stimulus_features']=[step_size,exp_duration,interval0]
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if t<"+str(interval0)+" else "+str(step_size)
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = neuron.stimulate_neuron(t,neuron.current_state,I_exp1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+    sim.fr=[t_sp,fr,i_n]
+    if compress !=[]:
+        res=Results(s_results, v_results,compress,sim.fr)
+        sim.d_Configuration['n_compressed_storage_factor']=compress
+        sim.fr=res.fr
+    ######And store them into the structure of results########
+    sim.a_Results=res
+
+    return sim
+
+def step_current_simulation_small_interval(neuron,step_size,exp_duration,curr0=0,interval0=500,intervalf=3500,compress=[]):
+    ####Very fast define the stimulation protocol
+    I_exp1 = lambda t: curr0 if t<interval0 or t>intervalf else step_size
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+
+    #####Initializa the structure to save results#########
+    sim=Simulation()
+    sim.c_neuron=neuron
+    sim.s_Description="Step current of magnitude "+str(step_size)+", to "+neuron.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Step_Current"
+    sim.d_Protocol['s_Stimulus_features']=["Step_size","Exp_Duration","Time_Before_Step","Time_Step_Ends","Baseline_Current"]
+    sim.d_Protocol['v_Stimulus_features']=[step_size,exp_duration,interval0,intervalf,curr0]
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: "+str(curr0)+" if t<"+str(interval0)+" or t>"+str(intervalf)+" else "+str(step_size)
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = neuron.stimulate_neuron(t,neuron.current_state,I_exp1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+    sim.fr=[t_sp,fr,i_n]
+    if compress !=[]:
+        res=Results(s_results, v_results,compress,sim.fr)
+        sim.d_Configuration['n_compressed_storage_factor']=compress
+    ######And store them into the structure of results########
+    sim.a_Results=res
+    return sim
+
+def step_current_simulation_fromNonZeroStart(neuron,ini_step_size,step_size,exp_duration,interval0=500):
+    ####Very fast define the stimulation protocol
+    I_exp1 = lambda t: ini_step_size if t<interval0 else step_size
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+
+    #####Initializa the structure to save results#########
+    sim=Simulation()
+    sim.c_neuron=neuron
+    sim.s_Description="Step current of magnitude "+str(step_size)+", to "+neuron.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Step_Current"
+    sim.d_Protocol['s_Stimulus_features']=["Step_size","Exp_Duration","Time_Before_Step"]
+    sim.d_Protocol['v_Stimulus_features']=[step_size,exp_duration,interval0]
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: "+str(ini_step_size)+" if t<"+str(interval0)+" else "+str(step_size)
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = neuron.stimulate_neuron(t,neuron.current_state,I_exp1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+    sim.fr=[t_sp,fr,i_n]
+    ######And store them into the structure of results########
+    sim.a_Results=res
+    return sim
+
+def ramp_current_simulation(neuron,step_size,exp_duration,interval0=500,interval_ramp=100,I0=[]):
+    ####Very fast define the stimulation protocol
+    if I0==[]:
+        I_exp1 = lambda t: 0 if t<interval0 else step_size*t/interval_ramp
+    else:
+        I_exp1 = lambda t: I0 if t<interval0 else I0+(step_size-I0)*(t-interval0)/interval_ramp
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+
+    #####Initializa the structure to save results#########
+    sim=Simulation()
+    sim.c_neuron=neuron
+    sim.s_Description="Step current of magnitude "+str(step_size)+", to "+neuron.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Step_Current"
+    sim.d_Protocol['s_Stimulus_features']=["Step_size","Exp_Duration","Time_Before_Step"]
+    sim.d_Protocol['v_Stimulus_features']=[step_size,exp_duration,interval0]
+    if I0==[]:
+        sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if t<"+str(interval0)+" else "+str(step_size)+"*t/"+str(interval_ramp)
+    else:
+        sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: "+str(I0)+" if t<"+str(interval0)+" else "+str(I0)+"+("+str(step_size)+"-"+str(I0)+")*(t-"+str(interval0)+")/"+str(interval_ramp)
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = neuron.stimulate_neuron(t,neuron.current_state,I_exp1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+    sim.fr=[t_sp,fr,i_n]
+    ######And store them into the structure of results########
+    sim.a_Results=res
+
+    return sim
+
+def ramp_current_simulation_discrete(neuron,step_size,exp_duration,interval0=500,interval_ramp=100,dt=20,I0=[]):
+    ####Very fast define the stimulation protocol
+    if I0==[]:
+        I_exp1 = lambda t: 0 if t<interval0 else dt*int(t/dt)*step_size/interval_ramp
+    else:
+        I_exp1 = lambda t: I0 if t<interval0 else I0+dt*int((t-interval0)/dt)*(step_size-I0)/interval_ramp*dt
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+
+    #####Initializa the structure to save results#########
+    sim=Simulation()
+    sim.c_neuron=neuron
+    sim.s_Description="Step current of magnitude "+str(step_size)+", to "+neuron.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Step_Current"
+    sim.d_Protocol['s_Stimulus_features']=["Step_size","Exp_Duration","Time_Before_Step"]
+    sim.d_Protocol['v_Stimulus_features']=[step_size,exp_duration,interval0]
+    if I0==[]:
+        sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if t<"+str(interval0)+" else "+str(dt)+"*int(t/"+str(dt)+")*"+str(step_size)+"/"+str(interval_ramp)
+    else:
+        sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: "+str(I0)+" if t<"+str(interval0)+" else "+str(I0)+"+"+str(dt)+"*int((t-"+str(interval0)+")/"+str(dt)+")*("+str(step_size)+"-"+str(I0)+")/"+str(interval_ramp)
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = neuron.stimulate_neuron(t,neuron.current_state,I_exp1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+    sim.fr=[t_sp,fr,i_n]
+    ######And store them into the structure of results########
+    sim.a_Results=res
+
+    return sim
+
+
+def sine_current_simulation(neuron,sine_base,sine_amp,sine_freq,exp_duration,interval0=500,compress=[]):
+    ####Very fast define the stimulation protocol
+    freq=sine_freq*0.001
+    I_exp1 = lambda t: 0 if t<interval0 else sine_base+sine_amp/2+sine_amp/2*np.sin(2 * np.pi * freq * t + (np.pi/2))
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+
+    #####Initializa the structure to save results#########
+    sim=Simulation()
+    sim.c_neuron=neuron
+    sim.s_Description="Sine current of base "+str(sine_base)+", amplitude "+str(sine_amp)+" and frequency"+" to "+neuron.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Step_Current"
+    sim.d_Protocol['s_Stimulus_features']=["sine_base","sine_amp","sine_freq","Exp_Duration","Time_Before_Step"]
+    sim.d_Protocol['v_Stimulus_features']=[sine_base,sine_amp,sine_freq,exp_duration,interval0]
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if t<"+str(interval0)+" else "+str(sine_base)+"+"+str(sine_amp/2)+"+"+str(sine_amp/2)+"*np.sin(2 * np.pi * "+str(freq)+"* t + (np.pi/2))"
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = neuron.stimulate_neuron(t,neuron.current_state,I_exp1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+    sim.fr=[t_sp,fr,i_n]
+    if compress !=[]:
+        res=Results(s_results, v_results,compress,sim.fr)
+        sim.d_Configuration['n_compressed_storage_factor']=compress
+        sim.fr=res.fr
+    ######And store them into the structure of results########
+    sim.a_Results=res
+    return sim
+
+def sine_current_simulation_small_interval(neuron,sine_base,sine_amp,sine_freq,exp_duration,interval0=500,intervalf=3500,compress=[]):
+    ####Very fast define the stimulation protocol
+    freq=sine_freq*0.001
+    I_exp1 = lambda t: 0 if t<interval0 or t>intervalf else sine_base+sine_amp/2+sine_amp/2*np.sin(2 * np.pi * freq * t + (np.pi/2))
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+
+    #####Initializa the structure to save results#########
+    sim=Simulation()
+    sim.c_neuron=neuron
+    sim.s_Description="small interval stimulation with Sine current of base "+str(sine_base)+", amplitude "+str(sine_amp)+" and frequency"+" to "+neuron.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Sine_Current"
+    sim.d_Protocol['s_Stimulus_features']=["sine_base","sine_amp","sine_freq","Exp_Duration","Time_Before_Step","intervalf"]
+    sim.d_Protocol['v_Stimulus_features']=[sine_base,sine_amp,sine_freq,exp_duration,interval0,intervalf]
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if t<"+str(interval0)+" or t>"+str(intervalf)+" else "+str(sine_base)+"+"+str(sine_amp/2)+"+"+str(sine_amp/2)+"*np.sin(2 * np.pi * "+str(freq)+"* t + (np.pi/2))"
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = neuron.stimulate_neuron(t,neuron.current_state,I_exp1)
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+    sim.fr=[t_sp,fr,i_n]
+    if compress !=[]:
+        res=Results(s_results, v_results,compress,sim.fr)
+        sim.d_Configuration['n_compressed_storage_factor']=compress
+        sim.fr=res.fr
+    ######And store them into the structure of results########
+    sim.a_Results=res
+    return sim
+
+def butter_lowpass(cutoff, fs, order=5):
+    nyq = 0.5 * fs
+    normal_cutoff = cutoff / nyq
+    b, a = butter(order, normal_cutoff, btype='low', analog=False)
+    return b, a
+
+def butter_lowpass_filter(data, cutoff, fs, order=5):
+    b, a = butter_lowpass(cutoff, fs, order=order)
+    y = lfilter(b, a, data)
+    return y
+
+def fftnoise(f):
+    f = np.array(f, dtype='complex')
+    Np = (len(f) - 1) // 2
+    phases = np.random.rand(Np) * 2 * np.pi
+    phases = np.cos(phases) + 1j * np.sin(phases)
+    f[1:Np+1] *= phases
+    f[-1:-1-Np:-1] = np.conj(f[1:Np+1])
+    return np.fft.ifft(f).real
+
+def band_limited_noise(min_freq, max_freq, samples=1024, samplerate=1):
+    freqs = np.abs(np.fft.fftfreq(samples, 1/samplerate))
+    f = np.zeros(samples)
+    idx = np.where(np.logical_and(freqs>=min_freq, freqs<=max_freq))[0]
+    f[idx] = 1
+    return fftnoise(f)
+
+
+def noisy_current_simulation_var(neuron_in,noise_mean,noise_amp,noise_maxfreq,exp_duration,interval0=500,compress=[],light_version=0):
+    ####Very fast define the stimulation protocol
+    sim=Simulation()
+    sim.c_neuron=copy(neuron_in)
+    time_signal=exp_duration
+    samples=int(np.ceil(time_signal*(1/resol)))
+    noisy_signal_raw=band_limited_noise(0.0001, noise_maxfreq,samples,1.0/resol*1000)
+
+    #noisy_signal=noise_mean+noise_amp/2.0/max(noisy_signal_raw)*noisy_signal_raw## This is wrong... the correct way is:
+    noisy_signal=noise_mean+np.sqrt(noise_amp/np.var(noisy_signal_raw))*noisy_signal_raw## to get a signal with var=noise_amp
+    sim.c_neuron.noisy_current=noisy_signal
+    try:
+        I_exp1 = lambda t: 0 if (t<=interval0 or t>time_signal) else sim.c_neuron.noisy_current[int((t-interval0)/resol)]
+    except:
+        print("Not working at t="+str(t)+" where int((t-interval0)/resol)="+str(int((t-interval0)/resol))+" for len(sim.c_neuron.noisy_current)="+str(len(sim.c_neuron.noisy_current)))
+        raise
+    #####Initializa the structure to save results#########
+
+    sim.s_Description="Noisy current of mean "+str(noise_mean)+", amplitude "+str(noise_amp)+" and max frequency "+str(noise_maxfreq)+"Hz to "+neuron_in.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Noisy_Current"
+    sim.d_Protocol['s_Stimulus_features']=["noise_mean","noise_amp","noise_maxfreq","Exp_Duration","Time_Before_Step"]
+    sim.d_Protocol['v_Stimulus_features']=[noise_mean,noise_amp,noise_maxfreq,exp_duration,interval0]
+    # I_exp1 = lambda t: 0 if t<interval0 else noisy_signal[int((t-interval0)/resol)]
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if (t<="+str(interval0)+" or t>"+str(time_signal)+") else sim.c_neuron.noisy_current[int((t-"+str(interval0)+")/"+str(resol)+")]"
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=sim.c_neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = sim.c_neuron.stimulate_neuron(t,sim.c_neuron.current_state,I_exp1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+
+    sim.fr=[t_sp,fr,i_n]
+    if compress !=[]:
+        res=Results(s_results, v_results,compress,sim.fr)
+        sim.d_Configuration['n_compressed_storage_factor']=compress
+        sim.fr=res.fr
+    ######And store them into the structure of results########
+    if light_version==0:
+        sim.a_Results=res
+    if light_version==1:
+        sim.a_Results=[]
+    return sim
+
+def noisy_current_simulation_var_plusDepPulses(neuron_in,noise_mean,noise_amp,noise_maxfreq,exp_duration,t_pulse=[700,1000],pulse_length=20.0,pulse_ampl=2.0,interval0=500,compress=[],light_version=0):
+    ####Very fast define the stimulation protocol
+    sim=Simulation()
+    sim.c_neuron=copy(neuron_in)
+    time_signal=exp_duration
+    samples=int(np.ceil(time_signal*(1/resol)))
+    noisy_signal_raw=band_limited_noise(0.0001, noise_maxfreq,samples,1.0/resol*1000)
+
+    noisy_signal=noise_mean+noise_amp/2.0/max(noisy_signal_raw)*noisy_signal_raw
+    for tt_pulse in t_pulse:
+        for i_tlp in range(0,int(pulse_length/resol)):
+            noisy_signal[int(tt_pulse/resol+i_tlp)]=noisy_signal[int(tt_pulse/resol+i_tlp)]+pulse_ampl
+
+    sim.c_neuron.noisy_current=noisy_signal
+    try:
+        I_exp1 = lambda t: 0 if (t<=interval0 or t>time_signal) else sim.c_neuron.noisy_current[int((t-interval0)/resol)]
+    except:
+        print("Not working at t="+str(t)+" where int((t-interval0)/resol)="+str(int((t-interval0)/resol))+" for len(sim.c_neuron.noisy_current)="+str(len(sim.c_neuron.noisy_current)))
+        raise
+    #####Initializa the structure to save results#########
+
+    sim.s_Description="Noisy current of mean "+str(noise_mean)+", amplitude "+str(noise_amp)+" and max frequency "+str(noise_maxfreq)+"Hz to "+neuron_in.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Noisy_Current"
+    sim.d_Protocol['s_Stimulus_features']=["noise_mean","noise_amp","noise_maxfreq","Exp_Duration","Time_Before_Step"]
+    sim.d_Protocol['v_Stimulus_features']=[noise_mean,noise_amp,noise_maxfreq,exp_duration,interval0]
+    # I_exp1 = lambda t: 0 if t<interval0 else noisy_signal[int((t-interval0)/resol)]
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if (t<="+str(interval0)+" or t>"+str(time_signal)+") else sim.c_neuron.noisy_current[int((t-"+str(interval0)+")/"+str(resol)+")]"
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=sim.c_neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = sim.c_neuron.stimulate_neuron(t,sim.c_neuron.current_state,I_exp1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+
+    sim.fr=[t_sp,fr,i_n]
+    if compress !=[]:
+        res=Results(s_results, v_results,compress,sim.fr)
+        sim.d_Configuration['n_compressed_storage_factor']=compress
+        sim.fr=res.fr
+    ######And store them into the structure of results########
+    if light_version==0:
+        sim.a_Results=res
+    if light_version==1:
+        sim.a_Results=[]
+    return sim
+
+def noisy_current_simulation_var_fKo(neuron_in,noise_mean,noise_amp,noise_maxfreq,exp_duration,interval0=500,Ko_freq=1.0,Ko_amp=3.0,Ko_mean=3.0,compress=[],light_version=0):
+    ####Very fast define the stimulation protocol
+    sim=Simulation()
+    sim.c_neuron=copy(neuron_in)
+    time_signal=exp_duration
+    samples=int(np.ceil(time_signal*(1/resol)))
+    noisy_signal_raw=band_limited_noise(0.0001, noise_maxfreq,samples,1.0/resol*1000)
+    noisy_signal=noise_mean+noise_amp/2.0/max(noisy_signal_raw)*noisy_signal_raw
+    sim.c_neuron.noisy_current=noisy_signal
+    try:
+        I_exp1 = lambda t: 0 if (t<=interval0 or t>time_signal) else sim.c_neuron.noisy_current[int((t-interval0)/resol)]
+    except:
+        print("Not working at t="+str(t)+" where int((t-interval0)/resol)="+str(int((t-interval0)/resol))+" for len(sim.c_neuron.noisy_current)="+str(len(sim.c_neuron.noisy_current)))
+        raise
+    ## Slow Ko wave
+    try:
+        f_Ko = lambda t: Ko_mean if (t<=interval0 or t>time_signal) else Ko_mean+Ko_amp*np.sin(2*np.pi*Ko_freq*(t-interval0)/1000.0-np.pi/2)
+    except:
+        print("Not working at t="+str(t)+" where int((t-interval0)/resol)="+str(int((t-interval0)/resol))+" for len(sim.c_neuron.noisy_current)="+str(len(sim.c_neuron.noisy_current)))
+        raise
+    #####Initializa the structure to save results#########
+    sim.s_Description="Noisy current of mean "+str(noise_mean)+", amplitude "+str(noise_amp)+" and max frequency "+str(noise_maxfreq)+"Hz to "+neuron_in.s_model_tag+" parameters can be found with sim.self.c_neuron.p"
+    sim.d_Protocol['s_Type']="Noisy_Current"
+    sim.d_Protocol['s_Stimulus_features']=["noise_mean","noise_amp","noise_maxfreq","Exp_Duration","Time_Before_Step"]
+    sim.d_Protocol['v_Stimulus_features']=[noise_mean,noise_amp,noise_maxfreq,exp_duration,interval0]
+    # I_exp1 = lambda t: 0 if t<interval0 else noisy_signal[int((t-interval0)/resol)]
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+    sim.c_neuron.Ko_wave=[f_Ko(ti) for ti in t]
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if (t<="+str(interval0)+" or t>"+str(time_signal)+") else sim.c_neuron.noisy_current[int((t-"+str(interval0)+")/"+str(resol)+")]"
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[t[0],t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=sim.c_neuron.ode_solver
+    ######Run simulation################
+    s_results, v_results = sim.c_neuron.stimulate_neuron(t,sim.c_neuron.current_state,I_exp1,i_p_f=[],K_o_f=f_Ko)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    try:
+        t_sp,fr,i_n=firing_rate(res.t,res.V)
+    except:
+        try:
+            t_sp,fr,i_n=firing_rate(res.t,res.v)
+        except:
+            pass
+
+    sim.fr=[t_sp,fr,i_n]
+    if compress !=[]:
+        res=Results(s_results, v_results,compress,sim.fr)
+        sim.d_Configuration['n_compressed_storage_factor']=compress
+        sim.fr=res.fr
+    ######And store them into the structure of results########
+    if light_version==0:
+        sim.a_Results=res
+    if light_version==1:
+        sim.a_Results=[]
+    return sim
+
+def organizing_experimentalData(neuron,v_t,v_V,interval0=0,resol=resol,compress=[],light_version=0,s_Description=[]):
+    I_exp1 = lambda t: 0 if t<=interval0 else neuron.noisy_current[int((t-interval0)/resol)-1]
+    #####Initializa the structure to save results#########
+    sim=Simulation()
+    sim.c_neuron=copy(neuron)
+    sim.c_neuron.s_state_vars=['V']
+    if s_Description==[]:
+        sim.s_Description="Input current from experiment to real cell"
+    else:
+        sim.s_Description=s_Description
+    sim.d_Protocol['s_Type']="Stimulation from experiment"
+    sim.d_Protocol['s_Stimulus_features']=["Exp_Duration","Time_Before_Step"]
+    exp_duration=v_t[-1]
+    sim.d_Protocol['v_Stimulus_features']=[exp_duration,interval0]
+
+
+    # I_exp1 = lambda t: 0 if t<interval0 else noisy_signal[int((t-interval0)/resol)]
+    t=np.linspace(0, exp_duration, exp_duration/resol)
+    sim.d_Protocol['s_Executable_stimulus']="I_exp1 = lambda t: 0 if t<="+str(interval0)+" else sim.c_neuron.noisy_current[int((t-"+str(interval0)+")/"+str(resol)+")-1]"
+    sim.d_Configuration['n_max_step_integration']=resol
+    sim.d_Configuration['v_time_integration']=[v_t[0],v_t[-1]]
+    sim.d_Configuration['s_ODE_Solver']=neuron.ode_solver
+    ######Run simulation################
+    s_results=['V','t']
+    v_results=np.append(v_V[...,None],v_t[...,None],1)
+    ######Organize Simulation results into a structure####
+    res=Results(s_results, v_results)
+    t_sp,fr,i_n=firing_rate(res.t,res.V)
+    sim.fr=[t_sp,fr,i_n]
+    if compress !=[]:
+        res=Results(s_results, v_results,compress,sim.fr)
+        sim.d_Configuration['n_compressed_storage_factor']=compress
+        sim.fr=res.fr
+    ######And store them into the structure of results########
+    if light_version==0:
+        sim.a_Results=res
+    if light_version==1:
+        sim.a_Results=[]
+    return sim
+

functions/m_generic_neuron_from_json.py

@@ -0,0 +1,279 @@
+import pylab, os, json, sympy, scipy, sys
+from sympy import S, symbols, lambdify
+from scipy import integrate
+from scipy.integrate import odeint
+import numpy as np
+from copy import copy
+
+def load_mod(modfile,bp):
+
+    ### Loads expressions from the json file
+    nakdic = json.load(open(modfile))
+    ### Auxiliary functions..
+    fundic = dict([(j,k.split(":")[0]) for j,k in nakdic["functions"].items()])
+    ### Definitions..
+    defdic = dict([(j,k.split(":")[0]) for j,k in nakdic["definitions"].items()])
+    ### Parameters
+    pardic = dict([(j,k) for j,k in nakdic["parameters"].items()])
+    ### Initial conditions
+    icdic0 =  [[j,k] for j,k in nakdic["init_states"].items()]
+    icdic = [(i,str(S(str(S(str(S(str(S(j).subs(defdic))).subs(fundic))).subs(fundic))).subs(pardic))) for i,j in icdic0]
+    ### parameters to be changed
+    bifpar = [(k,pardic.pop(k)) for k in bp]
+    ### Current expressions
+    if 'currents' in nakdic.keys():
+        currentlist1 =  [[j,k] for j,k in nakdic['currents'].items()]
+        currentlist = [(i,str(S(str(S(str(S(str(S(j).subs(defdic))).subs(fundic))).subs(fundic))).subs(pardic))) for i,j in currentlist1]
+
+    if 'currents' not in nakdic.keys():
+        currentlist={}
+
+    ### SS Current expressions
+    if 'ss_currents' in nakdic.keys():
+        ss_currents1=[[j,k] for j,k in nakdic['ss_currents'].items()]
+        ss_currents = [(i,str(S(str(S(str(S(str(S(j).subs(defdic))).subs(fundic))).subs(fundic))).subs(pardic))) for i,j in ss_currents1]
+
+    if 'ss_currents' not in nakdic.keys():
+        ss_currents={}
+
+    ### Membrane potential expressions
+    if 'resting_membrane_pot' in nakdic.keys():
+        membpotlist1 =  [[j,k] for j,k in nakdic['resting_membrane_pot'].items()]
+        membpotlist = [(i,str(S(str(S(str(S(str(S(j).subs(defdic))).subs(fundic))).subs(fundic))).subs(pardic))) for i,j in membpotlist1]
+    if 'resting_membrane_pot' not in nakdic.keys():
+        membpotlist={}
+
+    ### ODEs
+    statevar_list = list(nakdic["init_states"].keys())
+    time_derivative_list= ["d{}/dt".format(k) for k in statevar_list]
+    sdelist=[]
+    for ode in nakdic["ode"]:
+      odestr= "({})-({})".format(*ode.split("="))
+      odestr= odestr.replace(" ","")
+      time_derivative= [k for k in time_derivative_list if k in odestr][0] # bad style
+      state_variable= [k for k in statevar_list if k in time_derivative][0]
+      ode_rhs= sympy.S(sympy.solve(odestr,time_derivative)[0]) # [0] is bad style
+      sdelisti = [state_variable,str(S(str(S(str(S(str(S(ode_rhs).subs(defdic))).subs(fundic))).subs(fundic))).subs(pardic))]
+      sdelist+=[sdelisti]
+    ### Model Description
+    if 'source' in nakdic.keys():
+        s_description=nakdic['source']
+    if 'source' not in nakdic.keys():
+        s_description=''
+        if "comments" in nakdic.keys():
+            if "source" in nakdic['comments'].keys():
+                s_description=nakdic['comments']['source']
+
+    return sdelist, currentlist, membpotlist, s_description, icdic, pardic, ss_currents
+
+class generic_neuron_from_json():
+
+        def __init__(self, s_name_model,dir_file=[],strIapp='I_app'):
+            ### Creates a class generic_neuron_from_json() which contains all the expressions that describe the dynamics of [s_name_model]
+            ## s_name_model= file name efample: 'MTM_W_PNAs_Temp_snapshot_constleak.json'
+            ## dir_file= path of the file containning cfg/
+            # For practival reasons I_app is always a parameter.. (to be used for stimulationg with time dependent function)
+            bifpar={}
+            if bifpar=={}:
+                I_app=-5
+                bifpar = {
+                  strIapp : [str(I_app)+"* uA/cm2"]
+                  }
+            ## Definning location of file
+            import os
+            if dir_file==[]:
+                p = {"modfile"  : "cfg/" + s_name_model, #This the file of the model
+                     "workdir"  : os.getcwd(),
+                     "bifpar"   : bifpar
+                     }
+            else:
+                p = {"modfile"  : "cfg/" + s_name_model, #This the file of the model
+                     "workdir"  : os.getenv("HOME") + dir_file,
+                     "bifpar"   : bifpar
+                     }
+
+            ##### getting expressions for currents, membrane potentials and odes
+            sde, currents, mem_pot, s_description, d_inicond, d_pars, ss_currents = load_mod(p["modfile"],p['bifpar'].keys())
+            ##### Defining units as 1 to replace later (erase them..)
+            baseunits2 = [('mV', 1), ('ms', 1),('second', 1), ('cm2', 1), ('cm3', 1), ('uF', 1), ('psiemens', 1), ('um2', 1), ('msiemens', 1), ('cm', 1), ('kelvin', 1), ('mM', 1), ('mol', 1), ('uA', 1), ('mjoule', 1), ('coulomb',1), ('ufarad',1)]
+            ##### Replacing units for ones.. in the expressions..
+            varrhs = [(i,sympy.S(j).subs(baseunits2))
+                            for i,j in sde]
+            currwounits = [(i,sympy.S(j).subs(baseunits2))
+                            for i,j in currents]
+            sscurrwounits = [(i,sympy.S(j).subs(baseunits2))
+                            for i,j in ss_currents]
+            mempwounits = [(i,sympy.S(j).subs(baseunits2))
+                            for i,j in mem_pot]
+
+            inicondwounits = [(i,sympy.S(j).subs(baseunits2))
+                            for i,j in d_inicond]
+            #### Removing units from pars
+            parslist = dict([(i,float(str(S(j).subs(baseunits2)))) for i,j in d_pars.items()])
+            #### Separating definitions frpm expressions..
+            if mempwounits==[]:
+                s_mempot=[]
+                mempotexp=[]
+            if mempwounits!=[]:
+                s_mempot,mempotexp = zip(*mempwounits)
+            if currwounits==[]:
+                s_curr=[]
+                currexp=[]
+            if currwounits!=[]:
+                s_curr,currexp = zip(*currwounits)
+            if sscurrwounits==[]:
+                s_sscurr=[]
+                sscurrexp=[]
+            if sscurrwounits!=[]:
+                s_sscurr,sscurrexp = zip(*sscurrwounits)
+
+            s_svars,svarsexp = zip(*varrhs)
+
+            self.s_pars=[j for j,k in bifpar.items()]
+            self.svarsexp=svarsexp
+            self.currexp=currexp
+            self.sscurrexp=sscurrexp
+            self.mempotexp=mempotexp
+            self.baseunits=baseunits2
+            self.s_model_tag=s_description
+            self.s_model_reference_loc=p["workdir"]
+            self.s_model_reference=p["modfile"]
+            self.p = parslist
+            self.n_state_vars=len(s_svars)
+            self.s_state_vars=[j for j in s_svars]
+            self.s_currents=[j for j in s_curr]
+            self.s_sscurrents=[j for j in s_sscurr]
+            self.s_concentrations=[]
+            self.s_resting_membrane_potentials=[j for j in s_mempot]
+            self.current_state=[float(dict(inicondwounits)[j]) for j in self.s_state_vars]
+            self.ode_solver="odeint"
+            self.noisy_current=[]
+            for i_s in self.s_state_vars:
+                if '_i' in i_s:
+                    self.s_concentrations.append(i_s)
+                if '_o' in i_s:
+                    self.s_concentrations.append(i_s)
+
+        def changing_pars(self,bifpar,pars4auto=False,strIapp='I_app'):
+            ### Function that changes parameters (Only changes pars in the instance that calls it)
+            ## bifpar= dictionary of parameters to be changed and their values
+            if strIapp in bifpar:
+                pass
+            else:
+                bifpar[strIapp]=[str(-5)+"* uA/cm2"]
+
+            p = {"modfile"  : self.s_model_reference, #"cfg/wangBuzsaki_brian.json",      #This is our model
+                 "workdir"  : self.s_model_reference_loc,
+                 "bifpar"   : bifpar
+                 }
+            ##### Recalculating expressions for currents, membrane potentials and odes
+            sde, currents, mem_pot,s_description, d_inicond, d_pars, ss_currents = load_mod(p["modfile"],p['bifpar'].keys())
+            baseunits2=self.baseunits
+            ### Recalculating expressions
+            bifpar2=dict([(j,k[0]) for j,k in bifpar.items()])
+            if pars4auto:
+                for i_p in bifpar.keys():
+                    bifpar2.pop(i_p, None)
+
+            if strIapp in bifpar2:
+                bifpar2.pop(strIapp, None)
+
+            varrhs = [(i,sympy.S(str(sympy.S(j).subs(bifpar2))).subs(baseunits2))
+                            for i,j in sde]
+            currwounits = [(i,sympy.S(str(sympy.S(j).subs(bifpar2))).subs(baseunits2))
+                            for i,j in currents]
+            mempwounits = [(i,sympy.S(str(sympy.S(j).subs(bifpar2))).subs(baseunits2))
+                            for i,j in mem_pot]
+            sscurrwounits=[(i,sympy.S(str(sympy.S(j).subs(bifpar2))).subs(baseunits2))
+                            for i,j in ss_currents]
+            #### Separating definitions from expressions..
+            if mempwounits==[]:
+                s_mempot=[]
+                mempotexp=[]
+            if mempwounits!=[]:
+                s_mempot,mempotexp = zip(*mempwounits)
+            if currwounits==[]:
+                s_curr=[]
+                currexp=[]
+            if currwounits!=[]:
+                s_curr,currexp = zip(*currwounits)
+            if sscurrwounits==[]:
+                s_sscurr=[]
+                sscurrexp=[]
+            if sscurrwounits!=[]:
+                s_sscurr,sscurrexp = zip(*sscurrwounits)
+
+            s_svars,svarsexp = zip(*varrhs)
+
+            #### adding I_app as par again.. so that instance it can be stimulated afterwards
+            I_app=0
+            bifpar_temp={}
+            bifpar_temp= {
+              strIapp : [str(I_app)+"* uA/cm2"]
+              }
+            #### Removing the annoying u before string
+            s_pars=[j for j,k in bifpar_temp.items()]
+            s_svars=[j for j in s_svars]
+            s_curr=[j for j in s_curr]
+            ## Adding the changed parameters to the parameters list again..
+            for s_i,s_bi in bifpar2.items():
+                d_pars[s_i]=s_bi
+            ## Removing units again
+            parslist = dict([(i,float(str(S(j).subs(baseunits2)))) for i,j in d_pars.items()])
+            self.s_state_vars=s_svars
+            self.s_pars=s_pars
+            self.s_curr=s_curr
+            self.mempotexp=mempotexp
+            self.svarsexp=svarsexp
+            self.currexp=currexp
+            self.p = parslist
+
+        def resting_membrane_potentials_fun(self):
+            sym_svars=symbols(self.s_state_vars)
+            fun_mempot=lambdify(sym_svars,self.mempotexp)
+            return fun_mempot
+
+        def neuron_currents_fun(self):
+            sym_svars=symbols(self.s_state_vars)
+            fun_currents = lambdify(sym_svars,self.currexp)
+            return fun_currents
+
+        def neuron_sscurrents_fun(self):
+            sym_svars=symbols(self.s_state_vars)
+            fun_sscurrents = lambdify(sym_svars,self.sscurrexp)
+            return fun_sscurrents
+
+        def neuron_changing_state_fun(self):
+            # import pdb; pdb.set_trace()
+            sym_spars=symbols(self.s_pars)
+            sym_svars=symbols(self.s_state_vars)
+            fun_ode =lambdify([i for i in sym_svars]+[i for i in sym_spars],self.svarsexp)
+            return fun_ode
+
+        def resting_membrane_potentials(self,v_State_Vars):
+            fun_mempot=self.resting_membrane_potentials_fun()
+            return fun_mempot(*v_State_Vars)
+
+        def neuron_currents(self,v_State_Vars):
+            fun_currents=self.neuron_currents_fun()
+            return fun_currents(*v_State_Vars)
+
+        def neuron_sscurrents(self,v):
+            fun_sscurrents=self.neuron_sscurrents_fun()
+            return fun_sscurrents(*[0,0,0,v])
+
+        def neuron_changing_state(self,v_State_Vars,t,I_app,fun_ode):
+            n_I=I_app(t)
+            return fun_ode(* v_State_Vars.tolist()+[n_I])
+
+        def stimulate_neuron(self,t,c_ini,I_stim):
+            fun_ode=self.neuron_changing_state_fun()
+            sol= odeint(self.neuron_changing_state, c_ini, t, args=(I_stim,fun_ode,))
+            s_results=[]
+            s_results=copy(self.s_state_vars)
+            s_results.append("t")
+            t_prov=t[...,None]
+            v_results=np.append(sol,t_prov,1)
+            return s_results, v_results
+
+##### To call it neuron=generic_neuron_from_json('MTM_W_PNAs_Temp_snapshot_constleak'+'.json')