.. _neuron_tut:
`NEURON` simulator
******************

`NEURON <https://neuron.yale.edu/neuron/>`_ is a simulation environment to model neurons and networks thereof. `NEURON` itself is rather complex (neurons are complex things after all) and fairly low-level which results in lots of boiler plate code. There are some libraries(e.g. `NetPyNE <http://netpyne.org/>`_) that wrap `NEURON` and provide a higher-level interface to facilitate building models. In my experience these are typically geared towards creating models based on probabilities (e.g. "create 100 neurons with a 10% chance to be connected to another neuron") rather than the well defined morphology/connectivity you get from e.g. connectomes.

`navis` does *not* try to emulate a full simulation suite but tries to fill a gap for people wanting to use non-probabalistic data (e.g. from connectomics) by providing an entry point for you to get started with some simple models and take it from there. At this point there are two types of models:

1. :class:`~navis.interfaces.neuron.comp.CompartmentModel` for modeling individual neurons 
2. :class:`~navis.interfaces.neuron.network.PointNetwork` for modeling networks from point

Compartment models
------------------
A :class:`~navis.interfaces.neuron.comp.CompartmentModel` represents a single neuron (although you can connect multiple of these neurons) and is constructed from a skeleton (i.e. :class:`navis.TreeNeuron`, see also :func:`navis.conversion.mesh2skeleton`). 

.. code:: ipython3

    import navis
    import neuron
    
    import navis.interfaces.neuron as nrn

.. code:: ipython3

    # Load one of the example neurons (a Drosophila projection neuron from the hemibrain connectome)
    # Note the conversion to microns!
    n = navis.example_neurons(1).convert_units('um')
    
    # Here we manually corrected the soma 
    n.soma = 20
    
    # Reroot to the soma 
    n.reroot(n.soma, inplace=True)
    
    # Create the compartment model
    cmp = nrn.CompartmentModel(n, res=10)

`NEURON` compartment models are effectively collections of connected linear segments. Each of these segments can have its own set of properties and mechanisms. Here, we will use some biophysical properties from `Tobin et al. (2017) <https://elifesciences.org/articles/24838>`_:

.. code:: ipython3

    # Set the specific axial resistivity for the entire neuron in Ohm cm
    cmp.Ra = 266.1
    
    # Set the specific membran capacitance in mF / cm**2
    cmp.cm = 0.8
    
    # Add passive membran properties for the entire neuron
    cmp.insert('pas',
               g=1/20800,  # specific leakage conductance = 1/Rm; Rm = specific membran resistance in Ohm cm**2
               e=-60,      # leakage reverse potential
               )

At this point we already have a passive model of our *Drosophila* projection neuron. Next, we will determine what's axon and what's dendrite, and add a Hodgkins & Huxley mechanism so we have some actual action potentials to observe. 

This will also illustrate one of the features of ``CompartmentModels``: we keep a correspondence to the original skeleton via node (and connector) IDs. This allows you to modify, stimulate at or record from parts of the compartment model by using node IDs.

.. code:: ipython3

    # Label axon/dendrite
    navis.split_axon_dendrite(n, label_only=True, cellbodyfiber='soma')
    
    n.nodes.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>node_id</th>
          <th>label</th>
          <th>x</th>
          <th>y</th>
          <th>z</th>
          <th>radius</th>
          <th>parent_id</th>
          <th>type</th>
          <th>compartment</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>1</td>
          <td>0</td>
          <td>126.272003</td>
          <td>298.000000</td>
          <td>224.496017</td>
          <td>0.080000</td>
          <td>2</td>
          <td>end</td>
          <td>cellbodyfiber</td>
        </tr>
        <tr>
          <th>1</th>
          <td>2</td>
          <td>0</td>
          <td>126.112007</td>
          <td>297.840027</td>
          <td>224.656006</td>
          <td>0.146274</td>
          <td>3</td>
          <td>slab</td>
          <td>cellbodyfiber</td>
        </tr>
        <tr>
          <th>2</th>
          <td>3</td>
          <td>0</td>
          <td>125.952003</td>
          <td>297.520020</td>
          <td>224.976013</td>
          <td>0.277771</td>
          <td>4</td>
          <td>slab</td>
          <td>cellbodyfiber</td>
        </tr>
        <tr>
          <th>3</th>
          <td>4</td>
          <td>0</td>
          <td>125.952003</td>
          <td>297.200012</td>
          <td>225.616013</td>
          <td>0.277771</td>
          <td>5</td>
          <td>slab</td>
          <td>cellbodyfiber</td>
        </tr>
        <tr>
          <th>4</th>
          <td>5</td>
          <td>0</td>
          <td>125.632004</td>
          <td>297.040009</td>
          <td>225.936005</td>
          <td>0.277771</td>
          <td>6</td>
          <td>slab</td>
          <td>cellbodyfiber</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    # A quick visualization
    fig, ax = navis.plot2d(n, color_by='compartment', palette='tab10', lw=1.5)
    
    # Label compartments 
    import seaborn as sns
    for l in n.nodes.compartment.unique():
        loc = n.nodes.loc[n.nodes.compartment == l, ['x', 'y', 'z']].values[-1]
        ax.text(loc[0] + 10, loc[1], loc[2], l)    
    
    ax.elev = -110



.. image:: neuron_files/neuron_7_0.png


A little excursion about how `NEURON` represents neurons before we add the Hodgkins & Huxley (HH) mechanism: neurons consist of linear "*sections*" while individual (continuous) positions along each section are called "*segments*". That distinction is important because our skeleton nodes correspond to *segments* but a mechanism such as HH is applied for entire *sections*. 

.. code:: ipython3

    # Collect axon nodes
    axon_nodes = n.nodes.loc[n.nodes.compartment.isin(['axon', 'linker']), 'node_id'].values
    
    # Get the sections for the given nodes
    axon_secs = list(set(cmp.get_node_section(axon_nodes)))
    
    # Insert HH mechanism at the given sections
    cmp.insert('hh', subset=axon_secs)

Next we actually have to do something with our compartment model. For that we will add three voltage recorders: one at the soma, one at the base of the dendrites and one at the tip of the axon.

.. code:: ipython3

    # Recordings and stimulations work via segments (as opposed to sections for mechanisms) -> hence we can use node IDs directly here
    
    # Let's determine the tip of the axon and base of the dendrites programmatically using the 
    # geodesic distance
    dists = navis.geodesic_matrix(n, from_=n.soma)
    
    # Sort by distance from soma 
    dists = dists.iloc[0].sort_values()
    
    dists.head(10)




.. parsed-literal::

    20      0.000000
    21      0.455753
    19      0.633259
    2409    0.836538
    22      1.095767
    18      1.113258
    23      1.415775
    2410    1.474030
    24      1.692905
    17      2.032377
    Name: 20, dtype: float64



.. code:: ipython3

    # Find the closest "dendrite" and the most distal "axon" node
    dend = n.nodes[n.nodes.compartment == 'dendrite'].node_id.values
    dend_base = dists.index[dists.index.isin(dend)][0]
    print(f'Node at the base of the dendrites: {dend_base}')
    
    axo = n.nodes[n.nodes.compartment == 'axon'].node_id.values
    axo_tip = dists.index[dists.index.isin(axo)][-1]
    print(f'Node at the tip of the axon: {axo_tip}')


.. parsed-literal::

    Node at the base of the dendrites: 467
    Node at the tip of the axon: 4384


.. code:: ipython3

    # Add voltage recordings
    cmp.add_voltage_record(dend_base, label='dendrite_base')
    cmp.add_voltage_record(axo_tip, label='axon_tip')
    cmp.add_voltage_record(n.soma, label='soma')

Last but not least we need to provide some input to our neuron otherwise it will just sit there doing nothing. We can add a current injection, trigger some synaptic currents or add a leak current. Let's simulate some synaptic inputs:

.. code:: ipython3

    # Get dendritic postsynapses 
    post = n.postsynapses[n.postsynapses.compartment == 'dendrite']
    post.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>connector_id</th>
          <th>node_id</th>
          <th>type</th>
          <th>x</th>
          <th>y</th>
          <th>z</th>
          <th>roi</th>
          <th>confidence</th>
          <th>compartment</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>137</th>
          <td>137</td>
          <td>3736</td>
          <td>post</td>
          <td>135.232</td>
          <td>291.104</td>
          <td>210.840</td>
          <td>AL(R)</td>
          <td>0.906454</td>
          <td>dendrite</td>
        </tr>
        <tr>
          <th>138</th>
          <td>138</td>
          <td>635</td>
          <td>post</td>
          <td>132.872</td>
          <td>284.448</td>
          <td>196.056</td>
          <td>AL(R)</td>
          <td>0.928960</td>
          <td>dendrite</td>
        </tr>
        <tr>
          <th>139</th>
          <td>139</td>
          <td>4124</td>
          <td>post</td>
          <td>131.432</td>
          <td>277.048</td>
          <td>200.960</td>
          <td>AL(R)</td>
          <td>0.605084</td>
          <td>dendrite</td>
        </tr>
        <tr>
          <th>140</th>
          <td>140</td>
          <td>3176</td>
          <td>post</td>
          <td>124.984</td>
          <td>286.568</td>
          <td>216.480</td>
          <td>AL(R)</td>
          <td>0.922643</td>
          <td>dendrite</td>
        </tr>
        <tr>
          <th>219</th>
          <td>219</td>
          <td>3061</td>
          <td>post</td>
          <td>128.312</td>
          <td>276.848</td>
          <td>217.272</td>
          <td>AL(R)</td>
          <td>0.981077</td>
          <td>dendrite</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    # Here we will open successively more synapses over 5 stimulations
    for i in range(5):    
        # Onset for this stimulation
        start = 50 + i * 200
        # Number of synapses to activate
        n_syn = i * 5
        cmp.add_synaptic_current(where=post.node_id.unique()[0:n_syn],
                                 start=start,
                                 max_syn_cond=.1,
                                 rev_pot=-10)

.. code:: ipython3

    # Now we can run our simulation for 1000ms
    # (this is equivalent to neuron.h.finitialize + neuron.h.continuerun)
    cmp.run_simulation(1000, v_init=-60)

The compartment model has a quick & dirty way of plotting the results:

.. code:: ipython3

    # Plot the results
    axes = cmp.plot_results()



.. image:: neuron_files/neuron_19_0.png


.. code:: ipython3

    # Plot again and zoom in on one spike
    axes = cmp.plot_results()
    axes[0].set_xlim(240, 280)




.. parsed-literal::

    (240.0, 280.0)




.. image:: neuron_files/neuron_20_1.png


As you can see we get a nice depolarization at the base of the dendrites which elicits an action potential that we can measure in the tips of the axon. Because in our model the cell body fiber (i.e. the neurite that connects the soma to the base of the dendrites) is passive, the depolarization of a single spike attenuates before it reaches the soma.

Alternatively, you can access the recorded values directly like so: 

.. code:: ipython3

    cmp.records




.. parsed-literal::

    {'v': {'dendrite_base': Vector[0], 'axon_tip': Vector[1], 'soma': Vector[2]}}



.. code:: ipython3

    cmp.records['v']['dendrite_base'].as_numpy()




.. parsed-literal::

    array([-60.        , -60.04157854, -60.07494932, ..., -63.52188478,
           -63.52188478, -63.52188478])



Next, let's try to simulate some noisy input where the presynaptic neuron spikes multiple times:

.. code:: ipython3

    # First we need to reset our model (by re-assigning `cmp` the old model will be garbage-collected)
    cmp = nrn.CompartmentModel(n, res=10)
    
    # Set properties and mechanisms
    cmp.Ra, cmp.cm = 266.1, 0.8
    cmp.insert('pas', g=1/20800, e=-60)
    axon_secs = list(set(cmp.get_node_section(axon_nodes)))
    cmp.insert('hh', subset=axon_secs)
    
    # Add recording
    cmp.add_voltage_record(dend_base, label='dendrite_base')
    cmp.add_voltage_record(axo_tip, label='axon_tip')
    cmp.add_voltage_record(n.soma, label='soma')
    
    # Also add a spike counter at the axon 
    cmp.add_spike_detector(axo_tip, label='axon_tip')
    
    # Now add a noisy preinput to say 20 dendritic postsynapses
    post = n.connectors[(n.connectors.compartment == 'dendrite') & (n.connectors.type == 'post')]
    cmp.add_synaptic_input(post.node_id.unique()[0:20],
                           spike_no=20,  # produce 20 presynaptic spikes 
                           spike_int=50, # with an average interval of 50ms: 20 * 50ms = over 1s
                           spike_noise=1, # very noisy!
                           cn_weight=0.04)
    # Run for 1s
    cmp.run_simulation(1000, v_init=-60)

.. code:: ipython3

    # Plot results
    axes = cmp.plot_results()
    axes[1].set_ylabel('spikes [Hz]')




.. parsed-literal::

    Text(0, 0.5, 'spikes [Hz]')




.. image:: neuron_files/neuron_26_1.png


Note how we still don't see a depolarization in the soma? While that might be a genuine biological feature of this neuron, I suspect there is something wrong with the radii along the cell body fiber - perhaps a pinch point somewhere? This just illustrates that good skeletons are paramount and you should critically inspect the results of your models.


Many methods in :class:`~navis.interfaces.neuron.CompartmentModel` try to use sensible defaults to make sure that you get some sort of effect. That said, it's advisable that you adjust parameters as you fit your model to real world data. Check the help to see what you can do:

.. code:: ipython3

    # Try this for example
    cmp.add_synaptic_input?

Networks
--------
While you can link together multiple compartment models to simulate networks this quickly becomes prohibitively slow so run. For larger networks it can be sufficient to model each neuron as a single "point process". :class:`~navis.interfaces.neuron.PointNetwork` lets you quickly create such a network from an edge list. 

In this tutorial we will use toy data but it is just as straight forward to plugin real data:

.. code:: ipython3

    # First create a small 3 way network where one of the neurons (B) is inhibitory
    import pandas as pd 
    edges = pd.DataFrame([])
    edges['source'] = ['A', 'B']
    edges['target'] = ['C', 'C']
    edges['weight'] = [0.5, -1]
    edges




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>source</th>
          <th>target</th>
          <th>weight</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>A</td>
          <td>C</td>
          <td>0.5</td>
        </tr>
        <tr>
          <th>1</th>
          <td>B</td>
          <td>C</td>
          <td>-1.0</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    # Next initialize network from edge list 
    net = nrn.PointNetwork.from_edge_list(edges, model='IntFire4')
    net




.. parsed-literal::

    PointNetwork<neurons=3,edges=2>



So far, our network won't do anything because it doesn't have any input.
Let's add an input to neurons A and B, and try that out:

.. code:: ipython3

    # Add the stimulus
    net.add_stimulus('A', start=100, stop=1000, frequency=100, randomness=0)
    net.add_stimulus('B', start=600, stop=800, frequency=100, randomness=0)
    
    # Run simulation
    net.run_simulation(1000)

.. code:: ipython3

    # Plot 
    ax = net.plot_raster(backend='matplotlib', label=True)
    
    ax.set_xlim(0, 1000)




.. parsed-literal::

    (0.0, 1000.0)




.. image:: neuron_files/neuron_34_1.png


This toy example worked quite well: spikes in neuron `A` elicit occasional spikes in `C` via temporal summation. Activity in the inhibitiory neuron `B` hyperpolarizes `C` and it stops firing until well after activity in `B` has ceased. 

The `NEURON` interface is a very recent addition and it might well change in the future (or become its own package). Functionality is also still limited and while I don't intend to write a feature-complete wrapper for `NEURON`, I do welcome feature requests or contributions on Github.


Links
-----

- `Model DB <https://senselab.med.yale.edu/ModelDB/>`_ contains various published NEURON models and mechanisms
- `NetPyNE <http://netpyne.org/>`_ wraps NEURON and provides high-level syntax to create models