get_steady_state_Vm.py

import os.path as op
import hnn_core
from hnn_core import read_params, Network
import matplotlib.pyplot as plt
from hnn_core.network_builder import NetworkBuilder
from hnn_core.network_builder import _simulate_single_trial
from neuron import h
import numpy as np
from scipy.optimize import curve_fit


def fitfun(x, a, b, c):
    return a * np.exp(-b * x) + c


hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'default.json')
params = read_params(params_fname)

params.update({'N_pyr_x': 3,
               'N_pyr_y': 3,
               'tstop': 5000.})
n_pyrs = 9
net = Network(params)
net._instantiate_drives(n_trials=1)
L5_pyr_gid = net.gid_ranges['L5_pyramidal'][0]
L2_pyr_gid = net.gid_ranges['L2_pyramidal'][0]

neuron_net = NetworkBuilder(net, trial_idx=0)
dpl = _simulate_single_trial(neuron_net, 0)

L5_steady_state = dpl.data['L5'][-1] / n_pyrs
L2_steady_state = dpl.data['L2'][-1] / n_pyrs

fig, axs = plt.subplots(1, 2)
axs[0].plot(dpl.times, dpl.data['L5'] / n_pyrs)
axs[0].axhline(L5_steady_state, linestyle='--', color='red',
               label=f'{L5_steady_state:.5e}')
axs[0].legend()
axs[1].plot(dpl.times, dpl.data['L2'] / n_pyrs)
axs[1].axhline(L2_steady_state, linestyle='--', color='red',
               label=f'{L2_steady_state:.5e}')
axs[1].legend()

# plot the L5 initial segment
tinds = np.where(dpl.times < 100.)
L5_upclose = dpl.data['L5'][tinds]
plt.plot(dpl.times[tinds], L5_upclose)

# fit form 20 ms
p0 = (-0.3e3 / n_pyrs, 2e-3, L5_steady_state / n_pyrs)
tinds = np.where(dpl.times > 20)
ydata = dpl.data['L5'][tinds] / n_pyrs
popt, pcov = curve_fit(fitfun, dpl.times[tinds], ydata, p0=p0)
legstr = 'fit: {:.2e} x \n exp( {:.2e} t ) + {:.2e}'.format(*popt)

plt.figure()
plt.plot(dpl.times, dpl.data['L5'] / n_pyrs)
plt.plot(dpl.times[tinds], fitfun(dpl.times[tinds], *popt),
         linestyle='--', label=legstr)
plt.xlabel('Time (ms)')
plt.legend()
plt.ticklabel_format(style='sci', axis='y', scilimits=(-2, 2))
print('Fit params (L5): {:.5e} * np.exp({:.5e} * t) + {:.5e}'.format(*popt))
print(f'Fit params (L2): {L2_steady_state:.5e}')
# Fit params (L5): -4.60641e+00 * np.exp(2.51721e-03 * t) + -4.80496e+01
# Fit params (L2): 4.43005e-02

# This should be run after _simulate_single_trial

# get steady-state values for L5 pyramidals
seclist = h.SectionList()
seclist.wholetree(sec=neuron_net.cells[L5_pyr_gid].soma)
for sect in seclist:
    print(f"'{sect.name()}': [", end='')
    for seg in sect:
        print(f'{seg.v},')
    print('],')

# '<L5Pyr | soma: L 39.000000, diam 28.900000, Ra 200.000000, cm 0.850000>.L5Pyr_soma': [-71.99659698854643,
# ],
# 'L5Pyr_basal_1': [-72.02954647612275,
# ],
# 'L5Pyr_apical_trunk': [-71.99291957918324,
# -71.98551742625942,
# -71.97715674520684,
# ],
# 'L5Pyr_basal_3': [-72.06262357371888,
# -72.07093175190667,
# -72.07699267334765,
# -72.08091331425713,
# -72.08280972750579,
# ],
# 'L5Pyr_basal_2': [-72.06262357371888,
# -72.07093175190667,
# -72.07699267334765,
# -72.08091331425713,
# -72.08280972750579,
# ],
# 'L5Pyr_apical_oblique': [-71.98173167226447,
# -71.99604815550023,
# -72.00643562257369,
# -72.01311793566057,
# -72.01633184996282,
# ],
# 'L5Pyr_apical_1': [-71.95478360746081,
# -71.91629048542289,
# -71.87472632967865,
# -71.83008375033546,
# -71.78236437694262,
# -71.73158238928728,
# -71.67776729748056,
# -71.62096774371555,
# -71.56125603324949,
# -71.49873368682518,
# -71.43353829740335,
# -71.36584997180377,
# -71.29589980384837,
# ],
# 'L5Pyr_apical_2': [-71.17714911653867,
# -71.00881744257748,
# -70.8377556814502,
# -70.6641370708994,
# -70.48824928295832,
# -70.3105072818248,
# -70.13147531257948,
# -69.95186936752486,
# -69.77262211942187,
# -69.5949770788285,
# -69.42047576454539,
# -69.25102045904973,
# -69.08894300326746,
# ],
# 'L5Pyr_apical_tuft': [-68.86888820906418,
# -68.59550615150705,
# -68.3391034373857,
# -68.10239662633921,
# -67.88885381988956,
# -67.70294104630771,
# -67.55040156735672,
# -67.43851485243435,
# -67.37658802902182,
# ],

# Layer 2
seclist = h.SectionList()
seclist.wholetree(sec=neuron_net.cells[L2_pyr_gid].soma)
for sect in seclist:
    print(f"'{sect.name()}': [", end='')
    for seg in sect:
        print(f'{seg.v},')
    print('],')

# '<L2Pyr | soma: L 22.100000, diam 23.400000, Ra 200.000000, cm 0.619500>.L2Pyr_soma': [-71.46251994264288,
# ],
# 'L2Pyr_basal_1': [-71.46376735038517,
# ],
# 'L2Pyr_apical_trunk': [-71.46414752119715,
# ],
# 'L2Pyr_basal_3': [-71.46542691374754,
# -71.46653303307605,
# -71.46735309117791,
# -71.46789528516324,
# -71.4681650341793,
# ],
# 'L2Pyr_basal_2': [-71.46542691374754,
# -71.46653303307605,
# -71.46735309117791,
# -71.46789528516324,
# -71.4681650341793,
# ],
# 'L2Pyr_apical_oblique': [-71.46621193907679,
# -71.46721246231539,
# -71.46803677278407,
# -71.46869006901264,
# -71.46917647061586,
# -71.46949904455353,
# -71.46965982467452,
# ],
# 'L2Pyr_apical_1': [-71.46630566733937,
# -71.46753494562554,
# -71.46862901749326,
# -71.4695932385707,
# -71.47043232819412,
# -71.47115039275037,
# -71.47175094597775,
# ],
# 'L2Pyr_apical_tuft': [-71.4723749447141,
# -71.4729782895952,
# -71.47342716185227,
# -71.47372468566407,
# -71.47387293171057,
# ],