CA1 pyr cell: Inhibitory modulation of spatial selectivity+phase precession (Grienberger et al 2017)

 Download zip file 
Help downloading and running models
Accession:225080
Spatially uniform synaptic inhibition enhances spatial selectivity and temporal coding in CA1 place cells by suppressing broad out-of-field excitation.
Reference:
1 . Grienberger C, Milstein AD, Bittner KC, Romani S, Magee JC (2017) Inhibitory suppression of heterogeneously tuned excitation enhances spatial coding in CA1 place cells. Nat Neurosci [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell; Realistic Network;
Brain Region(s)/Organism: Hippocampus;
Cell Type(s): Hippocampus CA1 pyramidal cell;
Channel(s):
Gap Junctions:
Receptor(s): NMDA;
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; Python;
Model Concept(s): Active Dendrites; Detailed Neuronal Models; Place cell/field; Synaptic Integration; Short-term Synaptic Plasticity; Spatial Navigation; Feature selectivity;
Implementer(s): Milstein, Aaron D. [aaronmil at stanford.edu];
Search NeuronDB for information about:  Hippocampus CA1 pyramidal cell; NMDA;
/
GrienbergerEtAl2017
morphologies
readme.txt
ampa_kin.mod *
exp2EPSC.mod
exp2EPSG.mod
exp2EPSG_NMDA.mod
gaba_a_kin.mod *
h.mod
kad.mod *
kap.mod *
kdr.mod *
km2.mod
nas.mod
nax.mod
nmda_kin2.mod
nmda_kin3.mod
nmda_kin5.mod *
pr.mod *
vecevent.mod *
batch_EPSP_attenuation.sh
batch_place_cell_r_inp.sh
batch_place_cell_record_i_syn.sh
batch_place_cell_single_compartment.sh
batch_place_cell_subtr_inh.sh
batch_place_cell_subtr_inh_shifted.sh
batch_place_cell_subtr_inh_vclamp.sh
batch_process_i_syn_files.sh
batch_rinp.sh
batch_spine_attenuation_ratio.sh
build_expected_EPSP_reference.sh
build_expected_EPSP_reference_controller.py
build_expected_EPSP_reference_engine.py
consolidate_i_syn_files.py
consolidate_tracked_spine_data.py
fit_parameter_exponential_distribution.py
function_lib.py
optimize_AMPA_KIN.py
optimize_dendritic_excitability_020416.py
optimize_GABA_A_KIN.py
optimize_NMDA_KIN2.py
parallel_branch_cooperativity.sh
parallel_branch_cooperativity_no_nmda.sh
parallel_clustered_branch_cooperativity_nmda_controller_110315.py
parallel_clustered_branch_cooperativity_nmda_engine_110315.py
parallel_EPSP_attenuation_controller.py
parallel_EPSP_attenuation_engine.py
parallel_EPSP_i_attenuation_controller.py
parallel_EPSP_i_attenuation_engine.py
parallel_expected_EPSP_controller.py
parallel_expected_EPSP_engine.py
parallel_optimize_branch_cooperativity.sh
parallel_optimize_branch_cooperativity_nmda_kin3_controller.py
parallel_optimize_branch_cooperativity_nmda_kin3_engine.py
parallel_optimize_EPSP_amp_controller.py
parallel_optimize_EPSP_amp_engine.py
parallel_optimize_pr.sh
parallel_optimize_pr_controller_020116.py
parallel_optimize_pr_engine_020116.py
parallel_rinp_controller.py
parallel_rinp_engine.py
parallel_spine_attenuation_ratio_controller.py
parallel_spine_attenuation_ratio_engine.py
plot_channel_distributions.py
plot_NMDAR_kinetics.py
plot_results.py
plot_spine_traces.py
plot_synaptic_conductance_facilitation.py
process_i_syn_files.py
record_bAP_attenuation.py
simulate_place_cell_no_precession.py
simulate_place_cell_single_compartment.py
simulate_place_cell_single_compartment_no_nmda.py
simulate_place_cell_subtr_inh.py
simulate_place_cell_subtr_inh_add_noise.py
simulate_place_cell_subtr_inh_add_noise_no_na.py
simulate_place_cell_subtr_inh_no_na.py
simulate_place_cell_subtr_inh_no_nmda_no_na.py
simulate_place_cell_subtr_inh_r_inp.py
simulate_place_cell_subtr_inh_rec_i_syn.py
simulate_place_cell_subtr_inh_shifted.py
simulate_place_cell_subtr_inh_silent.py
simulate_place_cell_subtr_inh_vclamp.py
specify_cells.py
                            
__author__ = 'milsteina'
from specify_cells import *
from plot_results import *
import scipy.optimize as optimize
import os
"""

Aims for spike initiation at initial segment by increasing nax density and decreasing activation V1/2 relative to soma,
axon_hill, and axon compartments. Extend linear kap gradient into basals and obliques, aim for 60% spike attenuation
at bifurcation of trunk and tuft.

Hierarchical optimization:
I) optimize g_pas for target rinp at soma, trunk bifurcation, and tuft bifurcation [without h].
II) optimize ghbar_h for target rinp at soma, trunk bifurcation, and tuft bifurcation, while also optimizing for v_rest
offset between soma and tuft, and EPSP shape changes between proximal and distal synapses measured at the soma.
III) optimize gbar_nax/nas/sh/sha, gkabar_kap/d, gkdrbar for target na spike threshold, AHP amp, and vm stability

Goal is to acheive R_inp = 75, and lower f-I gain.
"""

morph_filename = 'EB2-late-bifurcation.swc'

mech_filename = '043016 Type A - km2_NMDA_KIN5_Pr'


def check_bounds(x, param_name):
    """
    For optimize_polish, based on simplex algorithm, check that the current set of parameters are within the bounds.
    :param x: array
    :param param_name: str
    :return: bool
    """
    for i in range(len(x)):
        if ((xmin[param_name][i] is not None and x[i] < xmin[param_name][i]) or
            (xmax[param_name][i] is not None and x[i] > xmax[param_name][i])):
            return False
    return True


def update_pas(x):
    """

    x0 = [2.28e-05, 1.58e-06, 58.4]
    :param x: array (soma.g_pas, trunk.g_pas slope, trunk.g_pas tau)
    """
    soma_e_pas = cell.mech_dict['soma']['pas']['e']['value']
    distance = cell.get_distance_to_node(cell.tree.root, distal_trunk, 1.)
    slope = (current_tuft_e_pas - soma_e_pas) / (np.exp(distance / x[2]) - 1.)
    cell.modify_mech_param('soma', 'pas', 'g', x[0])
    cell.modify_mech_param('trunk', 'pas', 'g', origin='soma', slope=x[1], tau=x[2])
    cell.modify_mech_param('trunk', 'pas', 'e', origin='soma', slope=slope, tau=x[2])
    for sec_type in ['apical', 'tuft']:
        cell.modify_mech_param(sec_type, 'pas', 'e', origin='trunk')
    for sec_type in ['basal', 'axon_hill', 'axon', 'ais', 'trunk', 'apical', 'tuft', 'spine_neck', 'spine_head']:
        cell.reinitialize_subset_mechanisms(sec_type, 'pas')


def update_h(x):
    """

    :param x: array [soma.ghbar_h, trunk.ghbar_h slope, trunk.ghbar_h tau, trunk.ghbar_h xhalf]
    """
    cell.modify_mech_param('soma', 'h', 'ghbar', x[0])
    cell.modify_mech_param('trunk', 'h', 'ghbar', origin='soma', slope=x[1], tau=x[2], xhalf=x[3])
    cell.modify_mech_param('basal', 'h', 'ghbar', origin='soma')
    for sec_type in ['apical', 'tuft']:
        cell.modify_mech_param(sec_type, 'h', 'ghbar', origin='trunk')


def update_v_rest(x):
    """

    :param x: array [soma.e_pas, tuft.e_pas]
    """
    distance = cell.get_distance_to_node(cell.tree.root, distal_trunk, 1.)
    tau = cell.mech_dict['trunk']['pas']['g']['tau']
    slope = (x[1] - x[0]) / (np.exp(distance / tau) - 1.)
    cell.modify_mech_param('soma', 'pas', 'e', x[0])
    cell.modify_mech_param('trunk', 'pas', 'e', origin='soma', slope=slope, tau=tau)
    for sec_type in ['apical', 'tuft']:
        cell.modify_mech_param(sec_type, 'pas', 'e', origin='trunk')


def update_sh_na(x):
    """

    :param x: array [ais.sh_nas]
    """
    cell.modify_mech_param('soma', 'nas', 'sh', x[0])
    for sec_type in ['basal', 'trunk', 'apical', 'tuft']:
        cell.modify_mech_param(sec_type, 'nas', 'sh', origin='soma')
    cell.modify_mech_param('axon_hill', 'nax', 'sh', x[0])
    for sec_type in ['ais', 'axon']:
        cell.modify_mech_param(sec_type, 'nax', 'sh', origin='axon_hill')


def update_ais_delay(x):
    """

    :param x: array [ais.sha_nas, ais.gbar_nax factor]
    """
    cell.modify_mech_param('ais', 'nax', 'sha', x[0])
    axon_gbar_nax = cell.axon[2].sec.gbar_nax
    cell.modify_mech_param('ais', 'nax', 'gbar', axon_gbar_nax * x[1])


def update_na_ka_stability(x):
    """

    :param x: array [soma.gkabar, soma.gkdrbar, axon.gkabar_kap factor, axon.gbar_nax factor]
    """
    global current_ka_trunk_factor
    global original_soma_kap
    cell.modify_mech_param('soma', 'kdr', 'gkdrbar', x[1])
    cell.modify_mech_param('soma', 'kap', 'gkabar', x[0])
    new_ka_trunk_factor = original_soma_kap * current_ka_trunk_factor / x[0]
    slope = (new_ka_trunk_factor - 1.) * x[0] / 300.
    cell.modify_mech_param('soma', 'nas', 'gbar', soma_na_gbar)
    for sec_type in ['basal', 'trunk', 'apical', 'tuft']:
        cell.reinitialize_subset_mechanisms(sec_type, 'nas')
        cell.modify_mech_param(sec_type, 'kap', 'gkabar', origin='soma', min_loc=75., value=0.)
        cell.modify_mech_param(sec_type, 'kap', 'gkabar', origin='soma', max_loc=75., slope=slope, replace=False)
        cell.modify_mech_param(sec_type, 'kad', 'gkabar', origin='soma', max_loc=75., value=0.)
        cell.modify_mech_param(sec_type, 'kad', 'gkabar', origin='soma', min_loc=75., max_loc=300., slope=slope,
                               value=x[0]+slope*75., replace=False)
        cell.modify_mech_param(sec_type, 'kad', 'gkabar', origin='soma', min_loc=300., value=x[0]+slope*300.,
                               replace=False)
        cell.modify_mech_param(sec_type, 'kdr', 'gkdrbar', origin='soma')
    cell.set_terminal_branch_na_gradient()
    cell.modify_mech_param('axon_hill', 'nax', 'gbar', soma_na_gbar)
    cell.modify_mech_param('ais', 'nax', 'gbar', soma_na_gbar * x[3] * current_ais_nax_factor)
    cell.modify_mech_param('axon', 'nax', 'gbar', soma_na_gbar * x[3])
    cell.modify_mech_param('axon_hill', 'kap', 'gkabar', origin='soma')
    cell.modify_mech_param('axon_hill', 'kdr', 'gkdrbar', origin='soma')
    for sec_type in ['ais', 'axon']:
        cell.modify_mech_param(sec_type, 'kdr', 'gkdrbar', origin='soma')
        cell.modify_mech_param(sec_type, 'kap', 'gkabar', x[0] * x[2])


def update_na_ka_dend(x):
    """

    :param x: array [soma.sh_nas, trunk.ka factor]
    """
    cell.modify_mech_param('soma', 'nas', 'sh', x[0])
    soma_gkabar = cell.mech_dict['soma']['kap']['gkabar']['value']
    slope = (x[1] - 1.) * soma_gkabar / 300.
    for sec_type in ['basal', 'trunk', 'apical', 'tuft']:
        cell.modify_mech_param(sec_type, 'kap', 'gkabar', origin='soma', min_loc=75., value=0.)
        cell.modify_mech_param(sec_type, 'kap', 'gkabar', origin='soma', max_loc=75., slope=slope, replace=False)
        cell.modify_mech_param(sec_type, 'kad', 'gkabar', origin='soma', max_loc=75., value=0.)
        cell.modify_mech_param(sec_type, 'kad', 'gkabar', origin='soma', min_loc=75., max_loc=300., slope=slope,
                               value=soma_gkabar+slope*75., replace=False)
        cell.modify_mech_param(sec_type, 'kad', 'gkabar', origin='soma', min_loc=300., value=soma_gkabar+slope*300.,
                               replace=False)
        cell.modify_mech_param(sec_type, 'nas', 'sh', origin='soma')
    cell.modify_mech_param('axon_hill', 'nax', 'sh', x[0])
    for sec_type in ['ais', 'axon']:
        cell.modify_mech_param(sec_type, 'nax', 'sh', origin='axon_hill')


def offset_vm(sec_type, vm_target=None):
    """

    :param sec_type: str
    """
    if vm_target is None:
        vm_target = v_init
    sim.modify_stim(0, amp=0.)
    if sec_type == 'soma':
        sim.modify_stim(1, node=cell.tree.root, loc=0., amp=0.)
        rec = sim.rec_list[0]['vec']
    elif sec_type == 'trunk':
        sim.modify_stim(1, node=trunk, loc=1., amp=0.)
        rec = sim.rec_list[1]['vec']
    elif sec_type == 'tuft':
        sim.modify_stim(1, node=tuft, loc=1., amp=0.)
        rec = sim.rec_list[2]['vec']
    offset = True
    sim.tstop = equilibrate
    sim.run(vm_target)
    t = np.arange(0., equilibrate, dt)
    vm = np.interp(t, sim.tvec, rec)
    v_rest = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
    initial_v_rest = v_rest
    sim.modify_stim(1, amp=i_holding[sec_type])
    sim.run(vm_target)
    vm = np.interp(t, sim.tvec, rec)
    v_rest = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
    if v_rest < vm_target - 0.5:
        i_holding[sec_type] += 0.01
        while offset:
            if sim.verbose:
                print 'increasing i_holding to %.3f (%s)' % (i_holding[sec_type], sec_type)
            sim.modify_stim(1, amp=i_holding[sec_type])
            sim.run(vm_target)
            vm = np.interp(t, sim.tvec, rec)
            v_rest = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
            if v_rest < vm_target - 0.5:
                i_holding[sec_type] += 0.01
            else:
                offset = False
    elif v_rest > vm_target + 0.5:
        i_holding[sec_type] -= 0.01
        while offset:
            if sim.verbose:
                print 'decreasing i_holding to %.3f (%s)' % (i_holding[sec_type], sec_type)
            sim.modify_stim(1, amp=i_holding[sec_type])
            sim.run(vm_target)
            vm = np.interp(t, sim.tvec, rec)
            v_rest = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
            if v_rest > vm_target + 0.5:
                i_holding[sec_type] -= 0.01
            else:
                offset = False
    sim.tstop = duration
    return initial_v_rest


def get_spike_shape(vm):
    """

    :param t: array
    :param vm: array
    :return: tuple of float: (v_peak, th_v, ADP, AHP)
    """
    vm = vm[int((equilibrate+0.4)/dt):]
    dvdt = np.gradient(vm, [dt])
    th_x = np.where(dvdt > th_dvdt)[0]
    if th_x.any():
        th_x = th_x[0] - int(1.6/dt)
    else:
        th_x = np.where(vm > -30.)[0][0] - int(2./dt)
    th_v = vm[th_x]
    v_before = np.mean(vm[th_x-int(0.1/dt):th_x])
    v_peak = np.max(vm[th_x:th_x+int(5./dt)])
    x_peak = np.where(vm[th_x:th_x+int(5./dt)] == v_peak)[0][0]
    v_AHP = np.min(vm[th_x+x_peak:th_x+int(12./dt)])
    x_AHP = np.where(vm[th_x+x_peak:th_x+int(12./dt)] == v_AHP)[0][0]
    AHP = v_before - v_AHP
    # if spike waveform includes an ADP before an AHP, return the value of the ADP in order to increase error function
    rising_x = np.where(dvdt[th_x+x_peak:th_x+x_peak+x_AHP] > 0.)[0]
    if rising_x.any():
        v_ADP = np.max(vm[th_x+x_peak+rising_x[0]:th_x+x_peak+x_AHP])
        ADP = v_ADP - v_AHP
    else:
        ADP = 0.
    return v_peak, th_v, ADP, AHP


def get_EPSP_dur(t, vm, equilibrate, duration):
    """

    :param t:
    :param vm:
    :param equilibrate:
    :param duration:
    :return:
    """
    dt = 0.001
    interp_t = np.arange(0., duration, dt)
    interp_vm = np.interp(interp_t, t, vm)
    left = int((equilibrate-3.)/dt)
    right = int((equilibrate-1.)/dt)
    baseline = np.mean(interp_vm[left:right])
    start = int(equilibrate/dt)
    end = int(duration/dt)
    interp_t = interp_t[start:end]
    interp_vm = interp_vm[start:end] - baseline
    amp = np.max(interp_vm)
    t_peak = np.where(interp_vm == amp)[0][0]
    interp_vm /= amp
    interp_t -= interp_t[0]
    rise_50 = np.where(interp_vm[:t_peak] >= 0.5)[0][0]
    decay_50 = np.where(interp_vm[t_peak:] <= 0.5)[0][0]
    decay_dur = interp_t[decay_50] + interp_t[t_peak] - interp_t[rise_50]
    return decay_dur


def pas_error(x, plot=0):
    """
    :param x: array (soma.g_pas, trunk.g_pas slope, trunk.g_pas tau)
    :param plot: int
    :return: float
    """
    if np.any(x < 0.):
        return 1e9
    start_time = time.time()
    sim.tstop = duration
    result = {}
    amp = -0.15
    cell.zero_h()
    cell.zero_na()
    update_pas(x)
    offset_vm('soma')
    sim.modify_stim(0, node=cell.tree.root, loc=0., amp=amp, dur=stim_dur)
    sim.run(v_init)
    result['soma'] = get_Rinp(np.array(sim.tvec), np.array(sim.rec_list[0]['vec']), equilibrate, duration, amp)[2]
    offset_vm('trunk')
    sim.modify_stim(0, node=trunk, loc=1., amp=amp)
    sim.run(v_init)
    result['trunk'] = get_Rinp(np.array(sim.tvec), np.array(sim.rec_list[1]['vec']), equilibrate, duration, amp)[2]
    offset_vm('tuft')
    sim.modify_stim(0, node=tuft, loc=1., amp=amp)
    sim.run(v_init)
    result['tuft'] = get_Rinp(np.array(sim.tvec), np.array(sim.rec_list[2]['vec']), equilibrate, duration, amp)[2]
    sim.modify_stim(1, amp=0.)
    Err = 0.
    for target in result:
        Err += ((target_val['pas'][target] - result[target])/target_range['pas'][target])**2.
    print('Simulation took %.3f s' % (time.time()-start_time))
    print 'Process %i: [soma g_pas, trunk slope, trunk tau]: [%.2E, %.2E, %.1f], soma R_inp: %.1f, ' \
          'trunk R_inp: %.1f, tuft R_inp: %.1f' % (os.getpid(), x[0], x[1], x[2], result['soma'],
                                                   result['trunk'], result['tuft'])
    print 'Process %i: Error: %.4E' % (os.getpid(), Err)
    if plot:
        sim.plot()
    global lowest_Err
    global best_x
    if Err < lowest_Err:
        lowest_Err = float(Err)
        best_x = list(x)
    return Err


def h_error(x, plot=0):
    """

    :param x: array (soma.ghbar_h, trunk.ghbar_h slope, trunk.ghbar_h tau, trunk.ghbar_h xhalf)
    :param plot: int
    :return: float
    """
    if np.any(x < 1.e-10) or x[2] < 25.:  # or x[2] + x[3] > 300.:
        return 1e9
    start_time = time.time()
    sim.tstop = duration
    result = {}
    amp = -0.15
    cell.zero_na()
    update_h(x)
    soma_vm = offset_vm('soma')
    result['v_rest'] = soma_vm
    sim.modify_stim(0, node=cell.tree.root, loc=0., amp=amp, dur=stim_dur)
    sim.run(v_init)
    result['soma'] = get_Rinp(np.array(sim.tvec), np.array(sim.rec_list[0]['vec']), equilibrate, duration, amp)[2]
    offset_vm('trunk')
    sim.modify_stim(0, node=trunk, loc=1., amp=amp)
    sim.run(v_init)
    result['trunk'] = get_Rinp(np.array(sim.tvec), np.array(sim.rec_list[1]['vec']), equilibrate, duration, amp)[2]
    sim.modify_stim(0, amp=0.)
    syn = trunk.spines[-1].synapses[0]
    syn.source.play(h.Vector([equilibrate]))
    sim.run(v_init)
    trunk_EPSP_dur = get_EPSP_dur(np.array(sim.tvec), np.array(sim.rec_list[0]['vec']), equilibrate, duration)
    syn.source.play(h.Vector())
    tuft_vm = offset_vm('tuft')
    sim.modify_stim(0, node=tuft, loc=1., amp=amp)
    sim.run(v_init)
    tuft_rinp = get_Rinp(np.array(sim.tvec), np.array(sim.rec_list[2]['vec']), equilibrate, duration, amp)[2]
    result['v_rest_offset'] = tuft_vm - soma_vm
    if plot:
        sim.plot()
    sim.modify_stim(0, amp=0.)
    syn = tuft.spines[-1].synapses[0]
    syn.source.play(h.Vector([equilibrate]))
    sim.run(v_init)
    tuft_EPSP_dur = get_EPSP_dur(np.array(sim.tvec), np.array(sim.rec_list[0]['vec']), equilibrate, duration)
    syn.source.play(h.Vector())
    result['delta_EPSP_dur'] = max(0., trunk_EPSP_dur - tuft_EPSP_dur)
    sim.modify_stim(1, amp=0.)
    Err = 0.
    for target in result:
        Err += ((target_val['h'][target] - result[target])/target_range['h'][target])**2.
    print 'Simulation took %i s' % (time.time()-start_time)
    print 'Process %i: [soma.ghbar, trunk ghbar slope, trunk ghbar tau, trunk ghbar xhalf]: [%.2E, %.2E, %.1f, %.1f],' \
          ' soma R_inp: %.1f, trunk R_inp: %.1f, tuft R_inp: %.1f, v_rest: %.1f, v_rest offset: %.1f, ' \
          'EPSP_dur delta: %.2f' % (os.getpid(), x[0], x[1], x[2], x[3], result['soma'], result['trunk'], tuft_rinp,
                                    soma_vm, result['v_rest_offset'], result['delta_EPSP_dur'])
    print 'Process %i: Error: %.4E' % (os.getpid(), Err)
    global lowest_Err
    global best_x
    if Err < lowest_Err:
        lowest_Err = float(Err)
        best_x = list(x)
    return Err


def v_rest_check(plot=0):
    """

    :param plot: int
    :return: float
    """
    start_time = time.time()
    sim.tstop = equilibrate
    result = {}
    sim.modify_stim(0, amp=0.)
    sim.modify_stim(1, amp=0.)
    sim.run(v_init)
    t = np.arange(0., equilibrate, dt)
    vm = np.interp(t, sim.tvec, sim.rec_list[0]['vec'])
    result['soma'] = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
    vm = np.interp(t, sim.tvec, sim.rec_list[1]['vec'])
    trunk_rest = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
    vm = np.interp(t, sim.tvec, sim.rec_list[2]['vec'])
    tuft_rest = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
    result['tuft_offset'] = tuft_rest - result['soma']
    sim.tstop = duration
    Err = 0.
    for target in result:
        Err += ((target_val['v_rest'][target] - result[target])/target_range['v_rest'][target])**2.
    print 'Simulation took %i s' % (time.time()-start_time)
    print 'Process %i: soma v_rest: %.1f, trunk v_rest: %.1f, tuft v_rest: %.1f, tuft offset: %.1f' % (os.getpid(),
                                                        result['soma'], trunk_rest, tuft_rest, result['tuft_offset'])
    print 'Process %i: Error: %.4E' % (os.getpid(), Err)
    if plot:
        sim.plot()


def v_rest_error(x, plot=0):
    """

    :param x: array (soma.e_pas, tuft.e_pas)
    :param plot: int
    :return: float
    """
    if x[1] > x[0] or not check_bounds(x, 'v_rest'):
        print 'Aborting: Invalid parameter values.'
        return 1e9
    start_time = time.time()
    update_v_rest(x)
    sim.tstop = equilibrate
    result = {}
    sim.modify_stim(0, amp=0.)
    sim.modify_stim(1, amp=0.)
    sim.run(v_init)
    t = np.arange(0., equilibrate, dt)
    vm = np.interp(t, sim.tvec, sim.rec_list[0]['vec'])
    result['soma'] = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
    vm = np.interp(t, sim.tvec, sim.rec_list[1]['vec'])
    trunk_rest = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
    vm = np.interp(t, sim.tvec, sim.rec_list[2]['vec'])
    tuft_rest = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
    result['tuft_offset'] = tuft_rest - result['soma']
    sim.tstop = duration
    Err = 0.
    for target in result:
        Err += ((target_val['v_rest'][target] - result[target])/target_range['v_rest'][target])**2.
    print 'Simulation took %i s' % (time.time()-start_time)
    print 'Process %i: [soma.e_pas, tuft.e_pas]: [%.1f, %.1f], soma v_rest: %.1f, trunk v_rest: %.1f, tuft v_rest: ' \
          '%.1f, tuft offset: %.1f' % (os.getpid(), x[0], x[1], result['soma'], trunk_rest, tuft_rest,
                                       result['tuft_offset'])
    print 'Process %i: Error: %.4E' % (os.getpid(), Err)
    if plot:
        sim.plot()
    return Err


def na_ka_threshold_error(x, plot=0):
    """

    :param x: array [sh_nax/s]
    :param plot: int
    :return: float
    """
    if x[0] < 0.:
        print 'Process %i: Aborting - Parameters outside optimization bounds.' % (os.getpid())
        return 1e9
    start_time = time.time()
    update_sh_na(x)
    sim.modify_stim(0, node=cell.tree.root, loc=0., dur=100.)
    duration = equilibrate + 200.
    sim.tstop = duration
    t = np.arange(0., duration, dt)
    spike = False
    amp = 0.05
    while not spike:
        sim.modify_stim(0, amp=amp)
        sim.run(v_init)
        vm = np.interp(t, sim.tvec, sim.rec_list[0]['vec'])
        if amp == 0.05 and np.any(vm[:int(equilibrate/dt)] > -30.):
            print 'Process %i: Aborting - spontaneous firing' % (os.getpid())
            return 1e9
        if np.any(vm[int(equilibrate/dt):int((equilibrate+50.)/dt)] > -30.):
            spike = True
        else:
            amp += 0.05
            if sim.verbose:
                print 'increasing amp to %.3f' % amp
    peak, threshold, ADP, AHP = get_spike_shape(vm)
    result = {'th_v': threshold}
    Err = 0.
    for target in result:
        Err += ((target_val['na_ka'][target] - result[target])/target_range['na_ka'][target])**2.
    print 'Simulation took %i s' % (time.time()-start_time)
    print 'Process %i: [sh.nax\s]: [%.3f], amp: %.3f, threshold: %.2f' % (os.getpid(), x[0], amp, threshold)
    print 'Process %i: Error: %.4E' % (os.getpid(), Err)
    global lowest_Err
    global best_x
    if Err < lowest_Err:
        lowest_Err = float(Err)
        best_x = list(x)
    return Err


def na_ka_stability_error(x, plot=0):
    """

    :param x: array [soma.gkabar, soma.gkdrbar, axon.gkabar_kap factor, axon.gbar_nax factor]
    :param plot: int
    :return: float
    """
    if not check_bounds(x, 'na_ka_stability'):
        print 'Process %i: Aborting - Parameters outside optimization bounds.' % (os.getpid())
        return 1e9
    start_time = time.time()
    update_na_ka_stability(x)
    soma_vm = offset_vm('soma', v_active)
    result = {'v_rest': soma_vm}
    sim.modify_stim(0, node=cell.tree.root, loc=0., dur=100.)
    duration = equilibrate + 200.
    sim.tstop = duration
    t = np.arange(0., duration, dt)
    spike = False
    amp = 0.05
    while not spike:
        sim.modify_stim(0, amp=amp)
        sim.run(v_active)
        vm = np.interp(t, sim.tvec, sim.rec_list[0]['vec'])
        if amp == 0.05 and np.any(vm[:int(equilibrate/dt)] > -30.):
            print 'Process %i: Aborting - spontaneous firing' % (os.getpid())
            return 1e9
        if np.any(vm[int(equilibrate/dt):int((equilibrate+50.)/dt)] > -30.):
            spike = True
        else:
            amp += 0.05
            if sim.verbose:
                print 'increasing amp to %.3f' % amp
    peak, threshold, ADP, AHP = get_spike_shape(vm)
    if plot:
        sim.plot()
    result['ADP'] = ADP
    result['AHP'] = AHP
    slow_depo = 0.
    stability = 0.
    for new_amp in [0.25, 0.5, 0.75]:
        print 'increasing amp to %.3f' % (new_amp + amp)
        sim.modify_stim(0, amp=(new_amp + amp))
        sim.run(v_active)
        if plot:
            sim.plot()
        vm = np.interp(t, sim.tvec, sim.rec_list[0]['vec'])
        v_rest = np.mean(vm[int((equilibrate - 3.)/dt):int((equilibrate - 1.)/dt)])
        v_before = np.max(vm[int((equilibrate - 50.)/dt):int((equilibrate - 1.)/dt)])
        v_after = np.max(vm[-int(50./dt):-1])
        stability += abs((v_before - v_rest) + (v_after - v_rest))
        if stability > 15.:
            stability += 100.
            break
        v_min_late = np.min(vm[int((equilibrate + 80.)/dt):int((equilibrate + 99.)/dt)])
        this_slow_depo = v_min_late - threshold
        if new_amp != 0.75 and v_min_late >= -35.:
            slow_depo += this_slow_depo + (0.75 - new_amp) * 1000.
            break
        elif this_slow_depo > 0.:
            slow_depo += this_slow_depo
    result['stability'] = stability
    result['slow_depo'] = slow_depo
    Err = 0.
    for target in result:
        # don't penalize AHP or slow_depo less than target
        if not ((target == 'AHP' and result[target] < target_val['na_ka'][target]) or
                (target == 'slow_depo' and result[target] < target_val['na_ka'][target])):
            Err += ((target_val['na_ka'][target] - result[target])/target_range['na_ka'][target])**2.
    print 'Simulation took %i s' % (time.time()-start_time)
    print 'Process %i: [soma.gkabar, soma.gkdrbar, axon.gkabar_kap factor, axon.gbar_nax factor]: ' \
          '[%.4f, %.4f, %.2f, %.2f], amp: %.3f, v_rest: %.1f, threshold: %.1f, ADP: %.1f, AHP: %.1f, ' \
          'stability: %.2f, slow_depo: %.2f' % (os.getpid(), x[0], x[1], x[2], x[3], amp, soma_vm,
                                                threshold, ADP, AHP, result['stability'], result['slow_depo'])
    print 'Process %i: Error: %.4E' % (os.getpid(), Err)
    global lowest_Err
    global best_x
    if Err < lowest_Err:
        lowest_Err = float(Err)
        best_x = list(x)
    return Err


def ais_delay_error(x, plot=0):
    """
    Part of a set of error functions designed to optimize one parameter at a time by changing as few variables as
    possible. This function just changes the V1/2 of activation of the axon initial segment and the relative density of
    the nax sodium channel mechanism in the axon initial segment to accomplish initiation of the action potential in the
    initial segment.
    :param x: array [ais.sha_nas, ais.gbar_nax factor]
    :param plot: int
    :return: float
    """
    if x[0] > 0. or x[0] < -5. or x[1] < 1. or x[1] > 5.:
        print 'Process %i: Aborting - Parameters outside optimization bounds.' % (os.getpid())
        return 1e9
    start_time = time.time()
    update_ais_delay(x)
    offset_vm('soma', v_active)
    sim.modify_stim(0, node=cell.tree.root, loc=0., dur=100.)
    duration = equilibrate + 200.
    sim.tstop = duration
    t = np.arange(0., duration, dt)
    spike = False
    amp = 0.05
    while not spike:
        sim.modify_stim(0, amp=amp)
        sim.run(v_active)
        vm = np.interp(t, sim.tvec, sim.rec_list[0]['vec'])
        if amp == 0.05 and np.any(vm[:int(equilibrate/dt)] > -30.):
            print 'Process %i: Aborting - spontaneous firing' % (os.getpid())
            return 1e9
        if np.any(vm[int(equilibrate/dt):int((equilibrate+50.)/dt)] > -30.):
            spike = True
        else:
            amp += 0.05
            if sim.verbose:
                print 'increasing amp to %.3f' % amp
    peak, threshold, ADP, AHP = get_spike_shape(vm)
    result = {}
    th_x = np.where(vm[int(equilibrate/dt):] >= threshold)[0][0] + int(equilibrate/dt)
    ais_vm = np.interp(t, sim.tvec, sim.rec_list[3]['vec'])
    ais_dvdt = np.gradient(ais_vm, [dt])
    axon_vm = np.interp(t, sim.tvec, sim.rec_list[4]['vec'])
    axon_dvdt = np.gradient(axon_vm, [dt])
    left = th_x - int(2./dt)
    right = th_x + int(2./dt)
    ais_peak = np.max(ais_dvdt[left:right])
    ais_peak_t = np.where(ais_dvdt[left:right] == ais_peak)[0][0] * dt
    axon_peak = np.max(axon_dvdt[left:right])
    axon_peak_t = np.where(axon_dvdt[left:right] == axon_peak)[0][0] * dt
    if axon_peak_t > ais_peak_t + dt:
        result['ais_delay'] = 0.
    else:
        result['ais_delay'] = ais_peak_t + dt - axon_peak_t
    if plot:
        sim.plot()
    Err = 0.
    for target in result:
        Err += ((target_val['na_ka'][target] - result[target])/target_range['na_ka'][target])**2.
    print 'Simulation took %i s' % (time.time()-start_time)
    print 'Process %i: [ais.sha_nas, ais.gbar_nas factor]: [%.2f, %.2f], ais_delay: %.3E' % (os.getpid(), x[0], x[1],
                                            result['ais_delay'])
    print 'Process %i: Error: %.4E' % (os.getpid(), Err)
    global lowest_Err
    global best_x
    if Err < lowest_Err:
        lowest_Err = float(Err)
        best_x = list(x)
    return Err


def na_ka_dend_error(x, plot=0):
    """

    :param x: array [soma.sh_nas, trunk.ka factor]
    :param plot: int
    :return: float
    """
    start_time = time.time()
    update_na_ka_dend(x)
    offset_vm('soma', v_active)
    sim.modify_stim(0, node=cell.tree.root, loc=0., dur=100.)
    sim.modify_rec(1, node=distal_trunk, loc=1.)
    duration = equilibrate + 200.
    sim.tstop = duration
    t = np.arange(0., duration, dt)
    spike = False
    amp = 0.05
    while not spike:
        sim.modify_stim(0, amp=amp)
        sim.run(v_active)
        vm = np.interp(t, sim.tvec, sim.rec_list[0]['vec'])
        if amp == 0.05 and np.any(vm[:int(equilibrate/dt)] > -30.):
            print 'Process %i: Aborting - spontaneous firing' % (os.getpid())
            return 1e9
        if np.any(vm[int(equilibrate/dt):int((equilibrate+50.)/dt)] > -30.):
            spike = True
        else:
            amp += 0.05
            if sim.verbose:
                print 'increasing amp to %.3f' % amp
    peak, threshold, ADP, AHP = get_spike_shape(vm)
    result = {'th_v': threshold}
    trunk_vm = np.interp(t, sim.tvec, sim.rec_list[1]['vec'])
    th_x = np.where(vm[int(equilibrate/dt):] >= threshold)[0][0] + int(equilibrate/dt)
    trunk_peak = np.max(trunk_vm[th_x:th_x+int(10./dt)])
    trunk_pre = np.mean(trunk_vm[th_x-int(0.2/dt):th_x-int(0.1/dt)])
    result['trunk_amp'] = (trunk_peak - trunk_pre) / (peak - threshold)
    if plot:
        sim.plot()
    Err = 0.
    for target in result:
        Err += ((target_val['na_ka'][target] - result[target])/target_range['na_ka'][target])**2.
    print 'Simulation took %i s' % (time.time()-start_time)
    print 'Process %i: [sh_nax/s, trunk.ka factor]: [%.2f, %.2f], amp: %.3f, threshold: %.1f, trunk_amp: %.3f' % \
          (os.getpid(), x[0], x[1], amp, threshold, result['trunk_amp'])
    print 'Process %i: Error: %.4E' % (os.getpid(), Err)
    sim.modify_rec(1, node=trunk, loc=1.)
    global lowest_Err
    global best_x
    if Err < lowest_Err:
        lowest_Err = float(Err)
        best_x = list(x)
    return Err


def optimize_polish((param_name, x), maxfev=None):
    """

    :param param_name: str
    :param x: array
    :param xmin: array
    :param xmax: array
    :param maxfev: int
    :return: array
    """
    if maxfev is None:
        maxfev = 200
    error_functions = {}
    error_functions['pas'] = pas_error
    error_functions['h'] = h_error
    error_functions['v_rest'] = v_rest_error
    error_functions['na_ka_stability'] = na_ka_stability_error
    error_functions['na_ka_dend'] = na_ka_dend_error

    result = optimize.minimize(error_functions[param_name], x, method='Nelder-Mead', options={'ftol': 1e-3,
                                                    'xtol': 1e-3, 'disp': True, 'maxiter': maxfev})
    formatted_x = '['+', '.join(['%.2E' % xi for xi in result.x])+']'
    print 'Process: %i optimized %s for %i iterations with Error: %.2f and x: %s' % (os.getpid(), param_name,
                                                                            result.nit, result.fun, formatted_x)
    return {param_name: {'x': result.x, 'Err': result.fun}}


def optimize_explore((param_name, x, xmin, xmax), maxfev=None):
    """

    :param param_name: str
    :param x: array
    :param xmin: array
    :param xmax: array
    :param maxfev: int
    :return: array
    """
    if maxfev is None:
        maxfev = 200
    error_functions = {}
    error_functions['pas'] = pas_error
    error_functions['h'] = h_error
    error_functions['v_rest'] = v_rest_error
    error_functions['na_ka_stability'] = na_ka_stability_error
    error_functions['na_ka_dend'] = na_ka_dend_error

    take_step = Normalized_Step(x, xmin, xmax)
    minimizer_kwargs = dict(method=null_minimizer)
    result = optimize.basinhopping(error_functions[param_name], x, niter=maxfev, niter_success=maxfev/2,
                                       disp=True, interval=20, minimizer_kwargs=minimizer_kwargs, take_step=take_step)
    formatted_x = '['+', '.join(['%.2E' % xi for xi in result.x])+']'
    print 'Process: %i optimized %s for %i iterations with Error: %.2f and x: %s' % (os.getpid(), param_name,
                                                                            result.nit, result.fun, formatted_x)
    return {param_name: {'x': result.x, 'Err': result.fun}}


equilibrate = 250.  # time to steady-state
stim_dur = 500.
duration = equilibrate + stim_dur
dt = 0.01
amp = 0.3
th_dvdt = 10.
v_init = -67.
v_active = -61.
spines = True  # False
i_holding = {'soma': 0., 'trunk': 0., 'tuft': 0.}

cell = CA1_Pyr(morph_filename, mech_filename, full_spines=spines)

# look for a trunk bifurcation
trunk_bifurcation = [trunk for trunk in cell.trunk if cell.is_bifurcation(trunk, 'trunk')]
if trunk_bifurcation:
    trunk_branches = [branch for branch in trunk_bifurcation[0].children if branch.type == 'trunk']
    # get where the thickest trunk branch gives rise to the tuft
    trunk = max(trunk_branches, key=lambda node: node.sec(0.).diam)
    trunk = (node for node in cell.trunk if cell.node_in_subtree(trunk, node) and 'tuft' in (child.type
                                                                            for child in node.children)).next()
else:
    trunk_bifurcation = [node for node in cell.trunk if 'tuft' in (child.type for child in node.children)]
    trunk = trunk_bifurcation[0]
tuft = (child for child in trunk.children if child.type == 'tuft').next()
distal_trunk = trunk
trunk = trunk_bifurcation[0]
tuft = distal_trunk  # testing optimizing pas and h at tuft bifurcation instead of in thin tuft

for branch in [trunk, tuft]:
    syn = Synapse(cell, branch.spines[-1], ['AMPA_KIN'], stochastic=0)
cell.init_synaptic_mechanisms()

soma_ek = -77.
cell.modify_mech_param('soma', 'ions', 'ek', soma_ek)
"""
cell.modify_mech_param('soma', 'nas', 'ar', 1.)
for sec_type in ['basal', 'trunk', 'apical', 'tuft']:
    cell.modify_mech_param(sec_type, 'nas', 'ar', 0.8)
"""
cell.modify_mech_param('soma', 'pas', 'e', -65.)
for sec_type in ['spine_neck', 'spine_head', 'axon_hill', 'ais', 'axon', 'basal']:
    cell.modify_mech_param(sec_type, 'pas', 'e', v_init)

soma_na_gbar = 0.04
"""
cell.modify_mech_param('ais', 'nax', 'sha', -2.8)
for sec_type in ['basal', 'trunk', 'apical', 'tuft']:
    cell.modify_mech_param(sec_type, 'nas', 'sha', 5.)
"""
cell.modify_mech_param('soma', 'km2', 'gkmbar', 0.0015)
cell.modify_mech_param('axon_hill', 'km2', 'gkmbar', origin='soma')
cell.modify_mech_param('ais', 'km2', 'gkmbar', 0.0075)
cell.modify_mech_param('axon', 'km2', 'gkmbar', origin='ais')

cell.reinit_mechanisms()
cell.set_terminal_branch_na_gradient()


sim = QuickSim(duration)  # , verbose=False)
sim.append_stim(cell, cell.tree.root, loc=0., amp=0., delay=equilibrate, dur=stim_dur)
sim.append_stim(cell, cell.tree.root, loc=0., amp=0., delay=0., dur=duration)
sim.append_rec(cell, cell.tree.root, loc=0.)
sim.append_rec(cell, trunk, loc=1.)
sim.append_rec(cell, tuft, loc=1.)  # loc=0.)

axon_seg_locs = [seg.x for seg in cell.axon[2].sec]

sim.append_rec(cell, cell.axon[1], loc=1., description='ais')
sim.append_rec(cell, cell.axon[2], loc=axon_seg_locs[0], description='axon')

#the target values and acceptable ranges
target_val = {}
target_range = {}
target_val['pas'] = {'soma': 120., 'trunk': 97., 'tuft': 121.}
target_range['pas'] = {'soma': 1., 'trunk': 1., 'tuft': 1.}
target_val['h'] = {'v_rest': v_init, 'soma': 75., 'trunk': 39., 'tuft': 42., 'v_rest_offset': 0., 'delta_EPSP_dur': 0.}
target_range['h'] = {'v_rest': 0.25, 'soma': 0.5, 'trunk': 1., 'tuft': 2., 'v_rest_offset': 0.1, 'delta_EPSP_dur': 0.05}
target_val['v_rest'] = {'soma': v_init, 'tuft_offset': 0.}
target_range['v_rest'] = {'soma': 0.25, 'tuft_offset': 0.1}
target_val['na_ka'] = {'v_rest': v_init, 'th_v': -51., 'soma_peak': 40., 'trunk_amp': 0.6, 'ADP': 0., 'AHP': 4.,
                       'stability': 0., 'ais_delay': 0., 'slow_depo': 25.}
target_range['na_ka'] = {'v_rest': 0.25, 'th_v': .2, 'soma_peak': 2., 'trunk_amp': 0.01, 'ADP': 0.01, 'AHP': .2,
                         'stability': 1., 'ais_delay': 0.001, 'slow_depo': 1.}

x0 = {}
x1 = {}
xmin = {}
xmax = {}

# x0['pas'] = [1.43E-06, 1.58E-05, 121.9]
# x0['pas'] = [1.46E-06, 1.56E-05, 124.2]  # Error: 9.5227E+00
x0['pas'] = [1.52E-06, 1.63E-05, 121.9]  # Error: 2.6862E+00
xmin['pas'] = [1.0E-6, 1.0E-06, 25.]
xmax['pas'] = [2.0E-5, 2.0E-05, 200.]

# x0['h'] = [5.93E-08, 4.88E-03, 114.4, 459.0]
# x0['h'] = [2.27E-09, 2.14E-03, 265.3, 484.0]  # Error: 4.0003E+02
# x0['h'] = [2.35E-09, 1.76E-03, 381.8, 280.2]  # Error: 4.0260E+02
x0['h'] = [2.43E-09, 1.85E-03, 383.1, 260.7]  # Error: 3.5202E+02
xmin['h'] = [1.0E-9, 1.e-4, 25., 50.]
xmax['h'] = [1.0E-7, 1.e-2, 200., 500.]

# [soma.gkabar, soma.gkdrbar, axon.gkabar_kap factor, axon.gbar_nax factor]
# x0['na_ka_stability'] = [0.03555, 0.04975, 4.51, 2.18]
# x0['na_ka_stability'] = [0.0305, 0.0478, 2.97, 2.23]
x0['na_ka_stability'] = [0.0305, 0.0478, 2.97, 2.34]  # Error: 4.2416E+02
xmin['na_ka_stability'] = [0.01, 0.01, 1., 2.]
xmax['na_ka_stability'] = [0.05, 0.05, 5., 5.]

# [soma.sh_nas, trunk.ka factor]
# x0['na_ka_dend'] = [0., 2.9]
x0['na_ka_dend'] = [1.7, 2.8]
xmin['na_ka_dend'] = [0., 1.5]
xmax['na_ka_dend'] = [4., 5.]

# [ais.sha_nas, ais.gbar_nax factor]
x0['ais_delay'] = [-3.6, 4.4]

# [soma.e_pas, tuft.e_pas]
x0['v_rest'] = [-63., -77.]
xmin['v_rest'] = [v_init, soma_ek]
xmax['v_rest'] = [-63., -63.]

x1 = dict(x0)
log = []
niter = 200

current_tuft_e_pas = x1['v_rest'][1]
update_v_rest(x1['v_rest'])
update_pas(x1['pas'])
update_h(x1['h'])

current_ais_nax_factor = x1['ais_delay'][1]
current_ka_trunk_factor = x1['na_ka_dend'][1]
original_soma_kap = x1['na_ka_stability'][0]

update_na_ka_stability(x1['na_ka_stability'])
update_na_ka_dend(x1['na_ka_dend'])
update_ais_delay(x1['ais_delay'])

best_x = []
lowest_Err = 1e9  # haven't gotten this to work with global variables for some reason


#na_ka_stability_error(x1['na_ka_stability'], 1)

"""
Err = 1e9
i = 0
x1_ais = [-2., 3.]
while Err > 10.:
    Err = ais_delay_error(x1_ais)
    if Err > 10:
        if i % 2 == 0:
            x1_ais[0] -= 0.2
        else:
            x1_ais[1] += 0.2
        i += 1
"""

"""
optimization 020416:
v_rest error:
[soma.e_pas, tuft.e_pas]: [-63.0, -77.0], soma v_rest: -68.8, trunk v_rest: -68.1, tuft v_rest: -68.0, tuft offset: 0.7
Process 81997: Error: 1.0624E+02

pas_error:
[soma g_pas, trunk slope, trunk tau]: [1.52E-06, 1.63E-05, 121.9],
soma R_inp: 120.4, trunk R_inp: 97.3, tuft R_inp: 120.0
Process 81966: Error: 1.2615E+00

h_error:
[soma.ghbar, trunk ghbar slope, trunk ghbar tau, trunk ghbar xhalf]: [2.43E-09, 1.85E-03, 383.1, 260.7],
soma R_inp: 73.4, trunk R_inp: 53.7, tuft R_inp: 74.6, v_rest: -68.8, v_rest offset: 0.8, EPSP_dur delta: 0.00
Process 81997: Error: 3.3479E+02

ais_delay_error:
[ais.sha_nas, ais.gbar_nas factor]: [-3.60, 4.40], ais_delay: 0.000E+00
Process 69150: Error: 0.0000E+00

na_ka_dend_error:
[sh_nax/s, trunk.ka factor]: [1.70, 2.80], amp: 0.150, threshold: -51.1, trunk_amp: 0.608
Process 81817: Error: 8.1919E-01

na_ka_stability_error:
[soma.gkabar, soma.gkdrbar, axon.gkabar_kap factor, axon.gbar_nax factor]:
[0.0305, 0.0478, 2.97, 2.34], amp: 0.150, v_rest: -68.8, threshold: -51.1, ADP: 0.0, AHP: 8.1, stability: 0.44,
slow_depo: 31.65
Process 82023: Error: 5.0614E+02
"""