Add models with exact and simplified NMDA dynamics according to Brunel and Wang (2001)

modelsets/full

@@ -60,6 +60,7 @@ iaf_psc_exp_multisynapse
 iaf_psc_exp_ps
 iaf_psc_exp_ps_lossless
 iaf_tum_2000
+iaf_wang_2002
 izhikevich
 jonke_synapse
 lin_rate

nestkernel/nest_names.cpp

@@ -106,6 +106,7 @@ const Name continuous( "continuous" );
 const Name count_covariance( "count_covariance" );
 const Name count_histogram( "count_histogram" );
 const Name covariance( "covariance" );
+const Name conc_Mg2( "conc_Mg2" );
 
 const Name Delta_T( "Delta_T" );
 const Name Delta_V( "Delta_V" );
@@ -178,10 +179,13 @@ const Name filenames( "filenames" );
 const Name frequency( "frequency" );
 const Name frozen( "frozen" );
 
+const Name GABA( "GABA" );
 const Name GABA_A( "GABA_A" );
 const Name GABA_B( "GABA_B" );
 const Name g( "g" );
 const Name g_AMPA( "g_AMPA" );
+const Name g_AMPA_ext( "g_AMPA_ext" );
+const Name g_GABA( "g_GABA" );
 const Name g_ahp( "g_ahp" );
 const Name g_C( "g_C" );
 const Name g_ex( "g_ex" );
@@ -308,10 +312,10 @@ const Name mult_coupling( "mult_coupling" );
 const Name music_channel( "music_channel" );
 
 const Name N( "N" );
+const Name NMDA( "NMDA" );
 const Name N_channels( "N_channels" );
 const Name N_NaP( "N_NaP" );
 const Name N_T( "N_T" );
-const Name NMDA( "NMDA" );
 const Name n( "n" );
 const Name n_events( "n_events" );
 const Name n_messages( "n_messages" );
@@ -414,6 +418,10 @@ const Name rule( "rule" );
 
 const Name S( "S" );
 const Name S_act_NMDA( "S_act_NMDA" );
+const Name s_NMDA( "s_NMDA" );
+const Name s_AMPA( "s_AMPA" );
+const Name s_AMPA_ext( "s_AMPA_ext" );
+const Name s_GABA( "s_GABA" );
 const Name SIC_scale( "SIC_scale" );
 const Name SIC_th( "SIC_th" );
 const Name sdev( "sdev" );
@@ -478,11 +486,17 @@ const Name targets( "targets" );
 const Name tau( "tau" );
 const Name tau_1( "tau_1" );
 const Name tau_2( "tau_2" );
-const Name tau_ahp( "tau_ahp" );
+const Name tau_AMPA( "tau_AMPA" );
 const Name tau_Ca( "tau_Ca" );
-const Name tau_c( "tau_c" );
 const Name tau_D_KNa( "tau_D_KNa" );
 const Name tau_Delta( "tau_Delta" );
+const Name tau_GABA( "tau_GABA" );
+const Name tau_Mg_fast_NMDA( "tau_Mg_fast_NMDA" );
+const Name tau_Mg_slow_NMDA( "tau_Mg_slow_NMDA" );
+const Name tau_P( "tau_P" );
+const Name tau_V_th( "tau_V_th" );
+const Name tau_ahp( "tau_ahp" );
+const Name tau_c( "tau_c" );
 const Name tau_decay( "tau_decay" );
 const Name tau_decay_AMPA( "tau_decay_AMPA" );
 const Name tau_decay_ex( "tau_decay_ex" );
@@ -493,15 +507,12 @@ const Name tau_decay_NMDA( "tau_decay_NMDA" );
 const Name tau_epsp( "tau_epsp" );
 const Name tau_fac( "tau_fac" );
 const Name tau_IP3( "tau_IP3" );
-const Name tau_Mg_fast_NMDA( "tau_Mg_fast_NMDA" );
-const Name tau_Mg_slow_NMDA( "tau_Mg_slow_NMDA" );
 const Name tau_m( "tau_m" );
 const Name tau_max( "tau_max" );
 const Name tau_minus( "tau_minus" );
 const Name tau_minus_stdp( "tau_minus_stdp" );
 const Name tau_minus_triplet( "tau_minus_triplet" );
 const Name tau_n( "tau_n" );
-const Name tau_P( "tau_P" );
 const Name tau_plus( "tau_plus" );
 const Name tau_plus_triplet( "tau_plus_triplet" );
 const Name tau_psc( "tau_psc" );
@@ -527,7 +538,6 @@ const Name tau_theta( "tau_theta" );
 const Name tau_u_bar_bar( "tau_u_bar_bar" );
 const Name tau_u_bar_minus( "tau_u_bar_minus" );
 const Name tau_u_bar_plus( "tau_u_bar_plus" );
-const Name tau_V_th( "tau_V_th" );
 const Name tau_v( "tau_v" );
 const Name tau_vacant( "tau_vacant" );
 const Name tau_w( "tau_w" );

nestkernel/nest_names.h

@@ -122,6 +122,7 @@ extern const Name clear;
 extern const Name comp_idx;
 extern const Name comparator;
 extern const Name compartments;
+extern const Name conc_Mg2;
 extern const Name configbit_0;
 extern const Name configbit_1;
 extern const Name connection_count;
@@ -204,12 +205,15 @@ extern const Name filenames;
 extern const Name frequency;
 extern const Name frozen;
 
+extern const Name GABA;
 extern const Name GABA_A;
 extern const Name GABA_B;
 extern const Name g;
 extern const Name g_AMPA;
+extern const Name g_AMPA_ext;
 extern const Name g_ahp;
 extern const Name g_C;
+extern const Name g_GABA;
 extern const Name g_ex;
 extern const Name g_GABA_A;
 extern const Name g_GABA_B;
@@ -334,6 +338,8 @@ extern const Name mult_coupling;
 extern const Name music_channel;
 
 extern const Name N;
+extern const Name NMDA;
+extern const Name NMDA_sum;
 extern const Name N_channels;
 extern const Name N_NaP;
 extern const Name N_T;
@@ -440,6 +446,10 @@ extern const Name rule;
 
 extern const Name S;
 extern const Name S_act_NMDA;
+extern const Name s_GABA;
+extern const Name s_AMPA;
+extern const Name s_AMPA_ext;
+extern const Name s_NMDA;
 extern const Name SIC_scale;
 extern const Name SIC_th;
 extern const Name sdev;
@@ -504,6 +514,8 @@ extern const Name targets;
 extern const Name tau;
 extern const Name tau_1;
 extern const Name tau_2;
+extern const Name tau_AMPA;
+extern const Name tau_GABA;
 extern const Name tau_ahp;
 extern const Name tau_Ca;
 extern const Name tau_c;

pynest/examples/wang_decision_making.py

@@ -0,0 +1,181 @@
+import nest
+import matplotlib.pyplot as plt
+import numpy as np
+
+np.random.seed(123)
+rng = np.random.default_rng()
+
+# Parameters from paper
+epop_params = {"g_AMPA_ext": 2.1,
+               "g_AMPA": 0.05,
+               "g_NMDA": 0.165,
+               "g_GABA": 1.3,
+               "tau_GABA": 5.0,
+               "tau_AMPA": 2.0,
+               "tau_decay_NMDA": 100.0,
+               "alpha": 0.5,
+               "conc_Mg2": 1.0,
+               "g_L": 25.,               # leak conductance
+               "E_L": -70.0,            # leak reversal potential
+               "E_ex": 0.0,             # excitatory reversal potential
+               "E_in": -70.0,           # inhibitory reversal potential
+               "V_reset": -55.0,        # reset potential
+               "V_th": -50.0,           # threshold
+               "C_m": 500.0,            # membrane capacitance
+               "t_ref": 2.0             # refreactory period
+               }
+
+
+ipop_params = {"g_AMPA_ext": 1.62,
+               "g_AMPA": 0.04,
+               "g_NMDA": 0.13,
+               "g_GABA": 1.0,
+               "tau_GABA": 5.0,
+               "tau_AMPA": 2.0,
+               "tau_decay_NMDA": 100.0,
+               "alpha": 0.5,
+               "conc_Mg2": 1.0,
+               "g_L": 20.,               # leak conductance
+               "E_L": -70.0,            # leak reversal potential
+               "E_ex": 0.0,             # excitatory reversal potential
+               "E_in": -70.0,           # inhibitory reversal potential
+               "V_reset": -55.0,        # reset potential
+               "V_th": -50.0,           # threshold
+               "C_m": 200.0,            # membrane capacitance
+               "t_ref": 1.0             # refreactory period
+               }
+
+simtime = 4000.
+signal_start = 1000.
+signal_duration = 2000.
+signal_update_interval = 50.
+f = 0.15 # proportion of neurons receiving signal inputs
+w_plus = 1.7
+w_minus = 1 - f * (w_plus - 1) / (1 - f)
+
+NE = 1600
+NI = 400
+
+selective_pop1 = nest.Create("iaf_wang_2002", int(0.15 * NE), params=epop_params)
+selective_pop2 = nest.Create("iaf_wang_2002", int(0.15 * NE), params=epop_params)
+nonselective_pop = nest.Create("iaf_wang_2002", int(0.7 * NE), params=epop_params)
+inhibitory_pop = nest.Create("iaf_wang_2002", NI, params=ipop_params)
+
+mu_0 = 40.
+rho_a = mu_0 / 100
+rho_b = rho_a
+c = 80.
+sigma = 4.
+mu_a = mu_0 + rho_a * c
+mu_b = mu_0 - rho_b * c
+
+num_updates = int(signal_duration / signal_update_interval)
+update_times = np.arange(0, signal_duration, signal_update_interval)
+update_times[0] = 0.1
+rates_a = np.random.normal(mu_a, sigma, size=num_updates)
+rates_b = np.random.normal(mu_b, sigma, size=num_updates)
+
+poisson_a = nest.Create("inhomogeneous_poisson_generator",
+                        params={"origin": signal_start-0.1,
+                                "start": 0.,
+                                "stop": signal_duration,
+                                "rate_times": update_times,
+                                "rate_values": rates_a})
+
+poisson_b = nest.Create("inhomogeneous_poisson_generator",
+                        params={"origin": signal_start-0.1,
+                                "start": 0.,
+                                "stop": signal_duration,
+                                "rate_times": update_times,
+                                "rate_values": rates_b})
+
+poisson_0 = nest.Create("poisson_generator", params={"rate": 2400.})
+
+syn_spec_pot = {"synapse_model": "static_synapse", "weight":w_plus, "delay":0.1, 'receptor_type': 0}
+syn_spec_default = {"synapse_model": "static_synapse", "weight":1.0, "delay":0.1, 'receptor_type': 0}
+syn_spec_dep = {"synapse_model": "static_synapse", "weight":w_minus, "delay":0.1, 'receptor_type': 0}
+syn_spec_inhibitory = {"synapse_model": "static_synapse", "weight":-1.0, "delay":0.1, 'receptor_type': 0}
+syn_spec_ext = {"synapse_model": "static_synapse", "weight":1., "delay":0.1, 'receptor_type': 1}
+
+sr_nonselective = nest.Create("spike_recorder")
+sr_selective1 = nest.Create("spike_recorder")
+sr_selective2 = nest.Create("spike_recorder")
+sr_inhibitory = nest.Create("spike_recorder")
+
+mm_selective1 = nest.Create("multimeter", {"record_from": ["V_m", "s_NMDA", "s_AMPA", "s_AMPA_ext", "s_GABA"]})
+mm_inhibitory = nest.Create("multimeter", {"record_from": ["V_m", "s_NMDA", "s_AMPA", "s_GABA"]})
+
+
+nest.Connect(poisson_0,
+             nonselective_pop + selective_pop1 + selective_pop2 + inhibitory_pop,
+             conn_spec="all_to_all",
+             syn_spec=syn_spec_ext)
+
+nest.Connect(nonselective_pop,
+             selective_pop1 + selective_pop2,
+             conn_spec="all_to_all",
+             syn_spec=syn_spec_dep)
+
+nest.Connect(nonselective_pop,
+             nonselective_pop + inhibitory_pop,
+             conn_spec="all_to_all",
+             syn_spec=syn_spec_default)
+
+nest.Connect(selective_pop1,
+             selective_pop2,
+             conn_spec="all_to_all",
+             syn_spec=syn_spec_dep)
+
+nest.Connect(selective_pop2,
+             selective_pop1,
+             conn_spec="all_to_all",
+             syn_spec=syn_spec_dep)
+
+nest.Connect(selective_pop1 + selective_pop2,
+             nonselective_pop + inhibitory_pop,
+             conn_spec="all_to_all",
+             syn_spec=syn_spec_default)
+
+nest.Connect(inhibitory_pop,
+             selective_pop1 + selective_pop2 + nonselective_pop + inhibitory_pop,
+             conn_spec="all_to_all",
+             syn_spec=syn_spec_inhibitory)
+
+nest.Connect(poisson_a,
+             selective_pop1,
+             conn_spec="all_to_all",
+             syn_spec=syn_spec_ext)
+
+nest.Connect(poisson_b,
+             selective_pop2,
+             conn_spec="all_to_all",
+             syn_spec=syn_spec_ext)
+
+nest.Connect(nonselective_pop, sr_nonselective)
+nest.Connect(selective_pop1, sr_selective1)
+nest.Connect(selective_pop2, sr_selective2)
+nest.Connect(inhibitory_pop, sr_inhibitory)
+
+
+nest.Connect(mm_selective1, selective_pop1)
+nest.Connect(mm_inhibitory, inhibitory_pop)
+
+nest.Simulate(4000.)
+
+spikes_nonselective = sr_nonselective.get("events", "times")
+spikes_selective1 = sr_selective1.get("events", "times")
+spikes_selective2 = sr_selective2.get("events", "times")
+spikes_inhibitory = sr_inhibitory.get("events", "times")
+
+
+
+
+
+senders = mm_selective1.get("events", "senders")
+inds = senders == 1
+vm = mm_selective1.get("events", "V_m")[inds]
+s_AMPA = mm_selective1.get("events", "s_AMPA")[inds]
+s_AMPA_ext = mm_selective1.get("events", "s_AMPA_ext")[inds]
+s_GABA = mm_selective1.get("events", "s_GABA")[inds]
+s_NMDA = mm_selective1.get("events", "s_NMDA")[inds]
+