Impedance spectrum in cortical tissue: implications for LFP signal propagation (Miceli et al. 2017)

 Download zip file 
Help downloading and running models
Accession:224923
" ... Here, we performed a detailed investigation of the frequency dependence of the conductivity within cortical tissue at microscopic distances using small current amplitudes within the typical (neuro)physiological micrometer and sub-nanoampere range. We investigated the propagation of LFPs, induced by extracellular electrical current injections via patch-pipettes, in acute rat brain slice preparations containing the somatosensory cortex in vitro using multielectrode arrays. Based on our data, we determined the cortical tissue conductivity over a 100-fold increase in signal frequency (5-500 Hz). Our results imply at most very weak frequency-dependent effects within the frequency range of physiological LFPs. Using biophysical modeling, we estimated the impact of different putative impedance spectra. Our results indicate that frequency dependencies of the order measured here and in most other studies have negligible impact on the typical analysis and modeling of LFP signals from extracellular brain recordings."
Reference:
1 . Miceli S, Ness TV, Einevoll GT, Schubert D (2017) Impedance Spectrum in Cortical Tissue: Implications for Propagation of LFP Signals on the Microscopic Level Eneuro 4:1-15
Model Information (Click on a link to find other models with that property)
Model Type: Extracellular; Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Neocortex V1 L6 pyramidal corticothalamic GLU cell;
Channel(s): I Na,p; I Na,t; I L high threshold; I A; I h; I K,Ca; I Calcium; I A, slow;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: Python; NEURON;
Model Concept(s): Extracellular Fields; Methods; Simplified Models; Detailed Neuronal Models;
Implementer(s):
Search NeuronDB for information about:  Neocortex V1 L6 pyramidal corticothalamic GLU cell; I Na,p; I Na,t; I L high threshold; I A; I h; I K,Ca; I Calcium; I A, slow;
/
tissue_impedance_impact
hay
sim_results
ampa.mod
Ca_HVA.mod *
Ca_LVAst.mod *
CaDynamics_E2.mod *
epsp.mod *
glutamate.mod
Ih.mod *
Im.mod *
K_Pst.mod *
K_Tst.mod *
Nap_Et2.mod
NaTa_t.mod *
NaTs2_t.mod
siniclamp.mod
SK_E2.mod *
SKv3_1.mod *
stim.mod
cell1.hoc
cell1.rot
custom_codes.hoc
hay_active_declarations.py
                            
__author__ = 'torbjone'
#!/usr/bin/env python

import os
from os.path import join
import sys
import numpy as np
import pylab as plt
import neuron
nrn = neuron.h
import LFPy

def make_cell_uniform(Vrest=-80):
    neuron.h.t = 0
    neuron.h.finitialize(Vrest)
    neuron.h.fcurrent()
    for sec in neuron.h.allsec():
        for seg in sec:
            seg.e_pas = seg.v
            if neuron.h.ismembrane("na_ion"):
                seg.e_pas += seg.ina/seg.g_pas
            if neuron.h.ismembrane("k_ion"):
                seg.e_pas += seg.ik/seg.g_pas
            if neuron.h.ismembrane("ca_ion"):
                seg.e_pas = seg.e_pas + seg.ica/seg.g_pas
            if neuron.h.ismembrane("Ih"):
                seg.e_pas += seg.ihcn_Ih/seg.g_pas
            if neuron.h.ismembrane("Ih_z"):
                seg.e_pas += seg.ih_Ih_z/seg.g_pas
            if neuron.h.ismembrane("Ih_frozen"):
                seg.e_pas += seg.ihcn_Ih_frozen/seg.g_pas
            if neuron.h.ismembrane("Ih_linearized_mod"):
                seg.e_pas = seg.e_pas + seg.ihcn_Ih_linearized_mod/seg.g_pas
            if neuron.h.ismembrane("Ih_linearized_v2"):
                seg.e_pas = seg.e_pas + seg.ihcn_Ih_linearized_v2/seg.g_pas
            if neuron.h.ismembrane("Ih_linearized_v2_frozen"):
                seg.e_pas = seg.e_pas + seg.ihcn_Ih_linearized_v2_frozen/seg.g_pas

def _get_longest_distance():

    nrn.distance(0, 0.5)
    max_dist = 0
    for sec in nrn.allsec():
        for seg in sec:
            max_dist = np.max([max_dist, nrn.distance(seg.x)])
    return max_dist

def _get_total_area():

    total_area = 0
    for sec in nrn.allsec():
        for seg in sec:
           # Never mind the units, as long as it is consistent
           total_area += nrn.area(seg.x)
    return total_area

def _get_linear_increase_factor(increase_factor, max_dist, total_conductance):
    normalization = 0
    for sec in neuron.h.allsec():
        for seg in sec:
            normalization += nrn.area(seg.x) * (1 + (increase_factor - 1) * nrn.distance(seg.x)/max_dist)
    return total_conductance / normalization

def _get_linear_decrease_factor(decrease_factor, max_dist, total_conductance):
    normalization = 0
    for sec in neuron.h.allsec():
        for seg in sec:
            normalization += nrn.area(seg.x) * (decrease_factor - (decrease_factor - 1) *
                                                nrn.distance(seg.x)/max_dist)
    return total_conductance / normalization


def biophys_zuchkova(**kwargs):
    nrn.distance(0, 0.5)
    for sec in neuron.h.allsec():
        sec.insert('pas')
        sec.e_pas = kwargs['hold_potential']
        sec.Ra = 200
        sec.cm = 1.0
        sec.g_pas = 0.09e-3
        sec.insert("Ih_z")
        for seg in sec:
            seg.gIhbar_Ih_z = 26.71e-6 * np.exp(0.0041*nrn.distance(seg.x))

    if 'hold_potential' in kwargs:
        make_cell_uniform(Vrest=kwargs['hold_potential'])

def biophys_generic(**kwargs):

    for sec in neuron.h.allsec():
        sec.insert("QA")
        sec.V_r_QA = kwargs['hold_potential']
        sec.tau_w_QA = kwargs['tau_w']
        sec.Ra = 100
        sec.cm = 1.0
        sec.g_pas_QA = kwargs['g_pas']

    total_w_conductance = kwargs['total_w_conductance']
    total_area = _get_total_area()
    max_dist = _get_longest_distance()

    if kwargs['distribution'] == 'uniform':
        for sec in neuron.h.allsec():
            for seg in sec:
                seg.g_w_QA = total_w_conductance / total_area

    elif kwargs['distribution'] == 'linear_increase':
        increase_factor = 100
        conductance_factor = _get_linear_increase_factor(increase_factor, max_dist, total_w_conductance)

        nrn.distance(0, 0.5)
        for sec in neuron.h.allsec():
            for seg in sec:
                seg.g_w_QA = conductance_factor * (1 + (increase_factor - 1) * nrn.distance(seg.x) / max_dist)

    elif kwargs['distribution'] == 'linear_decrease':
        decrease_factor = 100
        conductance_factor = _get_linear_decrease_factor(decrease_factor, max_dist, total_w_conductance)
        nrn.distance(0, 0.5)
        for sec in neuron.h.allsec():
            for seg in sec:
                seg.g_w_QA = (conductance_factor * (decrease_factor - (decrease_factor - 1) *
                                                    nrn.distance(seg.x) / max_dist))
    else:
        raise RuntimeError("Unknown distribution...")

    for sec in neuron.h.allsec():
        for seg in sec:
            seg.mu_QA = seg.g_w_QA / seg.g_pas_QA * kwargs['mu_factor']


def biophys_passive(**kwargs):
    Vrest = kwargs['hold_potential'] if 'hold_potential' in kwargs else -70
    for sec in neuron.h.allsec():
        sec.insert('pas')
        sec.cm = 1.0
        sec.Ra = 100.
        sec.e_pas = Vrest

    for sec in neuron.h.soma:
        sec.g_pas = 0.0000338

    for sec in neuron.h.apic:
        sec.cm = 2
        sec.g_pas = 0.0000589

    for sec in neuron.h.dend:
        sec.cm = 2
        sec.g_pas = 0.0000467

    for sec in neuron.h.axon:
        sec.g_pas = 0.0000325

    if 'hold_potential' in kwargs:
        make_cell_uniform(Vrest=kwargs['hold_potential'])

    print("Passive dynamics inserted.")


def biophys_Ih_linearized_frozen(**kwargs):

    Vrest = kwargs['hold_potential'] if 'hold_potential' in kwargs else -70

    for sec in neuron.h.allsec():
        sec.insert('pas')
        sec.cm = 1.0
        sec.Ra = 100.
        sec.e_pas = Vrest
    for sec in neuron.h.soma:
        sec.insert("Ih_linearized_v2_frozen")
        sec.gIhbar_Ih_linearized_v2_frozen = 0.0002
        sec.g_pas = 0.0000338
        sec.V_R_Ih_linearized_v2_frozen = Vrest
    for sec in neuron.h.apic:
        sec.insert("Ih_linearized_v2_frozen")
        sec.cm = 2
        sec.g_pas = 0.0000589
        sec.V_R_Ih_linearized_v2_frozen = Vrest

    nrn.distribute_channels("apic", "gIhbar_Ih_linearized_v2_frozen",
                            2, -0.8696, 3.6161, 0.0, 2.0870, 0.0002)
    for sec in neuron.h.dend:
        sec.insert("Ih_linearized_v2_frozen")
        sec.cm = 2
        sec.g_pas = 0.0000467
        sec.gIhbar_Ih_linearized_v2_frozen = 0.0002
        sec.V_R_Ih_linearized_v2_frozen = Vrest
    for sec in neuron.h.axon:
        sec.g_pas = 0.0000325
    if 'hold_potential' in kwargs:
        make_cell_uniform(Vrest=kwargs['hold_potential'])

    print("Frozen linearized Ih inserted.")

def biophys_Ih_linearized(**kwargs):

    Vrest = kwargs['hold_potential'] if 'hold_potential' in kwargs else -70

    for sec in neuron.h.allsec():
        sec.insert('pas')
        sec.cm = 1.0
        sec.Ra = 100.
        sec.e_pas = Vrest
    for sec in neuron.h.soma:
        sec.insert("Ih_linearized_v2")
        sec.gIhbar_Ih_linearized_v2 = 0.0002
        sec.g_pas = 0.0000338
        # sec.vss_Ih_linearized_mod = Vrest
        sec.V_R_Ih_linearized_v2 = Vrest
    for sec in neuron.h.apic:
        sec.insert("Ih_linearized_v2")
        sec.cm = 2
        sec.g_pas = 0.0000589
        # sec.vss_Ih_linearized_mod = Vrest
        sec.V_R_Ih_linearized_v2 = Vrest

    nrn.distribute_channels("apic", "gIhbar_Ih_linearized_v2",
                            2, -0.8696, 3.6161, 0.0, 2.0870, 0.00020000000)
    for sec in neuron.h.dend:
        sec.insert("Ih_linearized_v2")
        sec.cm = 2
        sec.g_pas = 0.0000467
        sec.gIhbar_Ih_linearized_v2 = 0.0002
        # sec.vss_Ih_linearized_v2 = Vrest
        sec.V_R_Ih_linearized_v2 = Vrest
    for sec in neuron.h.axon:
        sec.g_pas = 0.0000325
    if 'hold_potential' in kwargs:
        make_cell_uniform(Vrest=kwargs['hold_potential'])

    print("Linearized Ih inserted.")


def biophys_active(**kwargs):

    for sec in neuron.h.allsec():
        sec.insert('pas')
        sec.cm = 1.0
        sec.Ra = 100.
        sec.e_pas = -90.

    for sec in neuron.h.soma:
        sec.insert('Ca_LVAst')
        sec.insert('Ca_HVA')
        sec.insert('SKv3_1')
        sec.insert('SK_E2')
        sec.insert('K_Tst')
        sec.insert('K_Pst')
        sec.insert('Nap_Et2')
        sec.insert('NaTa_t')
        sec.insert('CaDynamics_E2')
        sec.insert('Ih')
        sec.ek = -85
        sec.ena = 50
        sec.gIhbar_Ih = 0.0002
        sec.g_pas = 0.0000338
        sec.decay_CaDynamics_E2 = 460.0
        sec.gamma_CaDynamics_E2 = 0.000501
        sec.gCa_LVAstbar_Ca_LVAst = 0.00343
        sec.gCa_HVAbar_Ca_HVA = 0.000992
        sec.gSKv3_1bar_SKv3_1 = 0.693
        sec.gSK_E2bar_SK_E2 = 0.0441
        sec.gK_Tstbar_K_Tst = 0.0812
        sec.gK_Pstbar_K_Pst = 0.00223
        sec.gNap_Et2bar_Nap_Et2 = 0.00172
        sec.gNaTa_tbar_NaTa_t = 2.04

    for sec in neuron.h.apic:
        sec.cm = 2
        sec.insert('Ih')
        sec.insert('SK_E2')
        sec.insert('Ca_LVAst')
        sec.insert('Ca_HVA')
        sec.insert('SKv3_1')
        sec.insert('NaTa_t')
        sec.insert('Im')
        sec.insert('CaDynamics_E2')
        sec.ek = -85
        sec.ena = 50
        sec.decay_CaDynamics_E2 = 122
        sec.gamma_CaDynamics_E2 = 0.000509
        sec.gSK_E2bar_SK_E2 = 0.0012
        sec.gSKv3_1bar_SKv3_1 = 0.000261
        sec.gNaTa_tbar_NaTa_t = 0.0213
        sec.gImbar_Im = 0.0000675
        sec.g_pas = 0.0000589

    nrn.distribute_channels("apic", "gIhbar_Ih", 2, -0.8696, 3.6161, 0.0, 2.087, 0.0002)
    nrn.distribute_channels("apic", "gCa_LVAstbar_Ca_LVAst", 3, 1.0, 0.010, 685.0, 885.0, 0.0187)
    nrn.distribute_channels("apic", "gCa_HVAbar_Ca_HVA", 3, 1.0, 0.10, 685.00, 885.0, 0.000555)

    for sec in neuron.h.dend:
        sec.cm = 2
        sec.insert('Ih')
        sec.gIhbar_Ih = 0.0002
        sec.g_pas = 0.0000467

    for sec in neuron.h.axon:
        sec.g_pas = 0.0000325

    if 'hold_potential' in kwargs:
        make_cell_uniform(Vrest=kwargs['hold_potential'])
    print("active ion-channels inserted.")


def biophys_active_frozen(**kwargs):

    for sec in neuron.h.allsec():
        sec.insert('pas')
        sec.cm = 1.0
        sec.Ra = 100.
        sec.e_pas = -90.

    for sec in neuron.h.soma:
        sec.insert('Ca_LVAst_frozen')
        sec.insert('Ca_HVA_frozen')
        sec.insert('SKv3_1_frozen')
        sec.insert('SK_E2_frozen')
        sec.insert('K_Tst_frozen')
        sec.insert('K_Pst_frozen')
        sec.insert('Nap_Et2_frozen')
        sec.insert('NaTa_t_frozen')
        sec.insert('CaDynamics_E2')
        sec.insert('Ih_frozen')
        sec.ek = -85
        sec.ena = 50
        sec.gIhbar_Ih_frozen = 0.0002
        sec.g_pas = 0.0000338
        sec.decay_CaDynamics_E2 = 460.0
        sec.gamma_CaDynamics_E2 = 0.000501
        sec.gCa_LVAstbar_Ca_LVAst_frozen = 0.00343
        sec.gCa_HVAbar_Ca_HVA_frozen = 0.000992
        sec.gSKv3_1bar_SKv3_1_frozen = 0.693
        sec.gSK_E2bar_SK_E2_frozen = 0.0441
        sec.gK_Tstbar_K_Tst_frozen = 0.0812
        sec.gK_Pstbar_K_Pst_frozen = 0.00223
        sec.gNap_Et2bar_Nap_Et2_frozen = 0.00172
        sec.gNaTa_tbar_NaTa_t_frozen = 2.04

    for sec in neuron.h.apic:
        sec.cm = 2
        sec.insert('Ih_frozen')
        sec.insert('SK_E2_frozen')
        sec.insert('Ca_LVAst_frozen')
        sec.insert('Ca_HVA_frozen')
        sec.insert('SKv3_1_frozen')
        sec.insert('NaTa_t_frozen')
        sec.insert('Im_frozen')
        sec.insert('CaDynamics_E2')
        sec.ek = -85
        sec.ena = 50
        sec.decay_CaDynamics_E2 = 122
        sec.gamma_CaDynamics_E2 = 0.000509
        sec.gSK_E2bar_SK_E2_frozen = 0.0012
        sec.gSKv3_1bar_SKv3_1_frozen = 0.000261
        sec.gNaTa_tbar_NaTa_t_frozen = 0.0213
        sec.gImbar_Im_frozen = 0.0000675
        sec.g_pas = 0.0000589

    nrn.distribute_channels("apic", "gIhbar_Ih_frozen", 2, -0.8696, 3.6161, 0.0, 2.0870, 0.00020000000)
    nrn.distribute_channels("apic", "gCa_LVAstbar_Ca_LVAst_frozen", 3, 1.000000, 0.010000, 685.000000, 885.000000, 0.0187000000)
    nrn.distribute_channels("apic", "gCa_HVAbar_Ca_HVA_frozen", 3, 1.000000, 0.100000, 685.000000, 885.000000, 0.0005550000)

    for sec in neuron.h.dend:
        sec.cm = 2
        sec.insert('Ih_frozen')
        sec.gIhbar_Ih_frozen = 0.0002
        sec.g_pas = 0.0000467

    for sec in neuron.h.axon:
        sec.g_pas = 0.0000325

    if 'hold_potential' in kwargs:
        make_cell_uniform(Vrest=kwargs['hold_potential'])
    print("Frozen active ion-channels inserted.")


def make_syaptic_stimuli(cell, input_idx):
    # Define synapse parameters
    synapse_parameters = {
        'idx': input_idx,
        'e': 0.,                   # reversal potential
        'syntype': 'ExpSyn',       # synapse type
        'tau': 10.,                # syn. time constant
        'weight': 0.001,            # syn. weight
        'record_current': True,
    }
    synapse = LFPy.Synapse(cell, **synapse_parameters)
    synapse.set_spike_times(np.array([5.]))
    return cell, synapse

def active_declarations(**kwargs):
    ''' set active conductances for Hay model 2011 '''
    nrn.delete_axon()
    nrn.geom_nseg()
    nrn.define_shape()
    exec('biophys_%s(**kwargs)' % kwargs['conductance_type'])


def simulate_synaptic_input(input_idx, holding_potential, conductance_type):

    timeres = 2**-4
    cut_off = 0
    tstopms = 150
    tstartms = -cut_off
    model_path = join('lfpy_version')
    neuron.load_mechanisms('mod')
    neuron.load_mechanisms('..')
    cell_params = {
        'morphology': join(model_path, 'morphologies', 'cell1.hoc'),
        #'rm' : 30000,               # membrane resistance
        #'cm' : 1.0,                 # membrane capacitance
        #'Ra' : 100,                 # axial resistance
        'v_init': holding_potential,             # initial crossmembrane potential
        'passive': False,           # switch on passive mechs
        'nsegs_method': 'lambda_f',  # method for setting number of segments,
        'lambda_f': 100,           # segments are isopotential at this frequency
        'timeres_NEURON': timeres,   # dt of LFP and NEURON simulation.
        'timeres_python': 1,
        'tstartms': tstartms,          # start time, recorders start at t=0
        'tstopms': tstopms,
        'custom_code': [join(model_path, 'custom_codes.hoc')],
        'custom_fun': [active_declarations],  # will execute this function
        'custom_fun_args': [{'conductance_type': conductance_type,
                             'hold_potential': holding_potential}],
    }

    cell = LFPy.Cell(**cell_params)
    make_syaptic_stimuli(cell, input_idx)
    cell.simulate(rec_vmem=True, rec_imem=True)

    plt.subplot(211, title='Soma')
    plt.plot(cell.tvec, cell.vmem[0, :], label='%d %s %d mV' % (input_idx, conductance_type,
                                                                holding_potential))

    plt.subplot(212, title='Input idx %d' % input_idx)
    plt.plot(cell.tvec, cell.vmem[input_idx, :], label='%d %s %d mV' % (input_idx, conductance_type,
                                                           holding_potential))

def test_frozen_currents(input_idx, holding_potential):

    # input_idx = 0
    # holding_potential = -70
    plt.close('all')
    simulate_synaptic_input(input_idx, holding_potential, 'active_frozen')
    simulate_synaptic_input(input_idx, holding_potential, 'active')
    simulate_synaptic_input(input_idx, holding_potential, 'passive')
    simulate_synaptic_input(input_idx, holding_potential, 'Ih_linearized')
    simulate_synaptic_input(input_idx, holding_potential, 'Ih_linearized_frozen')

    plt.legend(frameon=False)
    plt.savefig('frozen_test_%d_%d.png' % (input_idx, holding_potential))


def test_steady_state():

    timeres = 2**-4
    cut_off = 0
    tstopms = 100
    tstartms = -cut_off
    model_path = join('lfpy_version')

    cell_params = {
        'morphology': join(model_path, 'morphologies', 'cell1.hoc'),
        #'rm' : 30000,               # membrane resistance
        #'cm' : 1.0,                 # membrane capacitance
        #'Ra' : 100,                 # axial resistance
        'v_init': -80,             # initial crossmembrane potential
        'passive': False,           # switch on passive mechs
        'nsegs_method': 'lambda_f',  # method for setting number of segments,
        'lambda_f': 100,           # segments are isopotential at this frequency
        'timeres_NEURON': timeres,   # dt of LFP and NEURON simulation.
        'timeres_python': 1,
        'tstartms': tstartms,          # start time, recorders start at t=0
        'tstopms': tstopms,
        'custom_code': [join(model_path, 'custom_codes.hoc')],
        'custom_fun': [active_declarations],  # will execute this function
        'custom_fun_args': [{'conductance_type': 'passive',
                             'hold_potential': -80}],
    }

    cell = LFPy.Cell(**cell_params)

    cell.simulate(rec_vmem=True, rec_imem=True)

    plt.plot(cell.tvec, cell.somav)
    plt.show()

def return_frozen_gh(holding_potential):
    timeres = 2**-4
    cut_off = 0
    tstopms = 150
    tstartms = -cut_off
    model_path = join('lfpy_version')
    neuron.load_mechanisms('mod')
    neuron.load_mechanisms('..')
    cell_params = {
        'morphology': join(model_path, 'morphologies', 'cell1.hoc'),
        #'rm' : 30000,               # membrane resistance
        #'cm' : 1.0,                 # membrane capacitance
        #'Ra' : 100,                 # axial resistance
        'v_init': holding_potential,             # initial crossmembrane potential
        'passive': False,           # switch on passive mechs
        'nsegs_method': 'lambda_f',  # method for setting number of segments,
        'lambda_f': 100,           # segments are isopotential at this frequency
        'timeres_NEURON': timeres,   # dt of LFP and NEURON simulation.
        'timeres_python': 1,
        'tstartms': tstartms,          # start time, recorders start at t=0
        'tstopms': tstopms,
        'custom_code': [join(model_path, 'custom_codes.hoc')],
        'custom_fun': [active_declarations],  # will execute this function
        'custom_fun_args': [{'conductance_type': 'active_frozen',
                             'hold_potential': holding_potential}],
    }

    cell = LFPy.Cell(**cell_params)
    # make_syaptic_stimuli(cell, input_idx)
    # cell.simulate(rec_vmem=True, rec_imem=True)

    gIh = np.zeros(cell.totnsegs)
    gpas = np.zeros(cell.totnsegs)
    dist = np.zeros(cell.totnsegs)
    nrn.distance()
    comp = 0
    for sec in cell.allseclist:
        for seg in sec:
            gIh[comp] = seg.gIh_Ih_frozen if hasattr(seg, "ihcn_Ih_frozen") else 0
            gpas[comp] = seg.g_pas
            dist[comp] = nrn.distance(seg.x)
            comp += 1

    return dist, gIh, gpas


def plot_frozen_gh():
    dist, gIh_80, gpas = return_frozen_gh(-80)
    dist, gIh_70, gpas = return_frozen_gh(-70)
    dist, gIh_60, gpas = return_frozen_gh(-60)


    plt.plot(dist, gIh_80, 'bo', label='Frozen Ih at -80 mV')
    plt.plot(dist, gIh_70, 'go', label='Frozen Ih at -70 mV')
    plt.plot(dist, gIh_60, 'ro', label='Frozen Ih at -60 mV')
    plt.plot(dist, gpas, 'ko', label='g_pas')

    plt.legend(loc=2)
    plt.ylabel('Conductance [S/cm^2]')
    plt.xlabel('Distance from soma [$\mu m$]')
    plt.savefig('frozen_Ih_conductance.png')

if __name__ == '__main__':

    plot_frozen_gh()
    # test_frozen_currents(0, -80)
    # test_frozen_currents(0, -70)
    # test_frozen_currents(0, -60)
    # test_frozen_currents(750, -60)
    # test_frozen_currents(750, -70)
    # test_frozen_currents(750, -80)

Loading data, please wait...