Skip to content

Commit 9120b01

Browse files
committed
compare nml and netpy
1 parent aef326e commit 9120b01

File tree

7 files changed

+172
-33
lines changed

7 files changed

+172
-33
lines changed

NeuroML2/compare_MC/RE/RE_netpy.py

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -13,8 +13,8 @@
1313
soma['mechs']['pas'] = {'g': 5e-5, 'e': -77}
1414
soma['mechs']['cadad'] = {'cainf': 0.00024, 'depth': 1, 'kd': 0.0, 'kt': 0.0, 'taur': 5}
1515
# soma['mechs']['kl'] = {'gmax': 3e-06}
16-
soma['mechs']['itre'] = {"gmax": 0.002,"shift": 2.0}
17-
# soma['mechs']['hh2ad'] = {"gkbar": 0.01,"gnabar": 0.09,"vtraub": -50.0}
16+
# soma['mechs']['itre'] = {"gmax": 0.002,"shift": 2.0}
17+
soma['mechs']['hh2ad'] = {"gkbar": 0.01,"gnabar": 0.09,"vtraub": -50.0}
1818

1919

2020
RE_HH_reduced_dict = {'secs': {'soma': soma}}
@@ -25,30 +25,30 @@
2525

2626
netParams.popParams['RE'] = {'cellType': 'RE_HH_reduced', 'numCells': 1}
2727

28-
netParams.stimSourceParams['Input'] = {'type': 'IClamp', 'dur': 400, 'del': 500, 'amp': 0.02}
28+
netParams.stimSourceParams['Input'] = {'type': 'IClamp', 'dur': 500, 'del': 200, 'amp': 0.3}
2929
netParams.stimTargetParams['Input->RE'] = {'source': 'Input', 'sec': 'soma', 'loc': 0.5, 'conds': {'cellType': 'RE_HH_reduced'}}
3030
netParams.defaultThreshold = 5.0
3131
# Simulation options
3232
simConfig = specs.SimConfig() # object of class SimConfig to store simulation configuration
3333

3434
simConfig.recordCells = ['all']
3535
simConfig.hParams['celsius'] = 34
36-
simConfig.duration = 2000 # Duration of the simulation, in ms
36+
simConfig.duration = 1000 # Duration of the simulation, in ms
3737
simConfig.dt = 0.001
3838
# Internal integration timestep to use
3939
simConfig.verbose = True # Show detailed messages
4040

4141

4242
simConfig.recordTraces = {'V_soma':{'sec':'soma','loc':0.5,'var':'v'},
43-
'm_itre':{'sec': 'soma', 'loc': 0.5, 'mech': 'itre', 'var': 'm'},
44-
'h_itre':{'sec': 'soma', 'loc': 0.5, 'mech': 'itre', 'var': 'h'},
45-
'ica_itre':{'sec': 'soma', 'loc': 0.5, 'var': 'ica'},
46-
'carev_itre':{'sec': 'soma', 'loc': 0.5, 'mech': 'itre', 'var': 'carev'},
47-
'caConc_itre':{'sec': 'soma', 'loc': 0.5, 'var': 'cai'}} # Dict with traces to record
43+
'm_hh2_na':{'sec': 'soma', 'loc': 0.5, 'mech': 'hh2ad', 'var': 'm'},
44+
'h_hh2_na':{'sec': 'soma', 'loc': 0.5, 'mech': 'hh2ad', 'var': 'h'},
45+
'n_hh2_k':{'sec': 'soma', 'loc': 0.5, 'mech': 'hh2ad', 'var': 'n'},
46+
'ina':{'sec': 'soma', 'loc': 0.5, 'var': 'ina'},
47+
'ik':{'sec': 'soma', 'loc': 0.5, 'var': 'ik'}} # Dict with traces to record
4848

4949

50-
simConfig.recordStep = 0.01 # Step size in ms to save data (eg. V traces, LFP, etc)
51-
simConfig.filename = 'RE_reduced_itre' # Set file output name
50+
simConfig.recordStep = 0.001 # Step size in ms to save data (eg. V traces, LFP, etc)
51+
simConfig.filename = 'RE_reduced_hh2' # Set file output name
5252
simConfig.savePickle = False # Save params, network and sim output to pickle file
5353
simConfig.saveDataInclude = ['simData']
5454
simConfig.saveJson = True

NeuroML2/compare_MC/RE/RE_reduced_cell.nml

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -29,10 +29,11 @@
2929
<biophysicalProperties id="biophys">
3030
<membraneProperties>
3131

32-
<channelDensityNernst id="itre_soma" ionChannel="itre" condDensity="0.002 S_per_cm2" ion="ca"/>
33-
3432
<channelDensity ionChannel="pas" id="pas_soma" ion="non_specific" erev="-77.0 mV" condDensity="5e-5 S_per_cm2"/>
3533

34+
<channelDensity ionChannel="hh2_k" id="hh2_k_soma" ion="k" erev="-95.0 mV" condDensity="0.01 S_per_cm2"/>
35+
36+
<channelDensity ionChannel="hh2_na" id="hh2_na_soma" ion="na" erev="50.0 mV" condDensity="0.09 S_per_cm2"/>
3637

3738
<spikeThresh value="5mV"/>
3839
<specificCapacitance value="1 uF_per_cm2"/>
1010 KB
Loading
1.01 MB
Loading
Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
1+
import numpy as np
2+
import matplotlib.pyplot as plt
3+
import json
4+
import os
5+
6+
def compare_all():
7+
dat_file = "RE_reduced_cell_step_test.RE_reduced_cell_pop.v.dat"
8+
json_file = "RE_reduced_hh2_data.json"
9+
ina_hh2_dat_file = "hh2_na_iDensity.dat"
10+
ik_hh2_dat_file = "hh2_k_iDensity.dat"
11+
m_hh2_na_dat_file = "m_hh2_na_state.dat"
12+
h_hh2_na_dat_file = "h_hh2_na_state.dat"
13+
n_hh2_k_dat_file = "n_hh2_k_state.dat"
14+
15+
dat_data = np.loadtxt(dat_file)
16+
time_dat = dat_data[:, 0]
17+
voltage_dat = dat_data[:, 1] * 1000
18+
19+
ina_hh2_dat_data = np.loadtxt(ina_hh2_dat_file)
20+
time_ina_hh2_dat = ina_hh2_dat_data[:, 0]
21+
ina_hh2_dat = ina_hh2_dat_data[:, 1] * -0.1
22+
23+
ik_hh2_dat_data = np.loadtxt(ik_hh2_dat_file)
24+
time_ik_hh2_dat = ik_hh2_dat_data[:, 0]
25+
ik_hh2_dat = ik_hh2_dat_data[:, 1] * -0.1
26+
27+
m_hh2_na_dat_data = np.loadtxt(m_hh2_na_dat_file)
28+
time_m_hh2_na_dat = m_hh2_na_dat_data[:, 0]
29+
m_hh2_na_dat = m_hh2_na_dat_data[:, 1]
30+
31+
h_hh2_na_dat_data = np.loadtxt(h_hh2_na_dat_file)
32+
time_h_hh2_na_dat = h_hh2_na_dat_data[:, 0]
33+
h_hh2_na_dat = h_hh2_na_dat_data[:, 1]
34+
35+
n_hh2_k_dat_data = np.loadtxt(n_hh2_k_dat_file)
36+
time_n_hh2_k_dat = n_hh2_k_dat_data[:, 0]
37+
n_hh2_k_dat = n_hh2_k_dat_data[:, 1]
38+
39+
with open(json_file, 'r') as f:
40+
json_data = json.load(f)
41+
voltage_json = json_data['simData']['V_soma']['cell_0']
42+
ina_hh2_json = json_data['simData']['ina']['cell_0']
43+
ik_hh2_json = json_data['simData']['ik']['cell_0']
44+
m_hh2_na_json = json_data['simData']['m_hh2_na']['cell_0']
45+
h_hh2_na_json = json_data['simData']['h_hh2_na']['cell_0']
46+
n_hh2_k_json = json_data['simData']['n_hh2_k']['cell_0']
47+
48+
dt = 1e-6
49+
time_json = np.arange(len(voltage_json)) * dt
50+
51+
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
52+
53+
axes[0, 0].plot(time_dat, voltage_dat, 'b-', label='nml Voltage', linewidth=2)
54+
axes[0, 0].plot(time_json, voltage_json, 'r--', label='netpyne Voltage', linewidth=1.5)
55+
axes[0, 0].set_xlabel('Time (s)')
56+
axes[0, 0].set_ylabel('Voltage (mV)')
57+
axes[0, 0].set_title('Voltage Comparison')
58+
axes[0, 0].legend()
59+
axes[0, 0].grid(True, linestyle='--', alpha=0.7)
60+
61+
axes[0, 1].plot(time_ina_hh2_dat, ina_hh2_dat, 'b-', label='nml ina_hh2', linewidth=2)
62+
axes[0, 1].plot(time_json, ina_hh2_json, 'r--', label='netpyne ina_hh2', linewidth=1.5)
63+
axes[0, 1].set_xlabel('Time (s)')
64+
axes[0, 1].set_ylabel('ina_hh2')
65+
axes[0, 1].set_title('ina_hh2 Comparison')
66+
axes[0, 1].legend()
67+
axes[0, 1].grid(True, linestyle='--', alpha=0.7)
68+
69+
axes[0, 2].plot(time_ik_hh2_dat, ik_hh2_dat, 'b-', label='nml ik_hh2', linewidth=2)
70+
axes[0, 2].plot(time_json, ik_hh2_json, 'r--', label='netpyne ik_hh2', linewidth=1.5)
71+
axes[0, 2].set_xlabel('Time (s)')
72+
axes[0, 2].set_ylabel('ik_hh2')
73+
axes[0, 2].set_title('ik_hh2 Comparison')
74+
axes[0, 2].legend()
75+
axes[0, 2].grid(True, linestyle='--', alpha=0.7)
76+
77+
axes[1, 0].plot(time_m_hh2_na_dat, m_hh2_na_dat, 'b-', label='nml m_hh2_na', linewidth=2)
78+
axes[1, 0].plot(time_json, m_hh2_na_json, 'r--', label='netpyne m_hh2_na', linewidth=1.5)
79+
axes[1, 0].set_xlabel('Time (s)')
80+
axes[1, 0].set_ylabel('m_hh2_na')
81+
axes[1, 0].set_title('m_hh2_na Comparison')
82+
axes[1, 0].legend()
83+
axes[1, 0].grid(True, linestyle='--', alpha=0.7)
84+
85+
axes[1, 1].plot(time_h_hh2_na_dat, h_hh2_na_dat, 'b-', label='nml h_hh2_na', linewidth=2)
86+
axes[1, 1].plot(time_json, h_hh2_na_json, 'r--', label='netpyne h_hh2_na', linewidth=1.5)
87+
axes[1, 1].set_xlabel('Time (s)')
88+
axes[1, 1].set_ylabel('h_hh2_na')
89+
axes[1, 1].set_title('h_hh2_na Comparison')
90+
axes[1, 1].legend()
91+
axes[1, 1].grid(True, linestyle='--', alpha=0.7)
92+
93+
axes[1, 2].plot(time_n_hh2_k_dat, n_hh2_k_dat, 'b-', label='nml n_hh2_k', linewidth=2)
94+
axes[1, 2].plot(time_json, n_hh2_k_json, 'r--', label='netpyne n_hh2_k', linewidth=1.5)
95+
axes[1, 2].set_xlabel('Time (s)')
96+
axes[1, 2].set_ylabel('n_hh2_k')
97+
axes[1, 2].set_title('n_hh2_k Comparison')
98+
axes[1, 2].legend()
99+
axes[1, 2].grid(True, linestyle='--', alpha=0.7)
100+
101+
plt.tight_layout()
102+
current_dir = os.getcwd()
103+
save_path = os.path.join(current_dir, "compare_hh2.png")
104+
plt.savefig(save_path, dpi=300, bbox_inches='tight')
105+
plt.show()
106+
107+
if __name__ == "__main__":
108+
compare_all()

NeuroML2/compare_MC/RE/compare_itre.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@ def compare_all():
3838
ica_itre_json = json_data['simData']['ica_itre']['cell_0']
3939
caConc_json = json_data['simData']['caConc_itre']['cell_0']
4040

41-
dt = 1e-5
41+
dt = 1e-6
4242
time_json = np.arange(len(voltage_json)) * dt
4343

4444
fig, axes = plt.subplots(3, 2, figsize=(15, 15))

NeuroML2/compare_MC/RE/omv_test.py

Lines changed: 49 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -26,14 +26,14 @@
2626
from pyneuroml.io import write_neuroml2_file, read_neuroml2_file
2727
from pyneuroml.lems import generate_lems_file_for_neuroml
2828
from pyneuroml.runners import run_lems_with_jneuroml
29-
from pyneuroml.plot import generate_plot
29+
from pyneuroml.plot import generate_plot
3030
from pyneuroml.neuron.analysis.HHanalyse import get_states
3131
from pyneuroml.analysis.NML2ChannelAnalysis import get_channel_gates
32-
from pyneuroml.utils.plot import get_next_hex_color
32+
from pyneuroml.utils.plot import get_next_hex_color
3333
import neuron
3434
import datetime
3535
import random
36-
from pyneuroml.analysis import generate_current_vs_frequency_curve
36+
from pyneuroml.analysis import generate_current_vs_frequency_curve
3737
random.seed(1412)
3838
def step_current_omv():
3939
netdoc = read_neuroml2_file("RE_reduced_cell.nml")
@@ -43,7 +43,7 @@ def step_current_omv():
4343
pg1 = netdoc.add(
4444
neuroml.PulseGenerator(
4545
id="pg1", delay="200ms", duration="500ms",
46-
amplitude="20pA"
46+
amplitude="300pA"
4747
)
4848
)
4949
input_list1 = net.add(
@@ -59,31 +59,58 @@ def step_current_omv():
5959
)
6060
write_neuroml2_file(netdoc, "RE_reduced_cell.net.nml")
6161
recorder_dict = {}
62-
channel = "itre"
6362

64-
channel_doc = read_neuroml2_file(f"{channel}.channel.nml")
65-
ion_channel_Nernst = channel_doc.ion_channel[0]
66-
gates = get_channel_gates(ion_channel_Nernst)
63+
channel_doc1 = read_neuroml2_file(f"hh2_na.channel.nml")
64+
ion_channel1 = channel_doc1.ion_channel[0]
65+
66+
channel_doc2 = read_neuroml2_file(f"hh2_k.channel.nml")
67+
ion_channel2 = channel_doc2.ion_channel[0]
68+
69+
'''
70+
71+
recorder_dict[f"{channel}_erev.dat"] = [
72+
f"{pop.id}[0]/biophys/membraneProperties/itre_soma/erev"
73+
]
6774
68-
for gate in gates:
69-
recorder_dict[f"{gate}_{channel}_state.dat"] = [
70-
f"{pop.id}[0]/biophys/membraneProperties/itre_soma/{channel}/{gate}/q"
71-
]
75+
recorder_dict["caConcExt.dat"] = [
76+
f"{pop.id}[0]/caConcExt"
77+
]
78+
79+
80+
recorder_dict["caConc.dat"] = [
81+
f"{pop.id}[0]/caConc"
82+
]
7283
84+
85+
'''
86+
87+
'''
7388
recorder_dict[f"{channel}_iDensity.dat"] = [
7489
f"{pop.id}[0]/biophys/membraneProperties/itre_soma/iDensity"
7590
]
91+
'''
92+
93+
gates1 = get_channel_gates(ion_channel1)
7694

77-
recorder_dict[f"{channel}_erev.dat"] = [
78-
f"{pop.id}[0]/biophys/membraneProperties/itre_soma/erev"
79-
]
95+
for gate in gates1:
96+
recorder_dict[f"{gate}_hh2_na_state.dat"] = [
97+
f"{pop.id}[0]/biophys/membraneProperties/hh2_na_soma/hh2_na/{gate}/q"
98+
]
99+
100+
gates2 = get_channel_gates(ion_channel2)
80101

81-
recorder_dict["caConc.dat"] = [
82-
f"{pop.id}[0]/caConc"
102+
for gate in gates2:
103+
recorder_dict[f"{gate}_hh2_k_state.dat"] = [
104+
f"{pop.id}[0]/biophys/membraneProperties/hh2_k_soma/hh2_k/{gate}/q"
105+
]
106+
107+
108+
recorder_dict[f"hh2_k_iDensity.dat"] = [
109+
f"{pop.id}[0]/biophys/membraneProperties/hh2_k_soma/iDensity"
83110
]
84111

85-
recorder_dict["caConcExt.dat"] = [
86-
f"{pop.id}[0]/caConcExt"
112+
recorder_dict[f"hh2_na_iDensity.dat"] = [
113+
f"{pop.id}[0]/biophys/membraneProperties/hh2_na_soma/iDensity"
87114
]
88115

89116
generate_lems_file_for_neuroml(
@@ -103,6 +130,8 @@ def step_current_omv():
103130
"LEMS_RE_reduced_cell_step_test.xml", load_saved_data=True, nogui=True
104131
)
105132
print(data.keys())
133+
134+
'''
106135
generate_plot(
107136
xvalues=[data["t"]],
108137
yvalues=[data["RE_reduced_cell_pop[0]/v"]],
@@ -188,6 +217,7 @@ def step_current_omv():
188217
title_above_plot="NML",
189218
legend_position="outer right"
190219
)
220+
'''
191221

192222
if __name__ == "__main__":
193223
step_current_omv()

0 commit comments

Comments
 (0)