.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/workflows/plot_simulate_somato.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_workflows_plot_simulate_somato.py>`
        to download the full example code. or to run this example in your browser via Binder

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_workflows_plot_simulate_somato.py:


================================================
04. From MEG sensor-space data to HNN simulation
================================================

This example demonstrates how to calculate an inverse solution of the median
nerve evoked response potential (ERP) in S1 from the MNE somatosensory dataset,
and then simulate a biophysical model network that reproduces the observed
dynamics. Note that we do not expound on how we came up with the sequence of
evoked drives used in this example, rather, we only demonstrate its
implementation. For those who want more background on the HNN model and the
process used to articulate the proximal and distal drives needed to simulate
evoked responses, see the `HNN ERP tutorial`_. The sequence of evoked drives
presented here is not part of a current publication but is motivated by prior
studies [1]_, [2]_.

.. GENERATED FROM PYTHON SOURCE LINES 17-23

.. code-block:: Python


    # Authors: Mainak Jas <mainakjas@gmail.com>
    #          Ryan Thorpe <ryan_thorpe@brown.edu>

    # sphinx_gallery_thumbnail_number = 2








.. GENERATED FROM PYTHON SOURCE LINES 24-29

First, we will import the packages needed for computing the inverse solution
from the MNE somatosensory dataset. `MNE`_ (and its dependency `NiBabel`_)
can be installed with ``pip install mne nibabel``
The somatosensory dataset can be downloaded by
importing ``somato`` from ``mne.datasets``.

.. GENERATED FROM PYTHON SOURCE LINES 29-36

.. code-block:: Python

    import os.path as op
    import matplotlib.pyplot as plt

    import mne
    from mne.datasets import somato
    from mne.minimum_norm import apply_inverse, make_inverse_operator








.. GENERATED FROM PYTHON SOURCE LINES 37-38

Now we set the the path of the ``somato`` dataset for subject ``'01'``.

.. GENERATED FROM PYTHON SOURCE LINES 38-47

.. code-block:: Python

    data_path = somato.data_path()
    subject = '01'
    task = 'somato'
    raw_fname = op.join(data_path, 'sub-{}'.format(subject), 'meg',
                        'sub-{}_task-{}_meg.fif'.format(subject, task))
    fwd_fname = op.join(data_path, 'derivatives', 'sub-{}'.format(subject),
                        'sub-{}_task-{}-fwd.fif'.format(subject, task))
    subjects_dir = op.join(data_path, 'derivatives', 'freesurfer', 'subjects')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Using default location ~/mne_data for somato...
    Creating ~/mne_data
    Downloading file 'MNE-somato-data.tar.gz' from 'https://osf.io/tp4sg/download?version=7' to '/home/circleci/mne_data'.

      0%|                                               | 0.00/611M [00:00<?, ?B/s]
      1%|▎                                     | 4.61M/611M [00:00<00:13, 44.2MB/s]
      1%|▌                                     | 9.03M/611M [00:00<00:30, 19.9MB/s]
      3%|█                                     | 16.5M/611M [00:00<00:17, 34.8MB/s]
      4%|█▍                                    | 23.5M/611M [00:00<00:13, 44.7MB/s]
      5%|█▊                                    | 29.3M/611M [00:00<00:12, 48.3MB/s]
      6%|██▏                                   | 35.8M/611M [00:00<00:10, 53.2MB/s]
      7%|██▊                                   | 44.3M/611M [00:00<00:09, 62.4MB/s]
      9%|███▎                                  | 53.5M/611M [00:01<00:07, 71.2MB/s]
     10%|███▉                                  | 62.3M/611M [00:01<00:07, 76.1MB/s]
     12%|████▍                                 | 71.7M/611M [00:01<00:06, 81.5MB/s]
     13%|█████                                 | 80.9M/611M [00:01<00:06, 84.7MB/s]
     15%|█████▌                                | 89.8M/611M [00:01<00:06, 85.9MB/s]
     16%|██████▏                               | 98.5M/611M [00:01<00:05, 85.5MB/s]
     18%|██████▊                                | 107M/611M [00:01<00:06, 73.8MB/s]
     19%|███████▎                               | 115M/611M [00:01<00:06, 71.5MB/s]
     20%|███████▊                               | 122M/611M [00:01<00:07, 68.3MB/s]
     21%|████████▎                              | 129M/611M [00:02<00:07, 62.4MB/s]
     22%|████████▋                              | 136M/611M [00:02<00:07, 60.7MB/s]
     23%|█████████                              | 142M/611M [00:02<00:07, 61.0MB/s]
     24%|█████████▍                             | 148M/611M [00:02<00:08, 53.7MB/s]
     26%|█████████▉                             | 156M/611M [00:02<00:07, 60.1MB/s]
     27%|██████████▎                            | 162M/611M [00:02<00:07, 58.3MB/s]
     28%|██████████▊                            | 169M/611M [00:02<00:07, 60.8MB/s]
     29%|███████████▏                           | 175M/611M [00:02<00:07, 56.2MB/s]
     30%|███████████▋                           | 182M/611M [00:02<00:07, 60.0MB/s]
     31%|████████████                           | 188M/611M [00:03<00:10, 39.0MB/s]
     32%|████████████▍                          | 194M/611M [00:03<00:09, 42.4MB/s]
     33%|████████████▊                          | 200M/611M [00:03<00:08, 46.5MB/s]
     34%|█████████████                          | 205M/611M [00:03<00:08, 47.8MB/s]
     35%|█████████████▌                         | 212M/611M [00:03<00:07, 52.1MB/s]
     36%|█████████████▉                         | 218M/611M [00:03<00:07, 55.3MB/s]
     37%|██████████████▍                        | 226M/611M [00:03<00:06, 60.8MB/s]
     39%|███████████████                        | 236M/611M [00:03<00:05, 71.1MB/s]
     40%|███████████████▌                       | 243M/611M [00:04<00:05, 71.5MB/s]
     41%|████████████████                       | 251M/611M [00:04<00:04, 75.2MB/s]
     42%|████████████████▌                      | 259M/611M [00:04<00:05, 66.7MB/s]
     44%|█████████████████                      | 267M/611M [00:04<00:04, 69.4MB/s]
     45%|█████████████████▍                     | 274M/611M [00:04<00:05, 66.3MB/s]
     46%|█████████████████▉                     | 281M/611M [00:04<00:04, 67.4MB/s]
     47%|██████████████████▍                    | 288M/611M [00:04<00:04, 67.8MB/s]
     48%|██████████████████▉                    | 296M/611M [00:04<00:04, 71.6MB/s]
     50%|███████████████████▍                   | 304M/611M [00:04<00:04, 72.8MB/s]
     51%|███████████████████▊                   | 311M/611M [00:05<00:04, 61.6MB/s]
     52%|████████████████████▎                  | 318M/611M [00:05<00:05, 56.5MB/s]
     53%|████████████████████▋                  | 324M/611M [00:05<00:05, 57.0MB/s]
     54%|█████████████████████                  | 329M/611M [00:05<00:05, 55.7MB/s]
     55%|█████████████████████▍                 | 335M/611M [00:05<00:05, 54.7MB/s]
     56%|█████████████████████▊                 | 341M/611M [00:05<00:05, 49.2MB/s]
     57%|██████████████████████▏                | 347M/611M [00:05<00:04, 52.9MB/s]
     58%|██████████████████████▌                | 354M/611M [00:05<00:04, 55.6MB/s]
     59%|██████████████████████▉                | 360M/611M [00:06<00:04, 52.2MB/s]
     60%|███████████████████████▍               | 367M/611M [00:06<00:04, 57.9MB/s]
     61%|███████████████████████▉               | 374M/611M [00:06<00:03, 61.2MB/s]
     63%|████████████████████████▍              | 382M/611M [00:06<00:03, 65.3MB/s]
     64%|████████████████████████▊              | 389M/611M [00:06<00:03, 67.3MB/s]
     65%|█████████████████████████▎             | 397M/611M [00:06<00:03, 69.9MB/s]
     66%|█████████████████████████▊             | 404M/611M [00:06<00:04, 46.7MB/s]
     67%|██████████████████████████▏            | 409M/611M [00:06<00:04, 45.7MB/s]
     68%|██████████████████████████▋            | 417M/611M [00:07<00:03, 52.5MB/s]
     70%|███████████████████████████            | 425M/611M [00:07<00:03, 57.9MB/s]
     71%|███████████████████████████▌           | 431M/611M [00:07<00:03, 58.9MB/s]
     72%|███████████████████████████▉           | 437M/611M [00:07<00:03, 53.4MB/s]
     73%|████████████████████████████▎          | 444M/611M [00:07<00:02, 56.2MB/s]
     74%|████████████████████████████▊          | 451M/611M [00:07<00:02, 59.8MB/s]
     75%|█████████████████████████████▏         | 457M/611M [00:07<00:02, 58.1MB/s]
     76%|█████████████████████████████▋         | 465M/611M [00:07<00:02, 64.7MB/s]
     78%|██████████████████████████████▏        | 473M/611M [00:07<00:01, 69.1MB/s]
     79%|██████████████████████████████▋        | 481M/611M [00:08<00:01, 67.4MB/s]
     80%|███████████████████████████████▏       | 487M/611M [00:08<00:02, 60.9MB/s]
     81%|███████████████████████████████▌       | 494M/611M [00:08<00:01, 62.0MB/s]
     82%|███████████████████████████████▉       | 500M/611M [00:08<00:01, 60.2MB/s]
     83%|████████████████████████████████▍      | 507M/611M [00:08<00:01, 62.6MB/s]
     84%|████████████████████████████████▊      | 514M/611M [00:08<00:01, 63.8MB/s]
     85%|█████████████████████████████████▎     | 521M/611M [00:08<00:01, 66.5MB/s]
     87%|█████████████████████████████████▊     | 529M/611M [00:08<00:01, 68.4MB/s]
     88%|██████████████████████████████████▎    | 537M/611M [00:08<00:00, 74.1MB/s]
     89%|██████████████████████████████████▉    | 546M/611M [00:08<00:00, 78.2MB/s]
     91%|███████████████████████████████████▍   | 554M/611M [00:09<00:00, 75.2MB/s]
     92%|███████████████████████████████████▉   | 562M/611M [00:09<00:00, 72.0MB/s]
     93%|████████████████████████████████████▎  | 569M/611M [00:09<00:00, 70.0MB/s]
     94%|████████████████████████████████████▊  | 576M/611M [00:09<00:00, 70.7MB/s]
     96%|█████████████████████████████████████▎ | 584M/611M [00:09<00:00, 72.9MB/s]
     97%|█████████████████████████████████████▊ | 591M/611M [00:09<00:00, 54.9MB/s]
     98%|██████████████████████████████████████▏| 598M/611M [00:09<00:00, 52.4MB/s]
     99%|██████████████████████████████████████▌| 604M/611M [00:09<00:00, 56.2MB/s]
      0%|                                               | 0.00/611M [00:00<?, ?B/s]
    100%|███████████████████████████████████████| 611M/611M [00:00<00:00, 1.79TB/s]
    Untarring contents of '/home/circleci/mne_data/MNE-somato-data.tar.gz' to '/home/circleci/mne_data'
    Attempting to create new mne-python configuration file:
    /home/circleci/.mne/mne-python.json
    Download complete in 16s (582.2 MB)




.. GENERATED FROM PYTHON SOURCE LINES 48-49

Then, we load the raw data and estimate the inverse operator.

.. GENERATED FROM PYTHON SOURCE LINES 49-69

.. code-block:: Python


    # Read and band-pass filter the raw data
    raw = mne.io.read_raw_fif(raw_fname, preload=True)
    l_freq, h_freq = 1, 40
    raw.filter(l_freq, h_freq)

    # Identify stimulus events associated with MEG time series in the dataset
    events = mne.find_events(raw, stim_channel='STI 014')

    # Define epochs within the time series
    event_id, tmin, tmax = 1, -.2, .17
    baseline = None
    epochs = mne.Epochs(raw, events, event_id, tmin, tmax, baseline=baseline,
                        reject=dict(grad=4000e-13, eog=350e-6), preload=True)

    # Compute the inverse operator
    fwd = mne.read_forward_solution(fwd_fname)
    cov = mne.compute_covariance(epochs)
    inv = make_inverse_operator(epochs.info, fwd, cov)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Opening raw data file /home/circleci/mne_data/MNE-somato-data/sub-01/meg/sub-01_task-somato_meg.fif...
        Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
    Ready.
    Reading 0 ... 269399  =      0.000 ...   897.077 secs...
    Filtering raw data in 1 contiguous segment
    Setting up band-pass filter from 1 - 40 Hz

    FIR filter parameters
    ---------------------
    Designing a one-pass, zero-phase, non-causal bandpass filter:
    - Windowed time-domain design (firwin) method
    - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
    - Lower passband edge: 1.00
    - Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
    - Upper passband edge: 40.00 Hz
    - Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
    - Filter length: 993 samples (3.307 s)

    [Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.1s
    [Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.4s
    [Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    1.0s
    [Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    1.8s
    111 events found on stim channel STI 014
    Event IDs: [1]
    Not setting metadata
    111 matching events found
    No baseline correction applied
    0 projection items activated
    Using data from preloaded Raw for 111 events and 112 original time points ...
    0 bad epochs dropped
    Reading forward solution from /home/circleci/mne_data/MNE-somato-data/derivatives/sub-01/sub-01_task-somato-fwd.fif...
        Reading a source space...
        [done]
        Reading a source space...
        [done]
        2 source spaces read
        Desired named matrix (kind = 3523) not available
        Read MEG forward solution (8155 sources, 306 channels, free orientations)
        Source spaces transformed to the forward solution coordinate frame
    Computing rank from data with rank=None
        Using tolerance 1.5e-08 (2.2e-16 eps * 306 dim * 2.2e+05  max singular value)
        Estimated rank (mag + grad): 306
        MEG: rank 306 computed from 306 data channels with 0 projectors
    /home/circleci/project/examples/workflows/plot_simulate_somato.py:66: RuntimeWarning: Something went wrong in the data-driven estimation of the data rank as it exceeds the theoretical rank from the info (306 > 64). Consider setting rank to "auto" or setting it explicitly as an integer.
      cov = mne.compute_covariance(epochs)
    Reducing data rank from 306 -> 306
    Estimating covariance using EMPIRICAL
    Done.
    Number of samples used : 12432
    [done]
    Converting forward solution to surface orientation
        No patch info available. The standard source space normals will be employed in the rotation to the local surface coordinates....
        Converting to surface-based source orientations...
        [done]
    Computing inverse operator with 306 channels.
        306 out of 306 channels remain after picking
    Selected 306 channels
    Creating the depth weighting matrix...
        204 planar channels
        limit = 7615/8155 = 10.004172
        scale = 5.17919e-08 exp = 0.8
    Applying loose dipole orientations to surface source spaces: 0.2
    Whitening the forward solution.
    Computing rank from covariance with rank=None
        Using tolerance 2e-12 (2.2e-16 eps * 306 dim * 29  max singular value)
        Estimated rank (mag + grad): 64
        MEG: rank 64 computed from 306 data channels with 0 projectors
        Setting small MEG eigenvalues to zero (without PCA)
    Creating the source covariance matrix
    Adjusting source covariance matrix.
    Computing SVD of whitened and weighted lead field matrix.
        largest singular value = 2.42284
        scaling factor to adjust the trace = 3.86104e+18 (nchan = 306 nzero = 242)




.. GENERATED FROM PYTHON SOURCE LINES 70-78

There are several methods to do source reconstruction. Some of the methods
such as MNE are distributed source methods whereas dipole fitting will
estimate the location and amplitude of a single current dipole. At the
moment, we do not offer explicit recommendations on which source
reconstruction technique is best for HNN. However, we do want our users
to note that the dipole currents simulated with HNN are assumed to be normal
to the cortical surface. Hence, using the option ``pick_ori='normal'`` is
appropriate.

.. GENERATED FROM PYTHON SOURCE LINES 78-85

.. code-block:: Python

    snr = 3.
    lambda2 = 1. / snr ** 2
    evoked = epochs.average()
    stc = apply_inverse(evoked, inv, lambda2, method='MNE',
                        pick_ori="normal", return_residual=False,
                        verbose=True)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Preparing the inverse operator for use...
        Scaled noise and source covariance from nave = 1 to nave = 111
        Created the regularized inverter
        The projection vectors do not apply to these channels.
        Created the whitener using a noise covariance matrix with rank 64 (242 small eigenvalues omitted)
    Applying inverse operator to "1"...
        Picked 306 channels from the data
        Computing inverse...
        Eigenleads need to be weighted ...
        Computing residual...
        Explained  86.1% variance
    [done]




.. GENERATED FROM PYTHON SOURCE LINES 86-88

To extract the primary response in primary somatosensory cortex (S1), we
create a label for the postcentral gyrus (S1) in source-space

.. GENERATED FROM PYTHON SOURCE LINES 88-94

.. code-block:: Python

    hemi = 'rh'
    label_tag = 'G_postcentral'
    label_s1 = mne.read_labels_from_annot(subject, parc='aparc.a2009s', hemi=hemi,
                                          regexp=label_tag,
                                          subjects_dir=subjects_dir)[0]





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Reading labels from parcellation...
       read 1 labels from /home/circleci/mne_data/MNE-somato-data/derivatives/freesurfer/subjects/01/label/rh.aparc.a2009s.annot




.. GENERATED FROM PYTHON SOURCE LINES 95-106

Visualizing the distributed S1 activation in reference to the geometric
structure of the cortex (i.e., plotted on a structural MRI) can help us
figure out how to orient the dipole. Note that in the HNN framework,
positive and negative deflections of a current dipole source correspond to
upwards (from deep to superficial) and downwards (from superficial to deep)
current flow, respectively. Uncomment the following code to open an
interactive 3D render of the brain and its surface activation (requires the
``pyvista`` python library). You should get 2 plots, the first showing the
post-central gyrus label from which the dipole time course was extracted and
the second showing MNE activation at 0.040 sec that resemble the following
images.

.. GENERATED FROM PYTHON SOURCE LINES 106-116

.. code-block:: Python

    '''
    Brain = mne.viz.get_brain_class()
    brain_label = Brain(subject, hemi, 'white', subjects_dir=subjects_dir)
    brain_label.add_label(label_s1, color='green', alpha=0.9)
    stc_label = stc.in_label(label_s1)
    brain = stc_label.plot(subjects_dir=subjects_dir, hemi=hemi, surface='white',
                           view_layout='horizontal', initial_time=0.04,
                           backend='pyvista')
    '''





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    "\nBrain = mne.viz.get_brain_class()\nbrain_label = Brain(subject, hemi, 'white', subjects_dir=subjects_dir)\nbrain_label.add_label(label_s1, color='green', alpha=0.9)\nstc_label = stc.in_label(label_s1)\nbrain = stc_label.plot(subjects_dir=subjects_dir, hemi=hemi, surface='white',\n                       view_layout='horizontal', initial_time=0.04,\n                       backend='pyvista')\n"



.. GENERATED FROM PYTHON SOURCE LINES 117-120

|mne_label_fig|

|mne_activity_fig|

.. GENERATED FROM PYTHON SOURCE LINES 122-130

Now we extract the representative time course of dipole activation in our
labeled brain region using ``mode='pca_flip'`` (see `this MNE-python
example`_ for more details). Note that the most prominent component of the
median nerve response occurs in the posterior wall of the central sulcus at
~0.040 sec. Since the dipolar activity here is negative, we orient the
extracted waveform so that the deflection at ~0.040 sec is pointed downwards.
Thus, the ~0.040 sec deflection corresponds to current flow traveling from
superficial to deep layers of cortex.

.. GENERATED FROM PYTHON SOURCE LINES 130-142

.. code-block:: Python

    flip_data = stc.extract_label_time_course(label_s1, inv['src'],
                                              mode='pca_flip')
    dipole_tc = -flip_data[0] * 1e9

    plt.figure()
    plt.plot(1e3 * stc.times, dipole_tc, 'ro--')
    plt.xlabel('Time (ms)')
    plt.ylabel('Current Dipole (nAm)')
    plt.xlim((0, 170))
    plt.axhline(0, c='k', ls=':')
    plt.show()




.. image-sg:: /auto_examples/workflows/images/sphx_glr_plot_simulate_somato_001.png
   :alt: plot simulate somato
   :srcset: /auto_examples/workflows/images/sphx_glr_plot_simulate_somato_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Extracting time courses for 1 labels (mode: pca_flip)




.. GENERATED FROM PYTHON SOURCE LINES 143-145

Now, let us try to simulate the same with ``hnn-core``. We read in the
network parameters from ``N20.json`` and instantiate the network.

.. GENERATED FROM PYTHON SOURCE LINES 145-154

.. code-block:: Python


    import hnn_core
    from hnn_core import simulate_dipole, jones_2009_model
    from hnn_core import average_dipoles, JoblibBackend

    hnn_core_root = op.dirname(hnn_core.__file__)
    params_fname = op.join(hnn_core_root, 'param', 'N20.json')
    net = jones_2009_model(params_fname)








.. GENERATED FROM PYTHON SOURCE LINES 155-163

To simulate the source of the median nerve evoked response, we add a
sequence of synchronous evoked drives: 1 proximal, 2 distal, and 1 final
proximal drive. In order to understand the physiological implications of
proximal and distal drive as well as the general process used to articulate
a sequence of exogenous drive for simulating evoked responses, see the
`HNN ERP tutorial`_. Note that setting ``n_drive_cells=1`` and
``cell_specific=False`` creates a drive with synchronous input across cells
in the network.

.. GENERATED FROM PYTHON SOURCE LINES 163-220

.. code-block:: Python


    # Early proximal drive
    weights_ampa_p = {'L2_basket': 0.0036, 'L2_pyramidal': 0.0039,
                      'L5_basket': 0.0019, 'L5_pyramidal': 0.0020}
    weights_nmda_p = {'L2_basket': 0.0029, 'L2_pyramidal': 0.0005,
                      'L5_basket': 0.0030, 'L5_pyramidal': 0.0019}
    synaptic_delays_p = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
                         'L5_basket': 1.0, 'L5_pyramidal': 1.0}

    net.add_evoked_drive(
        'evprox1', mu=21., sigma=4., numspikes=1, location='proximal',
        n_drive_cells=1, cell_specific=False, weights_ampa=weights_ampa_p,
        weights_nmda=weights_nmda_p, synaptic_delays=synaptic_delays_p,
        event_seed=276)

    # Late proximal drive
    weights_ampa_p = {'L2_basket': 0.003, 'L2_pyramidal': 0.0039,
                      'L5_basket': 0.004, 'L5_pyramidal': 0.0020}
    weights_nmda_p = {'L2_basket': 0.001, 'L2_pyramidal': 0.0005,
                      'L5_basket': 0.002, 'L5_pyramidal': 0.0020}
    synaptic_delays_p = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
                         'L5_basket': 1.0, 'L5_pyramidal': 1.0}

    net.add_evoked_drive(
        'evprox2', mu=134., sigma=4.5, numspikes=1, location='proximal',
        n_drive_cells=1, cell_specific=False, weights_ampa=weights_ampa_p,
        weights_nmda=weights_nmda_p, synaptic_delays=synaptic_delays_p,
        event_seed=276)

    # Early distal drive
    weights_ampa_d = {'L2_basket': 0.0043, 'L2_pyramidal': 0.0032,
                      'L5_pyramidal': 0.0009}
    weights_nmda_d = {'L2_basket': 0.0029, 'L2_pyramidal': 0.0051,
                      'L5_pyramidal': 0.0010}
    synaptic_delays_d = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
                         'L5_pyramidal': 0.1}

    net.add_evoked_drive(
        'evdist1', mu=32., sigma=2.5, numspikes=1, location='distal',
        n_drive_cells=1, cell_specific=False, weights_ampa=weights_ampa_d,
        weights_nmda=weights_nmda_d, synaptic_delays=synaptic_delays_d,
        event_seed=277)

    # Late distal drive
    weights_ampa_d = {'L2_basket': 0.0041, 'L2_pyramidal': 0.0019,
                      'L5_pyramidal': 0.0018}
    weights_nmda_d = {'L2_basket': 0.0032, 'L2_pyramidal': 0.0018,
                      'L5_pyramidal': 0.0017}
    synaptic_delays_d = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
                         'L5_pyramidal': 0.1}

    net.add_evoked_drive(
        'evdist2', mu=84., sigma=4.5, numspikes=1, location='distal',
        n_drive_cells=1, cell_specific=False, weights_ampa=weights_ampa_d,
        weights_nmda=weights_nmda_d, synaptic_delays=synaptic_delays_d,
        event_seed=275)








.. GENERATED FROM PYTHON SOURCE LINES 221-224

Now we run the simulation over 2 trials so that we can plot the average
aggregate dipole. For a better match to the empirical waveform, set
``n_trials`` to be >=25.

.. GENERATED FROM PYTHON SOURCE LINES 224-229

.. code-block:: Python

    n_trials = 2
    # n_trials = 25
    with JoblibBackend(n_jobs=2):
        dpls = simulate_dipole(net, tstop=170., n_trials=n_trials)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Joblib will run 2 trial(s) in parallel by distributing trials over 2 jobs.




.. GENERATED FROM PYTHON SOURCE LINES 230-237

Since the model is a reduced representation of the larger network
contributing to the response, the model response is noisier than it would be
in the net activity from a larger network where these effects are averaged
out, and the dipole amplitude is smaller than the recorded data. The
post-processing steps of smoothing and scaling the simulated dipole response
allow us to more accurately approximate the true signal responsible for the
recorded macroscopic evoked response [1]_, [2]_.

.. GENERATED FROM PYTHON SOURCE LINES 237-243

.. code-block:: Python

    dpl_smooth_win = 20
    dpl_scalefctr = 12
    for dpl in dpls:
        dpl.smooth(dpl_smooth_win)
        dpl.scale(dpl_scalefctr)








.. GENERATED FROM PYTHON SOURCE LINES 244-246

Finally, we plot the driving spike histogram, empirical and simulated median
nerve evoked response waveforms, and output spike histogram.

.. GENERATED FROM PYTHON SOURCE LINES 246-258

.. code-block:: Python

    fig, axes = plt.subplots(3, 1, sharex=True, figsize=(6, 6),
                             constrained_layout=True)
    net.cell_response.plot_spikes_hist(ax=axes[0],
                                       spike_types=['evprox', 'evdist'],
                                       show=False)
    axes[1].axhline(0, c='k', ls=':', label='_nolegend_')
    axes[1].plot(1e3 * stc.times, dipole_tc, 'r--')
    average_dipoles(dpls).plot(ax=axes[1], show=False)
    axes[1].legend(['MNE label average', 'HNN simulation'])
    axes[1].set_ylabel('Current Dipole (nAm)')
    net.cell_response.plot_spikes_raster(ax=axes[2])




.. image-sg:: /auto_examples/workflows/images/sphx_glr_plot_simulate_somato_002.png
   :alt: Aggregate (L2/3 + L5)
   :srcset: /auto_examples/workflows/images/sphx_glr_plot_simulate_somato_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Figure size 600x600 with 3 Axes>



.. GENERATED FROM PYTHON SOURCE LINES 259-270

References
----------
.. [1] Jones, S. R., Pritchett, D. L., Stufflebeam, S. M., Hämäläinen, M.
   & Moore, C. I. Neural correlates of tactile detection: a combined
   magnetoencephalography and biophysically based computational modeling
   study. J. Neurosci. 27, 10751–10764 (2007).
.. [2] Neymotin SA, Daniels DS, Caldwell B, McDougal RA, Carnevale NT,
   Jas M, Moore CI, Hines ML, Hämäläinen M, Jones SR. Human Neocortical
   Neurosolver (HNN), a new software tool for interpreting the cellular and
   network origin of human MEG/EEG data. eLife 9, e51214 (2020).
   https://doi.org/10.7554/eLife.51214

.. GENERATED FROM PYTHON SOURCE LINES 272-282

.. LINKS

.. _MNE: https://mne.tools/
.. _NiBabel: https://nipy.org/nibabel/
.. _HNN ERP tutorial: https://jonescompneurolab.github.io/hnn-tutorials/erp/erp
.. _this MNE-python example: https://mne.tools/stable/auto_examples/inverse/plot_label_source_activations.html
.. |mne_label_fig| image:: https://user-images.githubusercontent.com/20212206/106524603-cfe75c80-64b0-11eb-9607-3415195c3e7a.png
  :width: 400
.. |mne_activity_fig| image:: https://user-images.githubusercontent.com/20212206/106524542-b514e800-64b0-11eb-835e-497454e75eb9.png
  :width: 400


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 55.951 seconds)


.. _sphx_glr_download_auto_examples_workflows_plot_simulate_somato.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: binder-badge

      .. image:: images/binder_badge_logo.svg
        :target: https://mybinder.org/v2/gh/jonescompneurolab/hnn-core/gh-pages?filepath=v0.4/notebooks/auto_examples/workflows/plot_simulate_somato.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_simulate_somato.ipynb <plot_simulate_somato.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_simulate_somato.py <plot_simulate_somato.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_simulate_somato.zip <plot_simulate_somato.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_