Generating biophysically detailed multi-compartmental models¶

from pathlib import Path
tutorial_output_dir = f"{Path.home()}/isf_tutorial_output"  # <-- Change this to your desired output directory

from biophysics_fitting.hay.default_setup import get_feasible_model_params

params = get_feasible_model_params().drop('x', axis=1)
params.index = [e.replace('CaDynamics_E2', 'CaDynamics_E2_v2') for e in params.index]
params.index = 'ephys.' + params.index

from biophysics_fitting.L5tt_parameter_setup import get_L5tt_template_v2
import biophysics_fitting.hay.default_setup as hay_setup

def scale_apical(cell_param, params):
    assert(len(params) == 1)
    cell_param.cell_modify_functions.scale_apical.scale = params['scale']
    return cell_param

def get_fixed_params(db_setup):
    """
    Configure the fixed params and return
    """
    fixed_params = db_setup['fixed_params']
    fixed_params['morphology.filename'] = db_setup['morphology'].get_file(
        'hoc')
    return fixed_params

def get_Simulator(db_setup, step=False):
    """
    Configure the Simulator object and return
    """
    fixed_params = db_setup['get_fixed_params'](db_setup)
    s = hay_setup.get_Simulator(
        I.pd.Series(fixed_params),
        step=step)
    s.setup.cell_param_generator = get_L5tt_template_v2
    s.setup.cell_param_modify_funs.append(
        ('scale_apical', scale_apical)
        )
    s.setup.params_modify_funs_after_cell_generation = []
    return s

def get_Evaluator(db_setup, step=False):
    """
    No additional configuration is needed for the Evaluator, simply return biophysics_fitting.L5tt_parameter_setup.get_Evaluator
    """
    return hay_setup.get_Evaluator(step=step)

def get_Combiner(db_setup, step=False):
    """
    No additional configuration is needed for the Combiner, simply return biophysics_fitting.L5tt_parameter_setup.get_Combiner
    """
    return hay_setup.get_Combiner(step=step)

def set_up_db_for_MOEA(db, morphology_id='89', morphology="", step=False):
    """
    Set up a DataBase for MOEA.

    Args:
        db: a DataBase object
        morphology_id: name of the morphology
        morphology: path to a .hoc morphology file
        step: whether or not to perform step current injections

    Returns:
        data_base.DataBase: a database containing:
            - fixed_params
            - get_fixed_params
            - get_Simulator
            - get_Evaluator
            - get_Combiner
            - the morphology file
    """
    from data_base.IO.LoaderDumper import pandas_to_pickle, to_cloudpickle

    db.create_sub_db(morphology_id)
    db[morphology_id].create_managed_folder('morphology')
    I.shutil.copy(
        I.os.path.join(
            morphology
        ), db[morphology_id]['morphology'].join(
            morphology.split(I.os.sep)[-1]
        ))

    db[morphology_id]['fixed_params'] = {
        'BAC.hay_measure.recSite': 294.8203371921156,   # recording site on the apical dendrite for BAC
        'BAC.stim.dist': 294.8203371921156,             # stimulus injection site on the apical dendrite for BAC
        'bAP.hay_measure.recSite1': 294.8203371921156,  # recording site 1 on the apical dendrite for bAP
        'bAP.hay_measure.recSite2': 474.8203371921156,  # recording site 2 on the apical dendrite for bAP
        'hot_zone.min_': 384.8203371921156,             # calcium zone start
        'hot_zone.max_': 584.8203371921156,             # calcium zone end
        'hot_zone.outsidescale_sections': [23,24,25,26,27,28,29,31,32,33,34,35,37,38,40,42,43,44,46,48,50,51,52,54,56,58,60],
        'morphology.filename': None
        }

    db[morphology_id]['get_fixed_params'] = get_fixed_params
    db[morphology_id].set('params', params, dumper=pandas_to_pickle)
    db[morphology_id].set('get_Simulator',
                      I.partial(get_Simulator, step=step),
                      dumper=to_cloudpickle)
    db[morphology_id].set('get_Evaluator',
                      I.partial(get_Evaluator, step=step),
                      dumper=to_cloudpickle)
    db[morphology_id].set('get_Combiner',
                      I.partial(get_Combiner, step=step),
                      dumper=to_cloudpickle)

    return db

morphology_id = "89"
morphology_path = I.os.path.join(
    example_data_dir,
    "anatomical_constraints",
    "89_C2.hoc")

if "morphology" in db[morphology_id].keys():
    # In case you are re-running this tutorial
    del db[morphology_id]["morphology"]

seedpoint = 42
if str(seedpoint) in db[morphology_id].keys():
    # In case you are re-running this tutorial
    del db[morphology_id][str(seedpoint)]

if not I.os.path.exists(I.os.path.join(db["89"].basedir, 'RW_exploration_example')):
    db["89"].create_managed_folder('RW_exploration_example')

evaluator = db["89"]["get_Evaluator"](db['89'])
simulator = db["89"]["get_Simulator"](db['89'])
biophysical_parameter_ranges = db['89']['params']
example_models = model_db['example_models']
biophysical_parameter_names = [e for e in example_models.columns if "ephys" in e or e == "scale_apical.scale"]

from biophysics_fitting.exploration_from_seedpoint.RW import RW
from biophysics_fitting.exploration_from_seedpoint.utils import evaluation_function_incremental_helper
from biophysics_fitting.hay.evaluation import hay_evaluate_bAP, hay_evaluate_BAC

evaluation_function = I.partial(
    evaluation_function_incremental_helper,  # this incremental helper stops the evaluation as soon as a stimulus protocol doesn't work.
    s=simulator,
    e=evaluator,
    stim_order=['bAP', 'BAC']
)

rw = RW(
    param_ranges = biophysical_parameter_ranges,
    df_seeds = example_models[biophysical_parameter_names],
    evaluation_function = evaluation_function,
    MAIN_DIRECTORY = db["89"]['RW_exploration_example'],
    min_step_size = 0.02,
    max_step_size = 0.02,
    checkpoint_every = 1  # This is a lot of IO. Increase this value if you are not merely doing an example.
)

from biophysics_fitting.exploration_from_seedpoint.RW_analysis import Load

outdir = db["89"]['RW_exploration_example'].join('0')
l = Load(
    client=I.get_client(),
    path=outdir,
    n_particles = 1)

def normalize(df, mn, mx):
    return (df - mn)/(mx-mn)

mn, mx = biophysical_parameter_ranges['min'], biophysical_parameter_ranges['max']
normalized_startpoint = normalize(example_models.iloc[0][biophysical_parameter_names], mn, mx)
normalized_explored_models = normalize(explored_models[biophysical_parameter_names], mn, mx)

# calc exploration relative to startpoint, in % of total allowed parameter limits
d = I.pd.concat([normalized_explored_models, I.pd.DataFrame(normalized_startpoint).T])
d -= normalized_startpoint
d[biophysical_parameter_names] *= 100
d = d.melt(var_name='Biophysical parameter', value_name='Normalized value (%)')

delayeds = [
    I.dask.delayed(simulator.run)(p, 'BAC')
    for _, p in explored_models[biophysical_parameter_names].iterrows()
    ]
f = I.get_client().compute(delayeds)

responses = [f_.result() for f_ in f]
BAC_responses = [response['BAC.hay_measure'] for response in responses]

soma_voltages = [
    e['vList'][0] for e in BAC_responses
]
dend_voltages = [
    e['vList'][1] for e in BAC_responses
]
time_points = [
    e['tVec'] for e in BAC_responses
]

def plot_traces(ax, time_points, voltages, highlight_id, color, title):
    for t, v in zip(time_points, voltages):
        ax.plot(t, v, c="silver")
        ax.plot(
            time_points[highlight_id],
            voltages[highlight_id],
            c=color,
            lw=2
            )
        ax.set_title(title)
        ax.set_xlabel("Time (ms)")
        ax.set_ylabel("Membrane voltage (mV)")
        ax.set_xlim(250, 450)

cell, _ = simulator.get_simulated_cell(
    explored_models[biophysical_parameter_names].iloc[random_model],
    "BAC")