Merge pull request #154 from BlueBrain/memodel-threshold

compute threshold in memodel example if not threshold-based protocols

bluepyemodel/emodel_pipeline/emodel_settings.py

@@ -79,6 +79,7 @@ def __init__(
         neuron_dt=None,
         cvode_minstep=0.0,
         use_params_for_seed=True,
+        current_precision=1e-2,
         max_threshold_voltage=-30,
         strict_holding_bounds=True,
         max_depth_holding_search=7,
@@ -268,6 +269,8 @@ def __init__(
             cvode_minstep (float): minimum time step allowed when using cvode.
             use_params_for_seed (bool): use a hashed version of the parameter
                 dictionary as a seed for the simulator
+            current_precision (float): size of search interval in current to stop the search
+                in SearchThresholdProtocol
             max_threshold_voltage (float): upper bound for the voltage during the
                 search for the threshold or rheobase current (see SearchThresholdProtocol).
             strict_holding_bounds (bool): if True, the minimum and maximum values for the current
@@ -305,6 +308,7 @@ def __init__(
         self.minimum_protocol_delay = minimum_protocol_delay
 
         # Settings related to the evaluator
+        self.current_precision = current_precision
         self.max_threshold_voltage = max_threshold_voltage
         self.threshold_efeature_std = threshold_efeature_std
         self.max_depth_holding_search = max_depth_holding_search

bluepyemodel/evaluation/evaluator.py

@@ -464,6 +464,7 @@ def define_holding_protocol(
 
 def define_threshold_protocol(
     efeatures,
+    current_precision=1e-2,
     max_threshold_voltage=-30,
     step_delay=500.0,
     step_duration=2000.0,
@@ -497,6 +498,7 @@ def define_threshold_protocol(
         name=thres_prot_name,
         location=location,
         target_threshold=target_current[0],
+        current_precision=current_precision,
         max_threshold_voltage=max_threshold_voltage,
         stimulus_delay=step_delay,
         stimulus_duration=step_duration,
@@ -626,6 +628,7 @@ def define_preprotocols(
     fitness_calculator_configuration,
     strict_holding_bounds=False,
     max_depth_holding_search=7,
+    current_precision=1e-2,
     max_threshold_voltage=-30,
     spikecount_timeout=50,
     max_depth_threshold_search=10,
@@ -685,6 +688,7 @@ def define_preprotocols(
     )
     threshold_prot = define_threshold_protocol(
         efeatures,
+        current_precision,
         max_threshold_voltage,
         fitness_calculator_configuration.search_threshold_step_delay,
         fitness_calculator_configuration.search_threshold_step_duration,
@@ -717,6 +721,7 @@ def define_threshold_based_optimisation_protocol(
     stochasticity=True,
     ais_recording=False,
     efel_settings=None,
+    current_precision=1e-2,
     max_threshold_voltage=-30,
     strict_holding_bounds=True,
     use_fixed_dt_recordings=False,
@@ -740,6 +745,7 @@ def define_threshold_based_optimisation_protocol(
             deterministic
         ais_recording (bool): if True all the soma recording will be at the first axonal section.
         efel_settings (dict): eFEl settings.
+        current_precision (float): size of search interval in current to stop the search
         max_threshold_voltage (float): maximum voltage at which the SearchThresholdProtocol
             will search for the rheobase.
         strict_holding_bounds (bool): to adaptively enlarge bounds if holding current is outside
@@ -789,6 +795,7 @@ def define_threshold_based_optimisation_protocol(
                             fitness_calculator_configuration=fitness_calculator_configuration,
                             strict_holding_bounds=strict_holding_bounds,
                             max_depth_holding_search=max_depth_holding_search,
+                            current_precision=current_precision,
                             max_threshold_voltage=max_threshold_voltage,
                             spikecount_timeout=spikecount_timeout,
                             max_depth_threshold_search=max_depth_threshold_search,
@@ -820,6 +827,7 @@ def define_threshold_based_optimisation_protocol(
                             fitness_calculator_configuration=fitness_calculator_configuration,
                             strict_holding_bounds=strict_holding_bounds,
                             max_depth_holding_search=max_depth_holding_search,
+                            current_precision=current_precision,
                             max_threshold_voltage=max_threshold_voltage,
                             spikecount_timeout=spikecount_timeout,
                             max_depth_threshold_search=max_depth_threshold_search,
@@ -847,6 +855,7 @@ def define_threshold_based_optimisation_protocol(
                         fitness_calculator_configuration=fitness_calculator_configuration,
                         strict_holding_bounds=strict_holding_bounds,
                         max_depth_holding_search=max_depth_holding_search,
+                        current_precision=current_precision,
                         max_threshold_voltage=max_threshold_voltage,
                         spikecount_timeout=spikecount_timeout,
                         max_depth_threshold_search=max_depth_threshold_search,
@@ -976,6 +985,7 @@ def create_evaluator(
             include_validation_protocols,
             stochasticity=stochasticity,
             efel_settings=pipeline_settings.efel_settings,
+            current_precision=pipeline_settings.current_precision,
             max_threshold_voltage=pipeline_settings.max_threshold_voltage,
             strict_holding_bounds=pipeline_settings.strict_holding_bounds,
             use_fixed_dt_recordings=use_fixed_dt_recordings,

examples/memodel/memodel.py

@@ -199,6 +199,75 @@ def plot_scores(access_point, cell_evaluator, mapper, figures_dir, seed):
     return emodel_score
 
 
+def get_default_threshold_search_protocol():
+    """Create a default protocol to use to search for the threshold with no holding current."""
+    from bluepyemodel.evaluation.evaluator import define_efeatures
+    from bluepyemodel.evaluation.evaluator import define_preprotocols
+    from bluepyemodel.evaluation.evaluator import soma_loc
+    from bluepyemodel.evaluation.fitness_calculator_configuration import FitnessCalculatorConfiguration
+    from bluepyemodel.evaluation.protocols import ProtocolRunner
+
+    rmp_prot_name = "RMPProtocol_noholding"
+    rin_prot_name = "RinProtocol_noholding"
+    efeatures_dict = [
+        {
+            "efel_feature_name": "steady_state_voltage_stimend",
+            "protocol_name": rmp_prot_name,
+            "recording_name": "soma.v",
+            "mean": 0,
+        },
+        {
+            "efel_feature_name": "ohmic_input_resistance_vb_ssse",
+            "protocol_name": rin_prot_name,
+            "recording_name": "soma.v",
+            "mean": 0,
+        },
+    ]
+    fcc = FitnessCalculatorConfiguration(
+        efeatures=efeatures_dict,
+        name_rmp_protocol=rmp_prot_name,
+        name_rin_protocol=rin_prot_name,
+    )
+    efeatures = define_efeatures(
+        fcc,
+        include_validation_protocols=False,
+        protocols={},
+        efel_settings={},
+    )
+
+    protocols = define_preprotocols(
+        efeatures=efeatures,
+        location=soma_loc,
+        fitness_calculator_configuration=fcc,
+        rmp_key="bpo_rmp_noholding",
+        hold_key="bpo_holding_current_noholding",
+        rin_key="bpo_rin_noholding",
+        thres_key="bpo_threshold_current_noholding",
+        rmp_prot_name=rmp_prot_name,
+        hold_prot_name="SearchHoldingCurrent_noholding",
+        rin_prot_name=rin_prot_name,
+        thres_prot_name="SearchThresholdCurrent_noholding",
+        recording_name="soma.v",
+        no_holding=True,
+    )
+
+    return ProtocolRunner(protocols)
+
+def get_threshold(cell_evaluator, access_point, mapper):
+    """Computes and returns threshold current."""
+    evaluator = copy.deepcopy(cell_evaluator)
+
+    protocol = get_default_threshold_search_protocol()
+    evaluator.fitness_protocols = {"main_protocol": protocol}
+    emodels = compute_responses(
+        access_point,
+        evaluator,
+        mapper,
+    )
+
+    return emodels[0].responses.get("bpo_threshold_current_noholding", None)
+
+
 def plot(access_point, seed, cell_evaluator, figures_dir, mapper):
     """Plot figures and return total fitness (sum of scores), holding and threshold currents"""
     # compute scores
@@ -383,6 +452,9 @@ def update_memodel(
         access_token=access_token,
     )
 
+    # update settings for better threshold precision
+    access_point.pipeline_settings.current_precision = 2e-3
+
     # get cell evaluator with 'new' morphology
     cell_evaluator = get_cell_evaluator(access_point, morph_name, morph_format, morph_id)
 
@@ -403,6 +475,14 @@ def update_memodel(
     emodel_score, emodel_holding, emodel_threshold = plot(access_point, seed, cell_evaluator, figures_dir, mapper)
     if not add_score:
         emodel_score = None
+
+    # if we have absolute amplitude protocols, and threshold current was not computed,
+    # then compute it
+    if emodel_threshold is None:
+        # assume holding = 0
+        emodel_holding = 0
+        emodel_threshold = get_threshold(cell_evaluator, access_point, mapper)
+
     # attention! after this step, do NOT push EModel again.
     # It has been modified and would overwrite the correct one on nexus