Recreating the Model

The model was created in stages: cell models, network models, odor input, LFP output, simulations, and experiments. While these stages are described in the dissertation, the sections below describe technical details and include links to specific files, tools, and resources. The sections are organized in a chronological order, which can be followed to re-create this or a similar model.

Mitral, tufted, and granule cell morphologies cleaned, morphology metrics computed, and representative morphologies chosen

Morphology Archives

Morphologies for each cell type (mitral, tufted, granule) were obtained from NeuroMorpho.org, by performing a ‘Metadata search’ of reconstruction archives. The archives were filtered to include cells with the following parameters: mouse, adult, olfactory bulb region, control condition, and mitral, tufted, or granule cell type.

Inspection and Editing of the SWC Files

Each resulting archive was downloaded and organized into folders which can be found under [repo]/morphology-data. The 3D renderings of the .SWC files under each archive were manually examined using the neuTube software. neuTube functions were used to fix simple and obvious flaws (e.g. parts of dendritic trees shifted laterally, or all sections mislabeled as ‘basal’ dendrites). Archives with cells missing important stereotypical features of each cell type were excluded. Besides, neuTube, other neuronal morphology/SWC editing software was evaluated, but ultimately not used.

Each SWC archive was labeled for whether it was a ‘thin’ (parts of lateral dendrites missing) or ‘full’ (the whole cell was reconstructed) slice. Whether the cells included a soma reconstruction, and whether the reconstructed dendrites included width/radius information (e.g. some archives used constant radius).

The fixed SWC files were saved back into the archive folders and archive categorizations were stored in the [repo]/notebooks/morphology.ipynb Jupyter Notebook.

Computing Reconstructed Cell Population Morphology Metrics

For each reconstruction, NeuroMorpho.org displays a set of whole-cell morphology metrics (e.g. see ‘Measurements’ section of this reconstruction). For cells in this model, those metrics were computed for soma, apical dendrites, and lateral dendrites of each reconstruction.

The metrics were computed using the pyLMeasure tool (which is a Python wrapper of the L-Measure tool). Each metric was formalized as tests of the NeuronUnit framework. The test definitions have been incorporated into NeuronUnit and can be found under [neuronunit]/tests/morphology.py

The morphology metrics of each cell were computed using suites of NeuronUnit tests (see: [repo]/notebooks/morphology.ipynb) and results stored in the model’s SQLite Database in the measurement table. The database tables and records can be viewed online with tools like SQLite online or offline with an SQLite editor like SQLite Expert.

The overall means and standard deviations of each morphology metric were stored in the database in the property table. These overall means and standard deviations were used to select five stereotypical morphologies of each cell type. Reconstructions whose morphology metrics were closest to the cell type population means were selected to be used for cell models in this model. They can be found under [repo]/prev_ob_models/Birgiolas2020/SWCs.

The validation of the representativeness of the selected morphologies and comparison with previously published models can be found in the morphology-validation.ipynb notebook

Morphology Standardization to Prepare for Network Embedding

To facilitate the embedding of the cell models within reconstructed network layers, the cell morphologies were further edited as follows. The cells were reoriented so that their apical dendrites were aligned with the Z-axis. The natural curvature of lateral dendrites of mitral and tufted cells was planarized so that the lateral dendrites were essentially perpendicular to the apical dendrites. This prepared the lateral dendrites for an algorithm to confine them to the appropriate cell-type-specific laminar regions in the network model. These standardized morphologies can be found under [repo]/prev_ob_models/Birgiolas2020/morphology.

Cell electrical properties and behavior identified and characterized

Once cell morphologies were selected, they were turned into active models by inserting ion channel mechanisms. To understand which ion channels should be included, a survey of electrical behaviors of each cell type was conducted via a literature search.

Generic Test and Publication Python Classes

Electrical properties of mitral, tufted, and granule cells from main olfactory bulbs of adult mice were aggregated from literature. These properties included passive (e.g. input resistance and resting voltage) and active membrane properties (e.g. spike half-width, spike accommodation, and inhibitory rebound). The protocols used to measure and compute each property were formalized using a set of NeuronUnit tests. Each publication used slightly different protocols and algorithms to compute the properties (e.g. different slice temperatures, different delays to measure input resistance, different ways to compute action potential threshold). To fairly evaluate how a model reproduces a property, these protocol differences were formalized into the NeuronUnit tests as well. A given property reported in a publication was formalized as a NeuronUnit test that consisted of 1) a generic electrophysiology property Python class (e.g. RestingVoltageTest) and 2) publication-specific modifications of the property measurement protocol (e.g. BurtonUrban2014).

A NeuronUnit test that computes a property using protocol variations from a specific publication is declared by combining the generic test class with a publication class using Python’s multiple inheritance (e.g. class InputResistanceUrbanBurton2014Test(UrbanBurton2014, InputResistanceTest) ). Due to the large number of property and publication combinations, these combined classes were created at run time using a helper method.

The measurements from each publication are recorded in the SQLite database, measurement table. The NeuronUnit class to use for each measurement is defined in the property table and test_class_generic column. At run-time, the helper method combined the generic test class with the publication class to create a class that computed the property using the protocol parameters specified in each publication.

Ion channels placed onto chosen cell morphologies and conductances fitted to experimental distributions

Previously Published Olfactory Models

A literature search to identify each cell type’s electrical properties revealed a set of ion channels likely responsible for the behavior. Furthermore, a separate literature search was performed to identify previous computational models of each cell type. ModelDB was utilized to locate runnable NEURON simulator models. Results from the initial search of all olfactory bulb or cell models were manually inspected by following the instructions to run the models and then inspecting which cell types were modeled. If a model included multi-compartment morphologies, the morphology was extracted into .SWC files using the hoc2swc tool. The set of previous olfactory models can be found under [repo]/prev_ob_models.

Selection of Ion Channel Models

In addition to morphologies, the previous models included ion channel .MOD files. Many of the .MOD files were nearly identical copies of files from earlier publications. After tracing the publication geneology of each .MOD file, the most recent versions were identified. By examining the reasons why the channels were included, and comparing the electrical behavior seen in experimental literature, channels hypothesized to be responsible for the experimentally observed behaviors were selected for inclusion in this model. If an ion channel model did not include temperature sensitivity, equations to utilize the Q10 coefficient were added to the .MOD file. The ion channel models used in this model can be found under [repo]/prev_ob_models/Birgiolas2020/Mechanisms.

NEURON Cell Builder feature was used to import morphology SWC files and then insert channel mechanisms. The initial, unoptimized cell models can be found under [repo]/prev_ob_models/Birgiolas2020/Cells.

Parameter Fitting and Validation

A genetic optimization algorithm was used to identify the combinations of ion channel parameters that best reproduced experimentally observed electrical behaviors (minimized the error between model electrical property values and experimentally observed property means). The algorithm is defined in fitting.py and Jupyter notebooks used for fitting can be found under [repo]/notebooks (files that start with fitting).

The validation of the cell models against experimental data and comparison to previously published models can be found in all-model-validation-results.ipynb notebook.

The final, optimized models and their parameters are defined in [repo]/prev_ob_models/Birgiolas2020/isolated_cells.py (see MC1..5, GC1..5, TC1..5 classes).

Olfactory bulb layers were reconstructed from sagittal and coronal slices

Bulb Slices and Limitations

To place olfactory bulb cell models in their appropriate anatomical locations, a reconstruction of olfactory bulb layers was necessary. Allen Brain Atlas contains labeled sagital and coronal slices of full adult mouse olfactory bulbs. However, the Allen Atlas API, exposed only partialy reconstructed bulbs (frontal extremes were missing), and did not have labeled cell layer information. But, images of the full bulb slices with layers labeld were available through the online viewer. Furthermore, the sagital and coronal slices were different thickness, and the thiner coronal slices did not contain the the frontal bulb extreme region. Thus the coronal slice images were used to reconstruct the 3D shape of the mid and posterior parts of the olfactory bulb at higher resolution, and the thicker sagital slices were used to reconstruct the frontal extreme region.

Conversion of 2D Slices to 3D Layer Meshes

The slice images were downloaded and converted into SVG vector format using Inkscape. The SVG files were imported into Blender and positioned apart to match the thickness of the slices. The corresponding points of each neighboring layer were joined together using Blender’s Make Face (from edges). This process was repeated for each layer separately, until the 3D mesh of all the layers of the bulb were reconstructed. The number of polygons of the final mesh were reduced by using Blender’s Decimate Modifier and the Dyntopo (dynamic topology) sculpting tool.

The Blender file containing the reconstructed layers can be found under [repo]/blender-files (bulb-layers.blend and ob-gloms-fast.blend).

Virtual Slices and Generating Potential Positions for Cell Somas within each Layer

To generate positions within each layer where cell somas could be located, the reconstructed layers were filled with approximately evenly distributed points. This was done to enable the use of ‘virtual slices’. A virtual slice can be any 3D shape, but in this model, it was chosen to be a rectangular box. To construct a virtual slice, the slice shape is positioned over a desired part of the reconstructed layers. Any potential-soma points that are located within the parts of the layers that are contained within the virtual slice are selected to be used as potential cell placement locations. Points outside the virtual slice are ignored.

The points for possible cell locations can be found in ob-gloms-fast.blend as Blender’s Particle Systems of reconstructed layer objects ending with “Particles” (e.g. 0 GL Particles for potential positions for glomeruli, and 4 GRL Particles for positions of cells within the granule cell layer). The number of points in each layer was based on the cell and glomeruli count estimates in adult mice bulbs. The points were created using the Blue Noise Particles addon for Blender.

Fitted cell somas were placed within each cell type’s stereotypical laminar location

When constructing a network model, after a virtual slice was positioned over a desired area of the bulbar layers, the potential locations of somas and glomeruli that were contained by the virtual slice mesh were selected for use as soma/glomeruli center locations.

This is started by the [repo]/buil-dslice.py and most of the work is performed by the SliceBuilderBlender class.

Depending on the number of cells of each type specified in SliceBuilderBlender, a random subset of the positions contained within the virtual slice are used to place cell somas/glomeruli.

Mitral cell somas are placed in locations contained within the mitral cell layer. Tufted cells are placed within the external plexiform layer, and granule cells are placed within the granule cell layer.

When a mitral cell position is chosen, the closest glomerulus location is used to determine the minimum length that a mitral cell’s apical dendrite must have to be placed at that position (to satisfy the condition that mitral cell apical dendrites must terminate within the glomerular layer). A model is chosen at random from this set of mitral cell models that meet the apical dendrite criteria. If no mitral cell has a long-enough apical dendrite, the mitral cell model with the longest apical dendrite is chosen. The above rules are also followed for tufted cells, whose somas are placed within the external plexiform layer.

Similarly, granule cell somas are placed at locations within the virtual slice, however, a specific model for each location is selected at random from a set of granule cell models that have apical dendrites that terminate within the external plexiform layer (reach beyond the closest point of the external plexiform but do not reach into the glomerular layer.

Apical dendrites were rotated towards their stereotypical terminations

Once their somas are placed, cells of all types are rotated so that their apical dendrites point towards glomeruli closest to them. The cells are then rotated around the apical dendrite axis by a random degree. At this point, mitral and tufted cell apical dendrites point towards the glomeruli that are closest to the cell somas, and due to the cell morphology standardization step described earlier, the lateral dendrites lie in a plane that is approximately orthogonal to the apical dendrites. However, the apical dendrites of the selected mitral or tufted cell model might extend beyond the closest glomerulus. For this reason, the apical dendrites are rotated independently from the soma, to point towards glomeruli that lie at approximately the same distance as the length of the apical dendrite. This guarantees that the mitral and tufted cell dendrites terminate within the glomerular layer.

Mitral and tufted cell lateral dendrites were aligned with the curvature of reconstructed olfactory bulb layers

In experimental literature, mitral and tufted cell lateral dendrites tend to be confined to the external plexiform layer, and generally follow the local curvature of the layer. To replicate this phenomenon, after mitral and tufted cells were placed, and apical dendrites oriented, the lateral dendrites of mitral and tufted cells were confined to the external plexiform layer. Mitral cell lateral dendrites were confined to the inner 60% portion of the external plexiform layer, while tufted cell lateral dendrites were confined to the remaining outer 40% portion.

The final cell positions and morphology were saved for later instantiation in NEURON. The morphology transform files can be found under [repo]/slices/DorsalColumnSlice ([X]C.json files).

In the above visualization of the virtual slice used in this model, mitral cells are yellow/green, tufted cells are red, and granule cells are blue.

Reciprocal synapses were formed based on dendritic proximity between principal and granule cell dendrites

After all the cell dendrites were in their appropriate locations, reciprocal synapses were formed between mitral and granule cells and tufted and granule cells. Synapses were formed where granule cell dendrites approached mitral or tufted cell dendrites within 5 microns. At each location, an excitatory AMPA/MNDA synapse was placed onto the granule cell membrane, and an inhibitory GABA synapse was placed onto the tufted/mitral cell membrane. Granule cells that did not form any synapses were removed from the simulation. The synapse information was saved for later instantiation in NEURON and can be seen under [repo]/slices/DorsalColumnSlice (GCs__MCs.json and GCs__TCs.json files).

Gap junctions were formed between glomerular sibling principal cell tufted dendrites

Gap junctions were added between the apical tufts of glomerular sibling mitral cells. The same was done with glomerular sibling tufted cells. NEURON .setup_transfer() function was used to share the continuous current flows between the connected principal compartments. Within a glomerulus, the sibling pricipal cells were connected with each other in a round-robin fashion, where each pricipal cell was connected by gap junctions to two other siblings of same cell type. This procedure can be seen in model.py add_gap_junctions().

Glomeruli stimulated during odor experiments were mapped onto the model glomeruli

To stimulate glomeruli with realistic activation patterns, glomerular activations visible in optical glomerular imaging experiments done in Vincis et. al. (2012) were used to stimulate the network. These activations were previously digitized in Migliore et. al. (2014) model. To identify the analogous glomeruli, the Migliore model was imported into Blender using BlenderNEURON and the model’s glomeruli were registered with the reconstructed glomerular layer of this model. The glomerular locations that were closest to those in the Migliore model were used for odor input. The odors from Vincis and the corresponding glomeruli and their peak intensities were stored in the SQLite database in odor and odor_glom tables.

A model of glomerular input spikes was created to stimulate the model

Each tufted dendrite of a pricipal cell received input stimulation in the form of a train of spikes. The spike times for each cell were picked from a Gaussian distribution whose symmetric 99% range was scaled to correspond to the chosen inhalation duration (e.g. 100 ms). The spikes were delivered to excitatory synapses placed in each pricipal cell’s apical dendrites. This procedure can be seen in model.py stim_glom_segments()

An extracellular local field potential electrode was placed into the granule cell layer

Earlier work by Parasuram et. al. (2016) was extended to create a Python-based, MPI-compatible version of LFPsim. The new version, LFPsimpy was utilized in this model. Using the library an extracellular electrode was placed in approximately the same location in the granule cell layer as in Manabe & Mori (2013) which demonstrated the gamma fingerprint. The code to insert the electrode can be seen in model.py create_lfp_electrode()

NEURON+MPI simulations were performed and LFP signal analyzed using wavelet transform

The model was designed to utilize multiple processes using MPI and Parallel NEURON. During model initializtion, cells were distributed by placing them onto least utilized ranks. The number of cell compartments (segments) was used a proxy for cell complexity. Simulations were performed on Arizona State University computers and on Amazon Web Services instances. Most simulations were 1,800 ms long. The full model initialization and simulation code can be seen in model.py OlfactoryBulb class

During simulations, input spike times, soma membrane voltages of pricipal cells, and the LFP electrode signal were collected and saved at the end of the simulation for off-line analysis. The LFP signal was first band-pass filtered to 30-120 Hz, and then the pywavelets package was used to create wavelet spectrogram of the filtered signal. The final spectrogam was the average spectrogram of all sniffs during the simulation. The LFP Wavelet Analysis.ipynb notebook contains the wavelet decomposition code.

Network parameters were explored until the two-cluster gamma fingerprint was reproduced

Computational experiments were performed to demonstrate the mechanisms underlying the gamma fingerprint

  • experiments
    • silent network

    • Only MCs or TCs

    • Added GCs

    • Added Gap Junctions