Modify the model and add more examples

examples/nest_gpu/Ca-AdEx_Example4.py

@@ -0,0 +1,243 @@
+#
+#  Ca-AdEx: example of usage of two-compartment neuron model
+#
+#  First version: 22/10/2025
+#  Author: Elena Pastorelli, INFN, Rome (IT)
+#
+#  Description: Example of usage of the Ca-AdEx neuron stimulated with spike generators
+#
+#  Copyright © 2025   Elena Pastorelli          <elena.pastorelli@roma1.infn.it>
+#  Copyright © 2025   Pier Stanislao Paolucci   <pier.paolucci@roma1.infn.it>
+#
+#  SPDX-License-Identifier: GPL-3.0-only
+#
+#  This program is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import matplotlib.pyplot as plt
+import numpy as np
+import nestgpu as ngpu
+
+default_param = {'C_m_d': 23.67372778891213, 'C_m_s': 246.7882968598874, 'Ca_0': 0.0001, 'Ca_th': 0.00043,
+                 'V_reset': -61.73952230767877, 'd_BAP': 0.1195980511869619, 'delta_T': 2.0, 'e_K': -90.0,
+                 'e_L_d': -55.000000000000014, 'e_L_s': -69.24596493128396, 'e_Na_Adex': -50.0, 'exp_K_Ca': 4.8,
+                 'g_C_d': 19.777320239615996, 'g_L_d': 3.377855016658499, 'g_L_s': 5.0, 'g_w': 1.1156385639067352,
+                 'gbar_Ca': 21.045506331690845, 'gbar_K_Ca': 13.199867205029523, 'h_half': -21.0, 'h_slope': -0.5,
+                 'm_half': -9.0, 'm_slope': 0.5, 'phi': 3.92830985228413e-08, 't_ref': 0.0,
+                 'tau_decay_Ca': 103.57233790866408, 'tau_h': 80.0, 'tau_m': 15.0, 'tau_m_K_Ca': 1.0,
+                 'w_BAP': 27.995561755479308}
+
+other_params = {'V_th': -40.,
+                'tau_w': 500.,
+                'a': 0.,
+                'b': 40.,
+                }
+
+w_BAP = 27.995561755479308
+d_BAP = 0.1195980511869619
+
+cm_params = default_param
+cm_params.update(other_params)
+
+# Following lines can be used to change specific parameters from defaul ones
+cm_params.update({
+    #'w_BAP': 0., 
+    #'g_C_d':5.,
+})
+
+print('')
+print('Updated Ca-AdEd params: ', cm_params)
+print('')
+
+SimTime = 6000
+brain_state = "awake"
+stimulusStart = 0.
+stimulusStop = SimTime
+countWindow = stimulusStop - stimulusStart
+
+I_s = 100
+I_d = 50
+
+res_cm = []
+
+# rec_list = ['V_m_s', 'V_m_d', 'w',
+#             'm_Ca', 'h_Ca', 'I_Ca', 'c_Ca', 'e_Ca',
+#             'm_K', 'I_K']
+rec_list = ['V_m_s', 'V_m_d', 'w',
+            'm_Ca', 'h_Ca', 'c_Ca',
+            'm_K']
+conn_spec = {"rule": "all_to_all"}
+
+# Creation of actual simulation starts from here
+cm = ngpu.Create("ca_adex_alt_nestml", 1)
+ngpu.SetStatus(cm, cm_params)
+params = ngpu.GetStatus(cm)
+print(dict(sorted(params[0].items())))
+ngpu.ActivateSpikeCount(cm)
+ngpu.ActivateRecSpikeTimes(cm, 100)
+
+receptor_types = {
+    "EXC_SPIKES_SOMA": 0,
+    "INH_SPIKES_SOMA": 1,
+    "EXC_SPIKES_DISTAL": 2,
+    "INH_SPIKES_DISTAL": 3,
+    "SPIKES_AMPA_D": 4,
+    "SPIKES_BETA_D": 5,
+    "SPIKES_AMPA_NMDA_D": 6
+}
+# connect soma to distal with delay to reproduce back propagation currents
+ngpu.Connect(cm, cm, conn_spec, syn_dict={'weight': w_BAP, 'delay': d_BAP,
+                                          'receptor': receptor_types["SPIKES_AMPA_D"]})
+
+sg0 = ngpu.Create('spike_generator', 1, status_dict={'spike_times': [500]})
+sg1 = ngpu.Create('spike_generator', 1, status_dict={'spike_times': [1500]})
+sg2 = ngpu.Create('spike_generator', 1, status_dict={'spike_times': [2500]})
+sg3 = ngpu.Create('spike_generator', 1, status_dict={'spike_times': [3500]})
+sg4 = ngpu.Create('spike_generator', 1, status_dict={'spike_times': [4500]})
+ngpu.Connect(sg0, cm, conn_spec, syn_dict={
+    'weight': 300, 'delay': 1., 'receptor': receptor_types['EXC_SPIKES_SOMA']})
+ngpu.Connect(sg1, cm, conn_spec, syn_dict={
+    'weight': 1500, 'delay': 1.,
+    'receptor': receptor_types['EXC_SPIKES_DISTAL']})
+ngpu.Connect(sg2, cm, conn_spec, syn_dict={
+    'weight': 500, 'delay': 1., 'receptor': receptor_types['EXC_SPIKES_SOMA']})
+ngpu.Connect(sg2, cm, conn_spec, syn_dict={
+    'weight': 300, 'delay': 1., 'receptor': receptor_types['SPIKES_AMPA_D']})
+ngpu.Connect(sg3, cm, conn_spec, syn_dict={
+    'weight': 300, 'delay': 1., 'receptor': receptor_types['EXC_SPIKES_SOMA']})
+ngpu.Connect(sg3, cm, conn_spec, syn_dict={
+    'weight': 50, 'delay': 1.,
+    'receptor': receptor_types['SPIKES_AMPA_NMDA_D']})
+ngpu.Connect(sg4, cm, conn_spec, syn_dict={
+    'weight': 100, 'delay': 1.,
+    'receptor': receptor_types['SPIKES_AMPA_NMDA_D']})
+
+record = ngpu.CreateRecord("", rec_list, [cm[0]]*len(rec_list), [0]*len(rec_list))
+
+ngpu.Simulate(SimTime)
+data = ngpu.GetRecordData(record)
+nest_gpu_data = np.array(data)
+
+spike_count = ngpu.GetStatus(cm, "spike_count")
+print("Spike count: ", spike_count)
+spike_times = ngpu.GetRecSpikeTimes(cm)
+print("Spike times: ", spike_times)
+
+# NEST Data
+nest_data = np.loadtxt("ca_adex_nest_data.txt", delimiter="\t")
+
+xlimStart = 0
+xlimEnd = SimTime
+
+plt.figure('Ca-AdEx', figsize=(10, 8))
+plt.subplot(411)
+###############################################################################
+plt.plot(nest_data[:, 0], nest_data[:, 1], c='b', label='v_m soma [NEST]')
+plt.plot(nest_gpu_data[:, 0], nest_gpu_data[:, 1], c='orange', ls="--", label='v_m soma [NEST-GPU]')
+
+plt.legend(loc='upper right')
+plt.xlim(xlimStart, xlimEnd)
+plt.ylabel('Vm soma [mV]')
+# plt.title('Voltage (soma)')
+
+plt.subplot(412)
+plt.plot(nest_data[:, 0], nest_data[:, 2], c='g', label='v_m dist cm [NEST]')
+plt.plot(nest_gpu_data[:, 0], nest_gpu_data[:, 2], c='orange', ls="--", label='v_m dist [NEST-GPU]')
+plt.legend(loc='upper right')
+plt.xlim(xlimStart, xlimEnd)
+plt.ylabel('Vm distal [mV]')
+# plt.title('Voltage (distal)')
+
+plt.subplot(413)
+plt.plot(nest_data[:, 0], nest_data[:, 4], c='b', lw=2., label='m [NEST]')
+plt.plot(nest_gpu_data[:, 0], nest_gpu_data[:, 4], c='orange', ls='--', lw=2., label='m [NEST-GPU]')
+
+plt.plot(nest_data[:, 0], nest_data[:, 5], c='r', lw=2., label='h [NEST]')
+plt.plot(nest_gpu_data[:, 0], nest_gpu_data[:, 5], c='orange', ls='--', lw=2., label='h [NEST-GPU]')
+
+# plt.plot(nest_data[:, 0], nest_data[:, 4] * nest_data[:, 5], c='k', lw=2., label='g [NEST]')
+# plt.plot(nest_gpu_data[:, 0], nest_gpu_data[:, 4] * nest_gpu_data[:, 5], c='orange', ls='--', lw=2., label='g [NEST-GPU]')
+
+plt.legend(loc='upper right')
+plt.xlim(xlimStart, xlimEnd)
+plt.xlim(xlimStart, xlimEnd)
+plt.ylabel('Ca')
+# plt.title('Distal Ca activation')
+
+plt.subplot(414)
+plt.plot(nest_data[:, 0], nest_data[:, 3], c='b', lw=2., label='W [NEST]')
+plt.plot(nest_gpu_data[:, 0], nest_gpu_data[:, 3], c='orange', ls='--', lw=2., label='W [NEST-GPU]')
+plt.legend(loc='upper right')
+plt.xlim(xlimStart, xlimEnd)
+plt.xlim(xlimStart, xlimEnd)
+plt.ylabel('W')
+plt.xlabel('Time [ms]')
+# plt.title('Adaptation')
+plt.show()
+# plt.savefig('Ca-AdEx-1.png')
+
+# plt.subplot(414)
+# plt.plot(res_cm[0]['times'], res_cm[0]['e_Ca'], c='b', lw=2., label='e_Ca [conv]')
+# plt.plot(res_cm[1]['times'], res_cm[1]['e_Ca'], c='orange', ls='--', lw=2., label='e_Ca')
+#
+# plt.ylabel('e_Ca')
+# plt.title('Ca reversal potential')
+# plt.xlim(xlimStart, xlimEnd)
+# plt.xlabel('Time [ms]')
+
+plt.figure('NESTML SC DIST', figsize=(10, 5))
+###############################################################################
+# plt.subplot(511)
+# plt.plot(res_cm[0]['times'], res_cm[0]['I_Ca'], c='b', label='I_Ca [conv]')
+# plt.plot(res_cm[1]['times'], res_cm[1]['I_Ca'], c='orange', ls='--', label='I_Ca')
+#
+# plt.legend()
+# plt.xlim(xlimStart, xlimEnd)
+
+# plt.subplot(512)
+# plt.plot(nest_data[:, 0], nest_data[0], c='b', lw=2., label='m_Ca [conv]')
+# plt.plot(nest_gpu_data[1]['times'], nest_gpu_data[1]['m_Ca'], c='orange', ls='--', lw=2., label='m_Ca')
+#
+# plt.plot(nest_data[0]['times'], nest_data[0]['h_Ca'], c='r', lw=2., label='h_Ca [conv]')
+# plt.plot(nest_gpu_data[1]['times'], nest_gpu_data[1]['h_Ca'], c='orange', ls='--', lw=2., label='h_Ca')
+#
+# plt.plot(nest_data[0]['times'], nest_data[0]['m_Ca'] * nest_data[0]['h_Ca'], c='k', lw=2., label='g_Ca [conv]')
+# plt.plot(nest_gpu_data[1]['times'], nest_gpu_data[1]['m_Ca'] * nest_gpu_data[1]['h_Ca'], c='orange', ls='--', lw=2., label='g_Ca')
+#
+# plt.legend()
+# plt.xlim(xlimStart, xlimEnd)
+
+plt.subplot(211)
+plt.plot(nest_data[:, 0], nest_data[:, 6], c='b', label='[Ca] [NEST]')
+plt.plot(nest_gpu_data[:, 0], nest_gpu_data[:, 6], c='orange', ls='--', label='[Ca] [NEST-GPU]')
+
+plt.legend()
+plt.xlim(xlimStart, xlimEnd)
+
+# plt.subplot(514)
+# plt.plot(res_cm[0]['times'], res_cm[0]['I_K'], c='b', label='I_K')
+# plt.plot(res_cm[1]['times'], res_cm[1]['I_K'], c='orange', ls='--', label='I_K')
+#
+# plt.legend()
+# plt.xlim(xlimStart, xlimEnd)
+
+plt.subplot(212)
+plt.plot(nest_data[:, 0], nest_data[:, 7], c='b', lw=2., label='m_K [NEST]')
+plt.plot(nest_gpu_data[:, 0], nest_gpu_data[:, 7], c='orange', ls='--', lw=2., label='m_K [NEST-GPU]')
+
+plt.legend()
+plt.xlim(xlimStart, xlimEnd)
+plt.xlabel('Time [ms]')
+
+# plt.show()
+plt.savefig('Ca-AdEx-2.png')

examples/nest_gpu/Ca-AdEx_Example5.py

@@ -0,0 +1,140 @@
+#
+#  Ca-AdEx: example of usage of two-compartment neuron model
+#
+#  First version: 22/10/2025
+#  Author: Elena Pastorelli, INFN, Rome (IT)
+#
+#  Description: Example of usage of the Ca-AdEx neuron stimulated with spike generators
+#
+#  Copyright © 2025   Elena Pastorelli          <elena.pastorelli@roma1.infn.it>
+#  Copyright © 2025   Pier Stanislao Paolucci   <pier.paolucci@roma1.infn.it>
+#
+#  SPDX-License-Identifier: GPL-3.0-only
+#
+#  This program is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import nest
+import nestgpu as ngpu
+
+default_param = {'C_m_d': 23.67372778891213, 'C_m_s': 246.7882968598874, 'Ca_0': 0.0001, 'Ca_th': 0.00043,
+                 'V_reset': -61.73952230767877, 'd_BAP': 0.1195980511869619, 'delta_T': 2.0, 'e_K': -90.0,
+                 'e_L_d': -55.000000000000014, 'e_L_s': -69.24596493128396, 'e_Na_Adex': -50.0, 'exp_K_Ca': 4.8,
+                 'g_C_d': 19.777320239615996, 'g_L_d': 3.377855016658499, 'g_L_s': 5.0, 'g_w': 1.1156385639067352,
+                 'gbar_Ca': 21.045506331690845, 'gbar_K_Ca': 13.199867205029523, 'h_half': -21.0, 'h_slope': -0.5,
+                 'm_half': -9.0, 'm_slope': 0.5, 'phi': 3.92830985228413e-08, 't_ref': 0.0,
+                 'tau_decay_Ca': 103.57233790866408, 'tau_h': 80.0, 'tau_m': 15.0, 'tau_m_K_Ca': 1.0,
+                 'w_BAP': 27.995561755479308}
+
+other_params = {'V_th': -40.,
+                'tau_w': 500.,
+                'a': 0.,
+                'b': 40.,
+                }
+
+w_BAP = 27.995561755479308
+d_BAP = 0.1195980511869619
+
+cm_params = default_param
+cm_params.update(other_params)
+
+SimTime = 5000
+brain_state = "awake"
+
+conn_spec = {"rule": "all_to_all"}
+
+Jnoise_ex = 30.0
+p_rate = 250.0
+
+
+# NEST simulation
+def simulate_nest():
+    nest.ResetKernel()
+
+    nest.Install('ca_adex_module')
+
+    cm = nest.Create("ca_adex_alt", params=cm_params)
+    receptor_types = cm.get("receptor_types")
+    print(" ")
+    print("cm neuron properties: ", dict(sorted(cm.get().items())))
+
+    # Poisson generators
+    noise1 = nest.Create("poisson_generator", params={"rate": p_rate})
+    nest.CopyModel("static_synapse", "noise_ex")
+    syn_spec_noise1 = {"synapse_model": "noise_ex", "weight": Jnoise_ex, "delay": 1.0}
+    syn_spec_noise1.update({"weight": Jnoise_ex, 'receptor_type': receptor_types['EXC_SPIKES_SOMA']})
+
+    # Spike recorder
+    sr = nest.Create("spike_recorder")
+
+    # connect soma to distal with delay to reproduce back propagation currents
+    nest.Connect(cm, cm, syn_spec={'synapse_model': 'static_synapse',
+                                   'weight': w_BAP, 'delay': d_BAP,
+                                   'receptor_type': receptor_types["SPIKES_AMPA_D"]})
+
+    nest.Connect(noise1, cm, conn_spec=conn_spec, syn_spec=syn_spec_noise1)
+    nest.Connect(cm, sr)
+
+    nest.Simulate(SimTime)
+
+    spike_count = sr.n_events
+    return spike_count
+
+
+def simulate_nest_gpu():
+    # Creation of actual simulation starts from here
+    cm = ngpu.Create("ca_adex_alt_nestml", 1)
+    ngpu.SetStatus(cm, cm_params)
+    ngpu.ActivateSpikeCount(cm)
+
+    receptor_types = {
+        "EXC_SPIKES_SOMA": 0,
+        "INH_SPIKES_SOMA": 1,
+        "EXC_SPIKES_DISTAL": 2,
+        "INH_SPIKES_DISTAL": 3,
+        "SPIKES_AMPA_D": 4,
+        "SPIKES_BETA_D": 5,
+        "SPIKES_AMPA_NMDA_D": 6
+    }
+
+    # Poisson generators
+    noise1 = ngpu.Create("poisson_generator")
+    ngpu.SetStatus(noise1, "rate", p_rate)
+
+    syn_spec_noise1 = {"weight": Jnoise_ex, "delay": 1.0, 'receptor': receptor_types['EXC_SPIKES_SOMA']}
+
+    # connect soma to distal with delay to reproduce back propagation currents
+    ngpu.Connect(cm, cm, conn_spec, syn_dict={'weight': w_BAP, 'delay': d_BAP,
+                                              'receptor': receptor_types["SPIKES_AMPA_D"]})
+    ngpu.Connect(noise1, cm, conn_dict=conn_spec, syn_dict=syn_spec_noise1)
+
+    ngpu.Simulate(SimTime)
+
+    spike_count = ngpu.GetStatus(cm, "spike_count")
+    return spike_count[0][0]
+
+
+# Simulate the neuron in NEST and NEST GPU using poisson generators
+spike_count_nest = simulate_nest()
+spike_count_nest_gpu = simulate_nest_gpu()
+
+rate_nest = spike_count_nest / SimTime * 1000.0
+rate_nest_gpu = spike_count_nest_gpu / SimTime * 1000.0
+
+print("######## NEST ########")
+print("Spike count: ", spike_count_nest)
+print("Firing rate: ", rate_nest)
+
+print("######## NEST GPU ########")
+print("Spike count: ", spike_count_nest_gpu)
+print("Firing rate: ", rate_nest_gpu)

nestml/build_ca_adex_gpu.py

@@ -0,0 +1,8 @@
+from pynestml.frontend.pynestml_frontend import generate_nest_gpu_target
+
+codegen_opts = {"solver": "numeric"}
+generate_nest_gpu_target(input_path="ca_adex_alt.nestml",
+                         target_path="target_gpu",
+                         logging_level="INFO",
+                         suffix="_nestml",
+                         codegen_opts=codegen_opts)