<jupyter_start><jupyter_text>Figure 2 - extracellular simulationsThis notebook performs a biophysical simulation of Cell 1 (neuromorpho ID: 561096006) to show how the model reproduces action potential propagation over axon branches.<jupyter_code>import neuron
from math import sin, cos
import numpy as np
import LFPy
import MEAutility as mu
import matplotlib.pyplot as plt
import neuroplotlib as npl
from pathlib import Path
import sys
import os
from pprint import pprint
from axon_velocity.models import insert_biophysics, insert_simple_biophysics, \
    get_default_biophysics_params, planarize_swc, save_cell, create_mea_probe, center_cell_xy
from axon_velocity import plot_amplitude_map, plot_peak_latency_map
%matplotlib widget
save_fig = False
fig_folder =  Path('figures') / "figure2"
fig_folder.mkdir(exist_ok=True, parents=True)
try:
    import neuron
except:
    print('NEURON is not installed.')

mechanism_folder = Path('..') / 'simulations' / 'mechanisms'

if not neuron.load_mechanisms(str(mechanism_folder)):
    print('Compile mod files in the mechanisms/ folder: from the mechanisms/ folder, run nrnivmodl')
# simple biophysiscs: dendrite - pas / soma/axon HH
# "complex" biophysics: dendrite - pas / soma - na + kv1 / axon - nax + kv1 
simple_biophysics = False
params_dict = get_default_biophysics_params()
pprint(params_dict)<jupyter_output><empty_output><jupyter_text>The `planar` variable decides wheter the z-axis is compressed (similar to a cell culture - `planar=True`) or the original morphology is used (`planar=False`).<jupyter_code>planar = True
z_offset = 10 # distance between cell plane and mea plane
zspan = 0
morphology_dir = Path('..') / 'simulations' / 'neuromorpho' / 'allen_cell_types'

morph_id = '561096006'
original_morphology_path = [m for m in morphology_dir.iterdir() if not 
                            m.name.startswith('.') and morph_id in str(m)][0]
if planar:
    morphology_path = planarize_swc(original_morphology_path, span_um=zspan)
else:
    morphology_path = original_morphology_path
cell = LFPy.Cell(str(morphology_path), v_init=params_dict['v_init'], celsius=params_dict['celsius'],
                 Ra=params_dict['ra'], cm=params_dict['cm'], pt3d=True)
# center in the xy plane
center_cell_xy(cell)
npl.plot_neuron(cell, plane='xy', color_axon='g')<jupyter_output><empty_output><jupyter_text>Insert cell biophysicsHere we make the cell active by inserting biophysical mechanisms.<jupyter_code>if simple_biophysics:
    insert_simple_biophysics(cell)
else:
    insert_biophysics(cell, params_dict)<jupyter_output><empty_output><jupyter_text>Stimulating the cellWe can now add some stimulation. The stimulation can be a current clamp `iclamp` or synaptic inputs `syn`. The `stim_point` is where the cell will be stimulated (the closest cell segment to the `stim_point` is used).<jupyter_code>stim = 'syn' # or syn
# stimulate on the soma
stim_idx = cell.somaidx

syn_input_times = np.arange(2, 5)

syn_params = {'idx' : stim_idx,
              'e' : 0,                                # reversal potential
              'syntype' : 'ExpSyn',                   # synapse type
              'tau' : 2,                              # syn. time constant ms
              'weight' : 0.05,                         # syn. weight
              'record_current' : True                 # syn. current record
    }
clamp_params = {'idx' : stim_idx,
                'pptype' : 'IClamp',                   # IClamp point process
                'dur' : 300,                            # dur in ms
                'amp' : 2,                             # amp in nA
                'delay' : 5                            # delay in ms
    }
if stim == 'syn':
    synapse = LFPy.Synapse(cell, **syn_params)
    synapse.set_spike_times(np.array(syn_input_times))
else:
    clamp = LFPy.StimIntElectrode(cell=cell, **clamp_params)
shift = z_offset
    
print(f"z-position of MEA: {shift}")<jupyter_output><empty_output><jupyter_text>Define extracellular electrodesLet's now define the extracellular electrodes using the [MEAutility](https://meautility.readthedocs.io/en/latest/) package.<jupyter_code>mea_dim = 100  # n rows x n cols
mea_pitch = 17.5  # rows and cols pitch
elec_size = 5

hdmea = create_mea_probe(pitch=mea_pitch, dim=mea_dim, elec_size=elec_size, z_offset=z_offset)

electrode = LFPy.RecExtElectrode(cell, probe=hdmea, n=10)<jupyter_output><empty_output><jupyter_text>Run the simulationBy passing the `electrode` argument `LFPy` also computes extracellular potentials. The `rec_vmem` argument allows to measure the membrane potenrtial at all segments.<jupyter_code>cell.simulate(probes=[electrode], rec_vmem=True)
eap = electrode.data * 1000  # mV --> uV
# cutout single template
fs = 1 / cell.dt
ms_before = 2
ms_after = 10

min_chan, min_idx = np.unravel_index(np.argmin(eap), eap.shape)
1 / cell.dt
eap_cut = eap[:, min_idx - int(ms_before * fs): min_idx + int(ms_after * fs)]
intra_tidxs = slice(min_idx - int(ms_before * fs), min_idx + int(ms_after * fs))<jupyter_output><empty_output><jupyter_text>Plots1. Plot MEA and overlaid morphology<jupyter_code>fig1, ax1 = plt.subplots(1, 1)

ax1 = mu.plot_probe(hdmea, type='planar', ax=ax1)
npl.plot_neuron(cell, plane='xy', color='k', ax=ax1)
ax1.axis('off')

# add scalebar
ax1.plot([-870, -770], [-920, -920], color='k', lw=2)
ax1.text(-880, -1000, "100 $\mu$m", fontsize=8)<jupyter_output><empty_output><jupyter_text>2. Plot amplitude map<jupyter_code>template = eap_cut
locations = hdmea.positions[:, :-1]  # save only x-y positions
fig2, ax2 = plt.subplots(1, 1)
ax2 = plot_amplitude_map(template, locations, ax=ax2, log=True, colorbar=True, colorbar_orientation="horizontal")<jupyter_output><empty_output><jupyter_text>3. Insets with intra and extra potential<jupyter_code>pos1 = [-200, -375]
pos2 = [-311, -132]
pos3 = [-400, 375]
pos4 = [-30, 685]

inset_positions = [pos1, pos2, pos3, pos4]
np.max(eap_cut)
vmem_ylim = [-100, 40]
vext_ylim = [-4, 1]
ms_samples = int(fs)
inset_figs = []
for i, pos in enumerate(inset_positions):
    cell_idx = cell.get_closest_idx(pos[0], pos[1], 0)
    cell_pos = [np.mean(cell.x[cell_idx]), np.mean(cell.y[cell_idx])]
    extra_el_idx = hdmea.get_closest_electrode_idx([cell_pos[0], cell_pos[1], 0])
    
    ax2.plot(np.mean(cell.x[cell_idx]), np.mean(cell.y[cell_idx]), color=f"C{i}", marker='o')
    
    fig_in, axs = plt.subplots(2, 1)
    axs[0].plot(cell.vmem[cell_idx, intra_tidxs], color=f"C{i}", ls='-')
    axs[1].plot(eap_cut[extra_el_idx], color=f"C{i}", ls='--')
    
    axs[0].axvline(min_idx, color='k', ls='--', alpha=0.3)
    axs[1].axvline(min_idx, color='k', ls='--', alpha=0.3)
    
    axs[0].axis('off')
    axs[1].axis('off')
    
    axs[0].set_ylim(vmem_ylim)
    axs[1].set_ylim(vext_ylim)
    
    axs[0].plot([0, 0], [0, 20], f"C{i}")
    axs[1].plot([0, 0], [-2, -1], f"C{i}")
    
    axs[0].text(2, 8, "20mV", color=f"C{i}", fontsize=18)
    axs[1].text(2, -1.8, "1$\mu$V",  color=f"C{i}", fontsize=18)
    
    if i == 0:
        axs[1].plot([min_idx + 5, min_idx + 5 + 2*ms_samples], [-3, -3], lw=2, color="k")
        axs[1].text(min_idx + 8, -4, "2ms",  color=f"k", fontsize=18)
    inset_figs.append(fig_in)<jupyter_output><empty_output><jupyter_text>Save figures<jupyter_code>save_fig = True
if save_fig:
    fig1.savefig(fig_folder / "panelA_mea-neuron.png", dpi=300)
    fig2.savefig(fig_folder / "panelB_amp-map.png", dpi=300)    
    
    for i, inset in enumerate(inset_figs):
        inset.savefig(fig_folder / f"inset-{i}.svg")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13224781572676728, 'fraction_numerical': 0.027535080751919513, 'mean_word_length': 3.883645765998707, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 42, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '138529', 'license': 'MIT', 'n_tokens_mistral': 3004, 'n_tokens_neox': 2782, 'n_words': 621, 'path': 'figure_notebooks/figure2_extracellular-simulations.ipynb', 'repo_name': 'alejoe91/axon_velocity', 'size': '7554'}