Merge pull request #349 from kaeldai/example/bio_all_active_sweep

Adding example for running all-active cell simulation with full axon

bmtk/simulator/bionet/__init__.py

@@ -20,7 +20,7 @@
 # WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #
-from bmtk.simulator.bionet.pyfunction_cache import synapse_model, synaptic_weight, cell_model, add_weight_function
+from bmtk.simulator.bionet.pyfunction_cache import synapse_model, synaptic_weight, cell_model, add_weight_function, model_processing
 from bmtk.simulator.bionet.config import Config
 from bmtk.simulator.bionet.bionetwork import BioNetwork
 from bmtk.simulator.bionet.biosimulator import BioSimulator

bmtk/simulator/core/pyfunction_cache.py

@@ -274,6 +274,31 @@ def func_wrapper(*args, **kwargs):
                 return func(*args, **kwargs)
             return func_wrapper
         return decorator
+
+
+def model_processing(*wargs, **wkwargs):
+    if len(wargs) == 1 and callable(wargs[0]):
+        # for the case without decorator arguments, grab the function object in wargs and create a decorator
+        func = wargs[0]
+        py_modules.add_cell_processor(func.__name__, func)  # add function assigned to its original name
+
+        @wraps(func)
+        def func_wrapper(*args, **kwargs):
+            return func(*args, **kwargs)
+        return func_wrapper
+    else:
+        # for the case with decorator arguments
+        assert(all(k in ['name'] for k in wkwargs.keys()))
+
+        def decorator(func):
+            # store the function in py_modules but under the name given in the decorator arguments
+            py_modules.add_cell_processor(wkwargs['name'], func)
+
+            @wraps(func)
+            def func_wrapper(*args, **kwargs):
+                return func(*args, **kwargs)
+            return func_wrapper
+        return decorator
 
 
 def add_weight_function(func, name=None, overwrite=True):

examples/bio_all_active_sweep/README.md

@@ -0,0 +1,50 @@
+# bio_all_active_sweep
+
+Contains an example of how to download and "Biophysical, all-active" models from the 
+(Allen Cell Types Database)[https://celltypes.brain-map.org/data]. Also includes examples of how to 
+alter the "model_processing" function so that you can run cells with their full axon rather than with
+the replacement stubs that is created by default (warning: model parameters were optimized by replacing
+morphologicial axon with stubs)
+
+## Directory Structure
+* ephys_inputs/ - Folder containing NWB files downloaded for the Allen Cell Types Database. These files contain different experimental stimulus sweeps (eg square-waves, ramps, noise). Used during simulation to create current clamp stimulus onto the cells.
+* network/ - Folder containing SONATA network files, created by the `build_network.py` script and used when running the 
+* models/ - Folder containing model files (especially parameters json file and morphology swc) downloaded from Allen Cell Types Database.
+* output\*/ - Folders containing results from simulations, created by running the run_bionet.py scripts. By default contains simulation spike-times and soma membrane voltage potentials.
+* build_network.py - Python script to (re)build the SONATA network/ files, can also be used to download models/ and ephys_inputs/ files.
+* run_bionet.py - Python script to execute a simulation.
+* config.simulation.\*.json - SONATA config file containg all the information bmtk will require to run a full simulation - including location of network files, download Cell Types DB files, simulation parameters, etc.
+
+
+
+## (Re)Building the network
+
+To download the Cell-Type files and create the SONATA network files required to run a simulation sweep you can run:
+```bash
+    $ python build_network.py
+```
+
+If you want to try running a simulation with a different cell model just replace the `specimen_id` value in the python file. The script will also attempt to automatically download the model files and ephys_data nwb files if not already in the folder (if you've manually download these files the script will not try to download it again).
+
+You can also change `axon_type` to either **full** or **stub**. The default behavior is a stub axon (eg. It removes the axon and replaces it with a simple small column). But you can change it full in which case when the simulation runs it will simulate the cell using full morphology reconstruction.
+
+
+
+## Running the simulation
+
+To run a simulation with a given `config.simulation.\*.json` file run the following python command:
+
+```bash
+    $ python run_bionet.py config.simulation.json
+```
+
+BMTK will automatically load in the configuration file, run a full simulation, and save the results into the output\*/ folder (as specified in the configuration json).
+
+If you want to run a simulation with a different cell model, sweep number, or axon type you can do the following:
+1. Copy one of the existing `config.simulation.\*.json` files and open it up with your preferred text editor.
+2. In the "manifest" section update the following params as needed.
+   1. **SWEEP_NUM** - Change to run with a different stimulus sweep.
+   2. **AXON_TYPE** - Set to either *fullaxon* or *stubaxon* depending on how to handle the model's axon.
+   3. **MODEL_ID** - Set the the model being used.
+3. For some stimulus sweeps you may also need to adjust **run/tstop** (the run time in ms) otherwise the simulation may stop before sweep has finished.
+

examples/bio_all_active_sweep/build_network.py

@@ -0,0 +1,33 @@
+from bmtk.builder.networks import NetworkBuilder
+from download_data import download_ephys_data, download_model, compile_mechanisms
+
+
+def build_network(specimen_id=488683425, network_dir='network', axon_type='full', fetch_data=True):
+    processing_func = 'aibs_allactive' if axon_type == 'stub' else 'aibs_allactive_fullaxon'
+
+    # Get data from Allen CellTypes database.
+    if fetch_data:
+        model_dir = download_model(specimen_id=specimen_id)
+        download_ephys_data(specimen_id=specimen_id)
+    else:
+        model_dir='models/neuronal_model_48868368425'
+
+    # Builds a SONATA network consisting of a single cell
+    net = NetworkBuilder(f'bio_{axon_type}axon')
+    net.add_nodes(
+        N=1,
+        model_type='biophysical',
+        model_template='ctdb:Biophys1.hoc',
+        model_processing=processing_func,
+        morphology=f'reconstruction.swc',
+        dynamics_params=f'fit_parameters.json',
+
+        specimen_id=specimen_id
+    )
+    net.build()
+    net.save(output_dir=network_dir)
+
+
+if __name__ == '__main__':
+    build_network(axon_type='full')
+    build_network(axon_type='stub')

examples/bio_all_active_sweep/config.simulation.491766131_fullaxon.sweep35.json

@@ -0,0 +1,70 @@
+{
+    "manifest": {
+      "$BASE_DIR": ".",
+      "$SWEEP_NUM": "35",
+      "$AXON_TYPE": "fullaxon",
+      "$MODEL_ID": "491766131",
+      "$OUTPUT_DIR": "$BASE_DIR/output_${MODEL_ID}_${AXON_TYPE}_sweep${SWEEP_NUM}",
+      "$INPUT_DIR": "$BASE_DIR/ephys_inputs",
+      "$NETWORK_DIR": "$BASE_DIR/network",
+      "$MODEL_DIR": "$BASE_DIR/models/neuronal_model_${MODEL_ID}"
+    },
+
+    "run": {
+      "tstop": 3000.0,
+      "dt": 0.1,
+      "dL": 20.0,
+      "spike_threshold": -15,
+      "nsteps_block": 5000
+    },
+
+    "target_simulator":"NEURON",
+
+    "conditions": {
+      "celsius": 34.0,
+      "v_init": -80
+    },
+
+    "inputs": {
+      "current_clamp_nwb": {
+        "module": "IClamp",
+        "input_type": "allen",
+        "node_set": "all",
+        "file": "${INPUT_DIR}/488683423_ephys.nwb",
+        "sweep_id": "${SWEEP_NUM}"
+      }
+    },
+
+    "output":{
+      "log_file": "log.txt",
+      "output_dir": "$OUTPUT_DIR",
+      "spikes_file": "spikes.h5",
+      "spikes_file_csv": "spikes.csv",
+      "spikes_sort_order": "time",
+      "overwrite_output_dir": true
+    },
+
+    "reports": {
+      "membrane_potential": {
+        "cells": [0],
+        "variable_name": "v",
+        "module": "membrane_report",
+        "sections": "soma"
+      }
+    },
+
+    "components": {
+      "morphologies_dir": "$MODEL_DIR",
+      "mechanisms_dir":"$MODEL_DIR",
+      "biophysical_neuron_models_dir": "$MODEL_DIR"
+    },
+
+    "networks": {
+      "nodes": [
+        {
+          "nodes_file": "$NETWORK_DIR/bio_${AXON_TYPE}_nodes.h5",
+          "node_types_file": "$NETWORK_DIR/bio_${AXON_TYPE}_node_types.csv"
+        }
+      ]
+    }
+  }

examples/bio_all_active_sweep/config.simulation.491766131_stubaxon.sweep35.json

@@ -0,0 +1,70 @@
+{
+    "manifest": {
+      "$BASE_DIR": ".",
+      "$SWEEP_NUM": "35",
+      "$AXON_TYPE": "stubaxon",
+      "$MODEL_ID": "491766131",
+      "$OUTPUT_DIR": "$BASE_DIR/output_${MODEL_ID}_${AXON_TYPE}_sweep${SWEEP_NUM}",
+      "$INPUT_DIR": "$BASE_DIR/ephys_inputs",
+      "$NETWORK_DIR": "$BASE_DIR/network",
+      "$MODEL_DIR": "$BASE_DIR/models/neuronal_model_${MODEL_ID}"
+    },
+
+    "run": {
+      "tstop": 3000.0,
+      "dt": 0.1,
+      "dL": 20.0,
+      "spike_threshold": -15,
+      "nsteps_block": 5000
+    },
+
+    "target_simulator":"NEURON",
+
+    "conditions": {
+      "celsius": 34.0,
+      "v_init": -80
+    },
+
+    "inputs": {
+      "current_clamp_nwb": {
+        "module": "IClamp",
+        "input_type": "allen",
+        "node_set": "all",
+        "file": "${INPUT_DIR}/488683423_ephys.nwb",
+        "sweep_id": 35
+      }
+    },
+
+    "output":{
+      "log_file": "log.txt",
+      "output_dir": "$OUTPUT_DIR",
+      "spikes_file": "spikes.h5",
+      "spikes_file_csv": "spikes.csv",
+      "spikes_sort_order": "time",
+      "overwrite_output_dir": true
+    },
+
+    "reports": {
+      "membrane_potential": {
+        "cells": [0],
+        "variable_name": "v",
+        "module": "membrane_report",
+        "sections": "soma"
+      }
+    },
+
+    "components": {
+      "morphologies_dir": "$MODEL_DIR",
+      "mechanisms_dir":"$MODEL_DIR",
+      "biophysical_neuron_models_dir": "$MODEL_DIR"
+    },
+
+    "networks": {
+      "nodes": [
+        {
+          "nodes_file": "$NETWORK_DIR/bio_${AXON_TYPE}_nodes.h5",
+          "node_types_file": "$NETWORK_DIR/bio_${AXON_TYPE}_node_types.csv"
+        }
+      ]
+    }
+  }

examples/bio_all_active_sweep/download_data.py

@@ -0,0 +1,96 @@
+from pathlib import Path
+import os
+import io
+
+import xml.etree.ElementTree as ET 
+import requests
+import zipfile
+
+
+def download_model(specimen_id, model_type='Biophysical - all active', base_dir='models', overwrite=False):
+    """Uses the Allen REST API to find and download the model files for a given cell. All relevant files (parameters json and morphology swc)
+    will be saved in a separate folder, the path to which will be returned.
+
+    Notes: Using the REST API because at the moment I don't think the AllenSDK is capable of downloading parameters_fit.json files, see 
+    https://community.brain-map.org/t/cell-types-database-api/3016 for more info on how to use the API.
+    """
+    # Set a request to get fetch model data available for a given speciment_id. It will return a xml string that needs to be parsed so
+    # that we can find the correct model_id.
+    api_url = f"http://api.brain-map.org/api/v2/data/query.xml?criteria=model::NeuronalModel,rma::critera,[specimen_id$eq{specimen_id}]"
+    response = requests.get(api_url)
+    xml_root = ET.fromstring(response.content)
+    model_id = None
+    for (model_name, model_id) in zip(xml_root.iter('name'), xml_root.iter('id')):
+        if 'Biophysical - all active' in model_name.text:
+            model_id = int(model_id.text)
+            break
+
+    if model_id is None:
+        raise ValueError(f'Could not find a "{model_type}" model for cell {specimen_id}')
+
+    # Now that we have the model_id for the given cell we can download and unzip the files into the correct directory. To prevent downloading
+    # the zip everytime we'll check to see if the directory already exists.
+    model_dir = Path(f'{base_dir}/neuronal_model_{model_id}')
+    if model_dir.exists() and not overwrite:
+        print(f'> {model_dir} already exits, skipping donwloadng data')
+        return model_dir
+
+    zip_uri = f'http://api.brain-map.org/neuronal_model/download/{model_id}'
+    zip_response = requests.get(zip_uri)
+    zip_file = zipfile.ZipFile(io.BytesIO(zip_response.content))
+    zip_file.extractall(model_dir)
+    return model_dir
+
+
+def download_ephys_data(specimen_id, base_dir='ephys_inputs'):
+    """Download nwb file containing sweeps."""
+    api_url = f'http://api.brain-map.org/api/v2/data/query.xml?criteria=model::Specimen,rma::criteria,[id$eq{specimen_id}],rma::include,ephys_result(well_known_files(well_known_file_type[name$eqNWBDownload]))'
+    response = requests.get(api_url)
+    # print(response.content)
+    download_uri = None
+    xml_root = ET.fromstring(response.content)
+    for dl_link in xml_root.iter('download-link'):
+        download_uri = dl_link.text
+        break
+
+    ephys_id = None
+    for erid in xml_root.iter('ephys-result-id'):
+        ephys_id = erid.text
+        break
+
+    nwb_path = Path(f'{base_dir}/{ephys_id}_ephys.nwb')
+    if nwb_path.exists():
+        print(f'> {nwb_path} already exits, skipping donwload.')
+        return
+
+    if not Path(base_dir).exists():
+        os.makedirs(base_dir)
+
+    url_req = f'https://celltypes.brain-map.org/{download_uri}'
+    nwb_req = requests.get(url_req, stream=True)
+    with open(nwb_path, 'wb') as fh:
+        for chunk in nwb_req.iter_content():
+            fh.write(chunk)
+
+
+def compile_mechanisms(model_dir, modfiles_dir='modfiles', overwrite=False):
+    from subprocess import call
+    import platform
+
+    print(model_dir)
+    if not Path(f'{model_dir}/{modfiles_dir}').is_dir():
+        print(f'> Could not find directory {model_dir}/{modfiles_dir}, skipping compiling.')
+
+    if Path(platform.machine()).is_dir():
+        print(f'> {Path(platform.machine())} already existing, skipping compiling')
+
+    cwd = os.getcwd()
+    os.chdir(model_dir)
+    call(['nrnivmodl', 'modfiles'])
+    os.chdir(cwd)
+
+
+if __name__ == '__main__':
+    model_dir = download_model(specimen_id=488683425)
+    download_ephys_data(specimen_id=488683425)
+    # compile_mechanisms(model_dir)

examples/bio_all_active_sweep/models/neuronal_model_491766131/README

@@ -0,0 +1 @@
+Please visit http://alleninstitute.github.io/AllenSDK/ for instructions on running this model.