New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Revise implementation of ht_neuron #491

Merged
merged 48 commits into from Dec 16, 2016

Conversation

@heplesser
Contributor

heplesser commented Sep 22, 2016

Implement the full NMDA model from the Hill-Tononi paper with instantaneous blocking and fast and slow unblocking. This PR will change model behavior.

As reviewers, I suggest @ingablundell and @suku248.

@Silmathoron

Ok, the changes make sense for me, I'll try to test the model at some point this week.
There are a few points that need discussing, though.

Show outdated Hide outdated models/ht_connection.h
, tau_P_( 50.0 )
, delta_P_( 0.2 )
, tau_P_( 500.0 )
, delta_P_( 0.5 )

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Since this could potentially have an influence on existing user scripts, could you justify the change in default parameters?

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Since this could potentially have an influence on existing user scripts, could you justify the change in default parameters?

This comment has been minimized.

@heplesser

heplesser Oct 18, 2016

Contributor

The actual values of these two parameters are not given in the Hill-Tononi (2005) paper. The new values are the default values in the original Synthesis simulator code by Sean Hill. I sent a mail to NEST User on 22 September warning of changes to the Hill-Tononi model and asking for comments from anyone who may be affected, but have not received a single answer, so I think we are safe here.

@heplesser

heplesser Oct 18, 2016

Contributor

The actual values of these two parameters are not given in the Hill-Tononi (2005) paper. The new values are the default values in the original Synthesis simulator code by Sean Hill. I sent a mail to NEST User on 22 September warning of changes to the Hill-Tononi model and asking for comments from anyone who may be affected, but have not received a single answer, so I think we are safe here.

This comment has been minimized.

@Silmathoron

Silmathoron Oct 18, 2016

Contributor

Ok, sounds fine!

@Silmathoron

Silmathoron Oct 18, 2016

Contributor

Ok, sounds fine!

Show outdated Hide outdated models/ht_neuron.cpp
@@ -49,8 +49,7 @@ RecordablesMap< ht_neuron >::create()
Name( "Theta" ), &ht_neuron::get_y_elem_< ht_neuron::State_::THETA > );
insert_(
Name( "g_AMPA" ), &ht_neuron::get_y_elem_< ht_neuron::State_::G_AMPA > );
insert_(
Name( "g_NMDA" ), &ht_neuron::get_y_elem_< ht_neuron::State_::G_NMDA > );
insert_( Name( "g_NMDA" ), &ht_neuron::get_g_NMDA_ );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Why is the recordable g_NMDA different from the state variable? (cf. get_g_NMDA)

EDIT: ok, I understand, the state variable G_NMDA is the g(t) from the docstring and the recordable is the real g_NMDA = g(t)*m(V, t), which cannot be stored directly as a state variable since it is not what's updated.
But then the naming looks misleading to me, I guess we should change G_NMDA to something else and update the docstring... though I must say I'm at a loss regarding the name. Maybe G_T would be the most sensible choice...

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

Why is the recordable g_NMDA different from the state variable? (cf. get_g_NMDA)

EDIT: ok, I understand, the state variable G_NMDA is the g(t) from the docstring and the recordable is the real g_NMDA = g(t)*m(V, t), which cannot be stored directly as a state variable since it is not what's updated.
But then the naming looks misleading to me, I guess we should change G_NMDA to something else and update the docstring... though I must say I'm at a loss regarding the name. Maybe G_T would be the most sensible choice...

This comment has been minimized.

@heplesser

heplesser Oct 18, 2016

Contributor

Naming is a bit of a mess here. In the paper, it is \tilde{g}_NMDA = m(V) g_NMDA(t), where g_NMDA(t) is the double exponential function.

I think that the end user will be most interested in the actual conductance of NMDA channels when Mg-blocking is taken into account (the m(V)). Therefore, I think the recordable should be g_NMDA and the getter get_g_NMDA_(), but I think I will change the name of the state variable, maybe to G_NMDA_TIMECOURSE.

@heplesser

heplesser Oct 18, 2016

Contributor

Naming is a bit of a mess here. In the paper, it is \tilde{g}_NMDA = m(V) g_NMDA(t), where g_NMDA(t) is the double exponential function.

I think that the end user will be most interested in the actual conductance of NMDA channels when Mg-blocking is taken into account (the m(V)). Therefore, I think the recordable should be g_NMDA and the getter get_g_NMDA_(), but I think I will change the name of the state variable, maybe to G_NMDA_TIMECOURSE.

This comment has been minimized.

@Silmathoron

Silmathoron Oct 18, 2016

Contributor

Yes, I agree that g_NMDA should be the variable recorded, and G_NMDA_TIMECOURSE looks good!

@Silmathoron

Silmathoron Oct 18, 2016

Contributor

Yes, I agree that g_NMDA should be the variable recorded, and G_NMDA_TIMECOURSE looks good!

Show outdated Hide outdated models/ht_neuron.cpp
* NMDA conductance
*
* We need to take care to handle instantaneous blocking correctly.
* If the unblock-variables Mg_{fast,slow} > Mg_ss, the steady-state

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

This looks unclear to me, you mean "if the unblock variables are greater than their respective steady-state values", right? I think the combination of equation and text sentence with commas is rather disturbing.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

This looks unclear to me, you mean "if the unblock variables are greater than their respective steady-state values", right? I think the combination of equation and text sentence with commas is rather disturbing.

This comment has been minimized.

@heplesser

heplesser Oct 18, 2016

Contributor

You are right, fixed.

@heplesser

heplesser Oct 18, 2016

Contributor

You are right, fixed.

Show outdated Hide outdated models/ht_neuron.cpp
, NMDA_Vact( -58.0 ) // mV
, NMDA_Sact( 2.5 ) // mV
, NMDA_Tau_1( 4.0 ) // ms
, NMDA_Tau_2( 40.0 ) // ms

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

If you think it's worth it (I'm not sure), having all taus in lowercase would be more coherent, though it would break existing codes... and Tau_m instead of tau_m is really too bad since it's supposed to be common to many models...

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

If you think it's worth it (I'm not sure), having all taus in lowercase would be more coherent, though it would break existing codes... and Tau_m instead of tau_m is really too bad since it's supposed to be common to many models...

@@ -280,7 +300,10 @@ class ht_neuron : public Archiving_Node
DG_GABA_A,

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

I would be in favour of changing G_NMDA to G_T but there should at least be a comment like
G_NMDA, // state variable corresponding to g(t), with g_NMDA = m(V, t) * g(t)
or something equivalent.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

I would be in favour of changing G_NMDA to G_T but there should at least be a comment like
G_NMDA, // state variable corresponding to g(t), with g_NMDA = m(V, t) * g(t)
or something equivalent.

Show outdated Hide outdated models/ht_neuron.h
@@ -409,6 +432,26 @@ class ht_neuron : public Archiving_Node
{
return S_.I_h_;
}
double_t

This comment has been minimized.

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

why are there double_t?

@Silmathoron

Silmathoron Oct 17, 2016

Contributor

why are there double_t?

This comment has been minimized.

@heplesser

heplesser Oct 18, 2016

Contributor

I suppose because the branch is a bit old .... Fixed.

@heplesser

heplesser Oct 18, 2016

Contributor

I suppose because the branch is a bit old .... Fixed.

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Oct 18, 2016

Contributor

@Silmathoron I now noticed that the model needs a bit more work. I will change all names to use names:: names, will add checks to parameter setting (currently, there are no checks at all) and change all _Tau_ and _Theta_ in names to lowercase equivalents.

But I probably won't get that done today.

Contributor

heplesser commented Oct 18, 2016

@Silmathoron I now noticed that the model needs a bit more work. I will change all names to use names:: names, will add checks to parameter setting (currently, there are no checks at all) and change all _Tau_ and _Theta_ in names to lowercase equivalents.

But I probably won't get that done today.

@Silmathoron

This comment has been minimized.

Show comment
Hide comment
@Silmathoron

Silmathoron Oct 18, 2016

Contributor

Ok, thanks for your hard work!
May the code be with you!

Contributor

Silmathoron commented Oct 18, 2016

Ok, thanks for your hard work!
May the code be with you!

Show outdated Hide outdated models/ht_neuron.cpp
* We need to take care to handle instantaneous blocking correctly.
* If the unblock-variables Mg_{fast,slow} > Mg_ss, the steady-state
* value for the present membrane potential, we cannot change those values
* in State_[], since the ODE Solver may call this function multiple times

This comment has been minimized.

@ingablundell

ingablundell Oct 21, 2016

I would prefer it if "those values" were specified

@ingablundell

ingablundell Oct 21, 2016

I would prefer it if "those values" were specified

This comment has been minimized.

@heplesser

heplesser Oct 24, 2016

Contributor

I will fix this with the next update to the PR.

@heplesser

heplesser Oct 24, 2016

Contributor

I will fix this with the next update to the PR.

@ingablundell

This comment has been minimized.

Show comment
Hide comment
@ingablundell

ingablundell Oct 21, 2016

The documentation of the NMDA current is very clear and helpful! Although I would add an explicit definition of g (what is g_NMDA in the paper). The calculation of the NMDA part of I_syn and propagation of Mg_fast and Mg_slow look correct to me! I approve this pull request.

The documentation of the NMDA current is very clear and helpful! Although I would add an explicit definition of g (what is g_NMDA in the paper). The calculation of the NMDA part of I_syn and propagation of Mg_fast and Mg_slow look correct to me! I approve this pull request.

Show outdated Hide outdated models/ht_neuron.h
m(V) = 1 / ( 1 + exp( - ( V - NMDA_Vact ) / NMDA_Sact ) )
where g(t) is a beta function (difference of exponentials).

This comment has been minimized.

@ingablundell

ingablundell Oct 21, 2016

I would prefer it if g was stated explicitly

@ingablundell

ingablundell Oct 21, 2016

I would prefer it if g was stated explicitly

This comment has been minimized.

@heplesser

heplesser Oct 24, 2016

Contributor

I will fix this with the next update to the PR.

@heplesser

heplesser Oct 24, 2016

Contributor

I will fix this with the next update to the PR.

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Oct 26, 2016

Contributor

@Silmathoron I have now completed my revisions, so you can continue reviewing.

Contributor

heplesser commented Oct 26, 2016

@Silmathoron I have now completed my revisions, so you can continue reviewing.

@Silmathoron

Silmathoron suggested changes Oct 27, 2016 edited

Nice restructuring!
Here is a first round of comments, I'll try to check the equations and test the model next week:

  • in ht_connection.h, line 44, it should be tau_P (or tau_rec, if possible, cf. comments below).
  • the THIS MODEL NEURON HAS NOT BEEN TESTED EXTENSIVELY! comment leaves me perplex... what do you mean exactly? Wouldn't a test file solve the problem?
def< double >( d, "tau_P", tau_P_ );
def< double >( d, "delta_P", delta_P_ );
def< double >( d, "P", p_ );
def< double >( d, names::tau_P, tau_P_ );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

Since we are trying to limit the number of names and given that this is clearly a recovery time, I think it would be a good idea to set it to names::tau_rec and specify in the docstring "the recovery time tau_P in the paper is called tau_rec here", or something equivalent

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

Since we are trying to limit the number of names and given that this is clearly a recovery time, I think it would be a good idea to set it to names::tau_rec and specify in the docstring "the recovery time tau_P in the paper is called tau_rec here", or something equivalent

This comment has been minimized.

@heplesser

heplesser Nov 23, 2016

Contributor

I see your point, but since everything here is about P, and I would like to stay reasonably close to the paper to ease comparisons, I would like to stick with tau_P.

@heplesser

heplesser Nov 23, 2016

Contributor

I see your point, but since everything here is about P, and I would like to stay reasonably close to the paper to ease comparisons, I would like to stick with tau_P.

This comment has been minimized.

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

Fair enough ;)

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

Fair enough ;)

Show outdated Hide outdated models/ht_neuron.cpp
def< double >( d, names::theta_eq, theta_eq );
def< double >( d, names::tau_theta, tau_theta );
def< double >( d, names::tau_spike, tau_spike );
def< double >( d, names::spike_duration, t_spike );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

t_spike is not detailed anywhere, could you add a line of explanation in the docstring?

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

t_spike is not detailed anywhere, could you add a line of explanation in the docstring?

Show outdated Hide outdated models/ht_neuron.h
//! Timer (counter) for potassium current.
int r_potassium_;
//! Timer (counter) for spike-activated repolarizing potassium current.
int r_spike_;
bool g_spike_; //!< active / not active

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

The name g_spike_ seems weird for a boolean telling whether the K current is in action or not. Also this variable is not documented though it's a recordable...

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

The name g_spike_ seems weird for a boolean telling whether the K current is in action or not. Also this variable is not documented though it's a recordable...

Show outdated Hide outdated models/ht_neuron.h
- Intrinsic currents I_h (pacemaker), I_T (low-threshold calcium),
I_Na(p) (persistent sodium), and I_KNa (depolarization-activated
potassium).
- Several apparent typographical errors in the descriptions of

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

I guess you mean errors in the paper [1], but maybe it should be more explicit...

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

I guess you mean errors in the paper [1], but maybe it should be more explicit...

Show outdated Hide outdated models/ht_neuron.h
t_peak = tau_2 tau_1 / ( tau_2 - tau_1 ) ln( tau_2 / tau_1 )
m(V, t) = a(V) m_fast*(V, t) + ( 1 - a(V) ) m_slow*(V, t)
a(V) = 0.51 - 0.0028 V

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

There is an apparent problem with the homogeneity, here... I suppose 0.0028 is in mV^{-1}, but could you specify it explicitely? (maybe introduce a named constant).
Shouldn't the user want to be able to set these two values? If not (which would be surprising to me), would you mind explaining their origin?

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

There is an apparent problem with the homogeneity, here... I suppose 0.0028 is in mV^{-1}, but could you specify it explicitely? (maybe introduce a named constant).
Shouldn't the user want to be able to set these two values? If not (which would be surprising to me), would you mind explaining their origin?

This comment has been minimized.

@heplesser

heplesser Nov 23, 2016

Contributor

I sure should document this better and provide details in a notebook. Given the number of parameters in this model that are fixed as numerical constants, making them all settable would make the code more complex (just inventing good names for all constants would be a major effort). Therefore, I would like to leave them as they are (seems to be quite usual for these more complex models). The particular values here are from Vargas-Caballero and Robinson, J Neurophysiol 89: 2778–2783, 2003. I will add the reference.

@heplesser

heplesser Nov 23, 2016

Contributor

I sure should document this better and provide details in a notebook. Given the number of parameters in this model that are fixed as numerical constants, making them all settable would make the code more complex (just inventing good names for all constants would be a major effort). Therefore, I would like to leave them as they are (seems to be quite usual for these more complex models). The particular values here are from Vargas-Caballero and Robinson, J Neurophysiol 89: 2778–2783, 2003. I will add the reference.

Show outdated Hide outdated nestkernel/nest_names.cpp
@@ -155,6 +175,8 @@ const Name growth_rate( "growth_rate" );
const Name gsl_error_tol( "gsl_error_tol" );
const Name h( "h" );
const Name h_E_rev( "h_E_rev" );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

same comment as above

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

same comment as above

Show outdated Hide outdated nestkernel/nest_names.cpp
@@ -179,6 +207,9 @@ const Name Interpol_Order( "Interpol_Order" );
const Name interval( "interval" );
const Name is_refractory( "is_refractory" );
const Name KNa_E_rev( "KNa_E_rev" );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

see above

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

see above

Show outdated Hide outdated nestkernel/nest_names.cpp
@@ -202,8 +233,18 @@ const Name N_channels( "N_channels" );
const Name n_events( "n_events" );
const Name n_proc( "n_proc" );
const Name n_receptors( "n_receptors" );
const Name NaP_E_rev( "NaP_E_rev" );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

see above

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

see above

Show outdated Hide outdated nestkernel/nest_names.cpp
const Name needs_prelim_update( "needs_prelim_update" );
const Name neuron( "neuron" );
const Name NMDA_E_rev( "NMDA_E_rev" );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

see above

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

see above

Show outdated Hide outdated nestkernel/nest_names.cpp
@@ -288,6 +332,8 @@ const Name synapse_model( "synapse_model" );
const Name synapse_modelid( "synapse_modelid" );
const Name synaptic_elements( "synaptic_elements" );
const Name T_E_rev( "T_E_rev" );

This comment has been minimized.

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

see above

@Silmathoron

Silmathoron Oct 27, 2016

Contributor

see above

- Revised parameter names in ht_neuron from AMPA_E_rev to E_rev_AMPA etc
- Fixed some other parameter names in ht_neuron: tau_1,2 -> tau_rise,decay; spike_duration -> t_ref
- Recordable names now g_GABA_A, g_GABA_B
- Corrected refratory handling: no spike emission during t_ref
@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Nov 25, 2016

Contributor

@Silmathoron The Notebook on the model is still work in progess. I have pushed it nonetheless for your information during the review.

Contributor

heplesser commented Nov 25, 2016

@Silmathoron The Notebook on the model is still work in progess. I have pushed it nonetheless for your information during the review.

@Silmathoron

Ok, almost everything looks good to me except some lingering t_spikes that need to disappear and a few remarks on the notebook.
Thanks a lot for the impressive restructuring, the code looks much clearer now, and the notebook really makes it a lot easier to understand too!
If you can just add axis names to the existent (and future) figures, it will be perfect ;)

double h = e.get_stamp().get_ms() - t_lastspike;
Node* target = get_target( t );
// t_lastspike_ = 0 initially
// propagation t_lastspike -> t_spike, t_lastspike_ = 0 initially, p_ = 1

This comment has been minimized.

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

change t_spike to t_ref

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

change t_spike to t_ref

This comment has been minimized.

@heplesser

heplesser Nov 25, 2016

Contributor

No, here t_spike is actually the time of the spike passing through the synapse.

@heplesser

heplesser Nov 25, 2016

Contributor

No, here t_spike is actually the time of the spike passing through the synapse.

This comment has been minimized.

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

Ok, right, I went a bit fast on that one ^^"

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

Ok, right, I went a bit fast on that one ^^"

Show outdated Hide outdated models/ht_neuron.cpp
}
if ( t_ref < 0 )
{
throw BadParameter( "t_spike >= 0 required." );

This comment has been minimized.

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

change t_spike to t_ref

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

change t_spike to t_ref

Show outdated Hide outdated models/ht_neuron.cpp
V_.PotassiumRefractoryCounts_ = Time( Time::ms( P_.t_spike ) ).get_steps();
V_.PotassiumRefractoryCounts_ = Time( Time::ms( P_.t_ref ) ).get_steps();
// since t_spike_ >= 0, this can only fail in error

This comment has been minimized.

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

change t_spike to t_ref

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

change t_spike to t_ref

Show outdated Hide outdated models/ht_neuron.cpp
S_.y_[ State_::Mg_slow ] = std::min( Mg_ss, S_.y_[ State_::Mg_slow ] );
S_.y_[ State_::Mg_fast ] = std::min( Mg_ss, S_.y_[ State_::Mg_fast ] );
// Deactivate potassium current after t_spike

This comment has been minimized.

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

change t_spike to t_ref

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

change t_spike to t_ref

Show outdated Hide outdated models/ht_neuron.h
THIS MODEL NEURON HAS NOT BEEN TESTED EXTENSIVELY!
Parameters:
V_m - membrane potential
spike_duration - duration of re-polarizing potassium current
Tau_m - membrane time constant applying to all currents but
t_spike - duration of re-polarizing potassium current

This comment has been minimized.

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

remove t_spike

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

remove t_spike

Show outdated Hide outdated doc/model_details/HillTononi.ipynb
"\n",
"Note the following:\n",
"- $N_T=2$ from Synthesis code\n",
"- In the equation for $\\tau_{m,T}$, the second exponential term must be added to the first (in the denominator) to make dimensional sense\n",

This comment has been minimized.

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

I guess this means that the constants 0.22 and 0.13 have the dimension of a time? could you precise this?
EDIT: ok, this appears just below... but still ^^

@Silmathoron

Silmathoron Nov 25, 2016

Contributor

I guess this means that the constants 0.22 and 0.13 have the dimension of a time? could you precise this?
EDIT: ok, this appears just below... but still ^^

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Nov 25, 2016

Contributor

All should be fixed now. I am continuing work on the notebook, next step is a systematic test of the intrinsic currents (voltage-clamping is coming to NEST ...)

Contributor

heplesser commented Nov 25, 2016

All should be fixed now. I am continuing work on the notebook, next step is a systematic test of the intrinsic currents (voltage-clamping is coming to NEST ...)

heplesser added some commits Nov 28, 2016

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Nov 29, 2016

Contributor

@Silmathoron @ingablundell I have just pushed a major revision of the model including a notebook exploring all aspects of the model (docs/model_details). Enjoy! Note that failing tests may be due to Travis trouble this afternoon.

Contributor

heplesser commented Nov 29, 2016

@Silmathoron @ingablundell I have just pushed a major revision of the model including a notebook exploring all aspects of the model (docs/model_details). Enjoy! Note that failing tests may be due to Travis trouble this afternoon.

@Silmathoron

Great work on that, the update function is much clearer now!
There are only a few points left that need some clarifications; also I realized that I forgot to ask, but I find it very surprising that, in such a detailed model, conductances cannot be taken directly from biological data... could you at least provide typical values that would allow for an approximate conversion?

Show outdated Hide outdated doc/model_details/HillTononiModels.ipynb
"\\end{equation}\n",
"This equation is also given in a comment in Synthesis, but is missing from the implementation.\n",
"\n",
"**Note: NEST implements the equation according to Compte et al (2003) with $N_{NaP}=3$, while Synthesis uses $N_{NaP}=1**\n",

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Missing $ closure after $N_{NaP}=1

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Missing $ closure after $N_{NaP}=1

This comment has been minimized.

@heplesser

heplesser Dec 1, 2016

Contributor

Done.

@heplesser

heplesser Dec 1, 2016

Contributor

Done.

Show outdated Hide outdated doc/model_details/HillTononiModels.ipynb
"\n",
"#### Synaptic \"minis\"\n",
"\n",
"Synaptic \"minis\" due to spontaneous release of neurotransmitter quanta [HT05, p 1679] are not included in the NEST implementation of the Hill-Tononi model. This due to the fact that the total mini input rate for a cell was just 2 Hz and they cause PSP changes by $0.5 \\pm 0.25$mV only and thus should have minimal effect.\n",

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

For people working on cultures (like me) this is a very important information, so maybe it would be good to make it appear clearly at the beginning of the notebook.
(also minis are usually proportional to the neuron's in-degree, so this constant value is a bit surprising to me... could you specify the conditions?)

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

For people working on cultures (like me) this is a very important information, so maybe it would be good to make it appear clearly at the beginning of the notebook.
(also minis are usually proportional to the neuron's in-degree, so this constant value is a bit surprising to me... could you specify the conditions?)

This comment has been minimized.

@heplesser

heplesser Dec 1, 2016

Contributor

From the paper: "In addition, we modeled “minis” ... as low-amplitude PSPs (mean 0.5+-0.25 mV) .... The mean frequency of these Poisson distributed synaptic minis was set to 2 Hz (total for an individual cell)."

@heplesser

heplesser Dec 1, 2016

Contributor

From the paper: "In addition, we modeled “minis” ... as low-amplitude PSPs (mean 0.5+-0.25 mV) .... The mean frequency of these Poisson distributed synaptic minis was set to 2 Hz (total for an individual cell)."

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Ok, thanks

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Ok, thanks

Show outdated Hide outdated models/ht_neuron.cpp
* function multiple times and in arbitrary temporal order. We thus need to
* use local variables for the values at the current time, and check the
* state variables once the ODE solver has completed the time step.
* If the unblock-variables m_NMDA_{fast,slow} are greater than the

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Naming is not in agreement with lines 87--88 (m_{fast,slow}_NMDA)

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Naming is not in agreement with lines 87--88 (m_{fast,slow}_NMDA)

This comment has been minimized.

@heplesser

heplesser Dec 1, 2016

Contributor

Done.

@heplesser

heplesser Dec 1, 2016

Contributor

Done.

Show outdated Hide outdated models/ht_neuron.cpp
const double D_influx =
D_influx_peak / ( 1.0 + std::exp( -( V - D_thresh ) / D_slope ) );
return P_.tau_D_KNa * D_influx + D_eq_0;

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

D_eq_0 was not present before and the value returned by the function is apparently not homogeneous... could you explain the role of this new variable?

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

D_eq_0 was not present before and the value returned by the function is apparently not homogeneous... could you explain the role of this new variable?

This comment has been minimized.

@heplesser

heplesser Dec 1, 2016

Contributor

Renamed D_eq as in the paper/equations.

@heplesser

heplesser Dec 1, 2016

Contributor

Renamed D_eq as in the paper/equations.

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

What about the homogeneity?
I noticed that, indeed, the equation for D does not seem to satisfy this requirement, unless D_influx has a dimension of $T^{-1}$, which does not seem to be the case...

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

What about the homogeneity?
I noticed that, indeed, the equation for D does not seem to satisfy this requirement, unless D_influx has a dimension of $T^{-1}$, which does not seem to be the case...

This comment has been minimized.

@heplesser

heplesser Dec 2, 2016

Contributor

I hadn't thought about this, but you are right. The dimension of D does not matter, as long as it is the same as of d_1/2, since the two are divided in the end, so dimensions drop out. I suggest to leave D and d_1/2 dimensionless, since they represent some "factor" which has no immediate physical counterpart (but models some sort of ion concentration). Then D_eq also is dimensionless, while D_influx and D_influx,peak must have units of inverse time, explicitly 1/ms. This makes sense, since these values essentially represent ion concentration change. I have updated the notebook.

@heplesser

heplesser Dec 2, 2016

Contributor

I hadn't thought about this, but you are right. The dimension of D does not matter, as long as it is the same as of d_1/2, since the two are divided in the end, so dimensions drop out. I suggest to leave D and d_1/2 dimensionless, since they represent some "factor" which has no immediate physical counterpart (but models some sort of ion concentration). Then D_eq also is dimensionless, while D_influx and D_influx,peak must have units of inverse time, explicitly 1/ms. This makes sense, since these values essentially represent ion concentration change. I have updated the notebook.

This comment has been minimized.

@Silmathoron

Silmathoron Dec 2, 2016

Contributor

Ok, good!

@Silmathoron

Silmathoron Dec 2, 2016

Contributor

Ok, good!

Show outdated Hide outdated models/ht_neuron.cpp
def< double >( d, names::theta, y_[ THETA ] ); // Threshold
}
void
nest::ht_neuron::State_::set( const DictionaryDatum& d, const Parameters_& )
nest::ht_neuron::State_::set( const DictionaryDatum& d,
const ht_neuron& node,

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

I don't understand this modification of the function

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

I don't understand this modification of the function

This comment has been minimized.

@heplesser

heplesser Dec 1, 2016

Contributor

Several of the intrinsic channels have slow dynamics, O(1s) for I_h, e.g. If you change parameters, especially the membrane potential using SetStatus/SetDefaults, it will often be useful to set the state variables for these channels to their new steady-state values. Otherwise, one would have to simulate several seconds of biological time to get to equilibrium. Especially when setting model parameters at the beginning of a simulation, the new equilibration capability should be useful. When changing parameters during the simulation, on the other hand, you probably want to see how the system responds to the change (eg my voltage stepping tests and the wake-sleep transitions in HT05), so then one will not want to equilibrate.

@heplesser

heplesser Dec 1, 2016

Contributor

Several of the intrinsic channels have slow dynamics, O(1s) for I_h, e.g. If you change parameters, especially the membrane potential using SetStatus/SetDefaults, it will often be useful to set the state variables for these channels to their new steady-state values. Otherwise, one would have to simulate several seconds of biological time to get to equilibrium. Especially when setting model parameters at the beginning of a simulation, the new equilibration capability should be useful. When changing parameters during the simulation, on the other hand, you probably want to see how the system responds to the change (eg my voltage stepping tests and the wake-sleep transitions in HT05), so then one will not want to equilibrate.

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

I'm sorry, my previous statement was unclear! I don't understand why you added const ht_neuron& node as parameter, could you not simply call the functions?

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

I'm sorry, my previous statement was unclear! I don't understand why you added const ht_neuron& node as parameter, could you not simply call the functions?

This comment has been minimized.

@heplesser

heplesser Dec 2, 2016

Contributor

set() is a member function of State_, not of ht_neuron. Therefore, I cannot directly call member functions of ht_neuron. But since I now call methods in ht_neuron to get the equilibrium values, I don't need the Parameters_ explicitly and will remove that argument.

@heplesser

heplesser Dec 2, 2016

Contributor

set() is a member function of State_, not of ht_neuron. Therefore, I cannot directly call member functions of ht_neuron. But since I now call methods in ht_neuron to get the equilibrium values, I don't need the Parameters_ explicitly and will remove that argument.

This comment has been minimized.

@Silmathoron

Silmathoron Dec 2, 2016

Contributor

Of course, I read too fast, sorry...

@Silmathoron

Silmathoron Dec 2, 2016

Contributor

Of course, I read too fast, sorry...

@@ -577,7 +624,7 @@ nest::ht_neuron::Buffers_::Buffers_( const Buffers_&, ht_neuron& n )
nest::ht_neuron::ht_neuron()
: Archiving_Node()
, P_()
, S_( P_ )
, S_( *this, P_ )

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

See above

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

See above

This comment has been minimized.

@heplesser

heplesser Dec 1, 2016

Contributor

... and same answer, essentially: the state should be initialized to the correct steady-state values.

@heplesser

heplesser Dec 1, 2016

Contributor

... and same answer, essentially: the state should be initialized to the correct steady-state values.

This comment has been minimized.

@heplesser

heplesser Dec 2, 2016

Contributor

Here I need both parameters, since I compute some equilibrium values explicitly in the constructor.

@heplesser

heplesser Dec 2, 2016

Contributor

Here I need both parameters, since I compute some equilibrium values explicitly in the constructor.

{E_rev,g_peak}_{h,T,NaP,KNa} - reversal potential and peak conductance for
intrinsic currents
receptor_types - dictionary mapping synapse names to ports on
neuron model
recordables - list of recordable quantities.
recordables - list of recordable quantities
equilibrate - if given and true, time-dependent activation

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Maybe a precision that values are otherwise initialized at zero would be good

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Maybe a precision that values are otherwise initialized at zero would be good

This comment has been minimized.

@heplesser

heplesser Dec 1, 2016

Contributor

I added and "otherwise" part.

@heplesser

heplesser Dec 1, 2016

Contributor

I added and "otherwise" part.

Show outdated Hide outdated models/ht_neuron.h
@@ -285,12 +288,15 @@ class ht_neuron : public Archiving_Node
double g_peak_KNa;
double E_rev_KNa; // mV
double tau_D_KNa; // ms; in P_ for technical reasons, not meant to be public

This comment has been minimized.

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Why not?

@Silmathoron

Silmathoron Dec 1, 2016

Contributor

Why not?

This comment has been minimized.

@heplesser

heplesser Dec 1, 2016

Contributor

I turned it into a normal, accessible parameter now.

@heplesser

heplesser Dec 1, 2016

Contributor

I turned it into a normal, accessible parameter now.

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Dec 1, 2016

Contributor

I have just pushed a new version addressing @Silmathoron's comments and also

  • integrating insights from #547, especially moving threshold check into the ODE-integration loop and handling refractoriness in continuous time ("inside" approach)
  • shortened the documentation in the code, since the details are now in the notebook; the details on NMDA-unblocking seemed out of balance, as no other, more basic equations were shown.
Contributor

heplesser commented Dec 1, 2016

I have just pushed a new version addressing @Silmathoron's comments and also

  • integrating insights from #547, especially moving threshold check into the ODE-integration loop and handling refractoriness in continuous time ("inside" approach)
  • shortened the documentation in the code, since the details are now in the notebook; the details on NMDA-unblocking seemed out of balance, as no other, more basic equations were shown.
@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Dec 2, 2016

Contributor

@ingablundell Does your "thumbs up" still stand after the changes I made in response to @Silmathoron's comments?

Contributor

heplesser commented Dec 2, 2016

@ingablundell Does your "thumbs up" still stand after the changes I made in response to @Silmathoron's comments?

Show outdated Hide outdated models/ht_neuron.h
@@ -337,7 +300,7 @@ class ht_neuron : public Archiving_Node
/** Timer (counter) for spike-activated repolarizing potassium current.
* Neuron is absolutely refractory during this period.
*/
int ref_steps_;
long ref_steps_;

This comment has been minimized.

@Silmathoron

Silmathoron Dec 2, 2016

Contributor

I never thought about this before, but it could indeed be a good idea to generalize this change...

@Silmathoron

Silmathoron Dec 2, 2016

Contributor

I never thought about this before, but it could indeed be a good idea to generalize this change...

@Silmathoron

Silmathoron approved these changes Dec 2, 2016 edited

All clear for me, thanks for the good work! 👍

@heplesser heplesser changed the title from Implement full NMDA dynamics in ht_neuron to Revise implementation of ht_neuron Dec 5, 2016

@Silmathoron

Two minor details left ;)

Show outdated Hide outdated doc/model_details/HillTononiModels.ipynb
"- $m_{DK}^{\\infty}$ is a steep sigmoid which is almost 0 or 1 except for a narrow window around $d_{1/2}$.\n",
"- To the left of this window, $I_{DK}\\approx 0$.\n",
"- To the right of this window, $I_{DK}\\sim -(V-E_{DK})$.\n",
"The equations implemented in NEST are thus\n",

This comment has been minimized.

@Silmathoron

Silmathoron Dec 6, 2016

Contributor

Either some other equation are missing or this is a wrong copy-paste

@Silmathoron

Silmathoron Dec 6, 2016

Contributor

Either some other equation are missing or this is a wrong copy-paste

This comment has been minimized.

@heplesser

heplesser Dec 7, 2016

Contributor

Yes, that was a copy-paste leftover, fixed.

@heplesser

heplesser Dec 7, 2016

Contributor

Yes, that was a copy-paste leftover, fixed.

" m_{DK}^{\\infty} &= \\frac{1}{1 + \\left(\\frac{d_{1/2}}{D}\\right)^{3.5}}\\\\\n",
" \\frac{dD}{dt} &= D_{\\text{influx}}(V) - \\frac{D-D_{\\text{eq}}}{\\tau_D} = \\frac{D_{\\infty}(V)-D}{\\tau_D} \\\\\n",
" D_{\\infty}(V) &= \\tau_D D_{\\text{influx}}(V) + {D_{\\text{eq}}}\\\\\n",
" D_{\\text{influx}} &= \\frac{D_{\\text{influx,peak}}}{1+ \\exp\\left(-\\frac{V-D_{\\theta}}{\\sigma_D}\\right)} \n",

This comment has been minimized.

@Silmathoron

Silmathoron Dec 6, 2016

Contributor

Could you add in the "Note the following" list what you told me about the fact that, despite their ambiguous names, D_{influx} and D_{influx,peak}

represent ion concentration change

If you could elaborate a little on this here, I think it would be a good thing.

@Silmathoron

Silmathoron Dec 6, 2016

Contributor

Could you add in the "Note the following" list what you told me about the fact that, despite their ambiguous names, D_{influx} and D_{influx,peak}

represent ion concentration change

If you could elaborate a little on this here, I think it would be a good thing.

This comment has been minimized.

@heplesser

heplesser Dec 7, 2016

Contributor

Done.

@heplesser

heplesser Dec 7, 2016

Contributor

Done.

@Silmathoron

👍

@ingablundell

ingablundell approved these changes Dec 15, 2016 edited

I really like the notebook! I approve this pull request!

Show outdated Hide outdated doc/model_details/HillTononi.ipynb
"### Integration \n",
"\n",
"The original Synthesis implementation of the model uses Runge-Kutta integration with fixed 0.25 ms step size, and integrates channels dynamics first, followed by integration of membrane potential and threshold.\n",
"\n",

This comment has been minimized.

@ingablundell

ingablundell Dec 15, 2016

Maybe you could state which RK solver you are using?

@ingablundell

ingablundell Dec 15, 2016

Maybe you could state which RK solver you are using?

This comment has been minimized.

@heplesser

heplesser Dec 15, 2016

Contributor

Done.

@heplesser

heplesser Dec 15, 2016

Contributor

Done.

Show outdated Hide outdated doc/model_details/HillTononi.ipynb
"\n",
"- The equation does not contain membrane capacitance. As a side-effect, all conductances are dimensionless.\n",
"- Na and K leak conductances $g_{\\text{NaL}}$ and $g_{\\text{KL}}$ are constant, although $g_{\\text{KL}}$ may be adjusted on slow time scales to mimic neuromodulatory effects.\n",
"- Reversal potentials are assumed constant.\n",

This comment has been minimized.

@ingablundell

ingablundell Dec 15, 2016

I would prefer: Reversal potentials $E_{NA}$ and $E_{K}$ are assumed constant.

@ingablundell

ingablundell Dec 15, 2016

I would prefer: Reversal potentials $E_{NA}$ and $E_{K}$ are assumed constant.

This comment has been minimized.

@heplesser

heplesser Dec 15, 2016

Contributor

Done.

@heplesser

heplesser Dec 15, 2016

Contributor

Done.

Show outdated Hide outdated doc/model_details/HillTononi.ipynb
"Synthesis: `ItChannel`\n",
"\n",
"##### Equations given in paper \n",
"\n",

This comment has been minimized.

@ingablundell

ingablundell Dec 15, 2016

Why are you calling it "equations" and not "definitions"?

@ingablundell

ingablundell Dec 15, 2016

Why are you calling it "equations" and not "definitions"?

This comment has been minimized.

@heplesser

heplesser Dec 15, 2016

Contributor

Because physicist are a sloppy bunch and call most things with a "=" in the middle an equation. Fixing this here may introduce inconsistencies elsewhere in the notebook, so I prefer to leave it as it is.

@heplesser

heplesser Dec 15, 2016

Contributor

Because physicist are a sloppy bunch and call most things with a "=" in the middle an equation. Fixing this here may introduce inconsistencies elsewhere in the notebook, so I prefer to leave it as it is.

Show outdated Hide outdated doc/model_details/HillTononi.ipynb
"Intrinsic currents are based on the Hodgkin-Huxley description, i.e.,\n",
"\n",
"\\begin{align}\n",
"I_X &= g_{\\text{peak}, X} m_X(V, t)^N_X h_X(V, t)(V-E_X) \\\\\n",

This comment has been minimized.

@ingablundell

ingablundell Dec 15, 2016

I would prefer it if the meaning of the variables were stated explicitly. But maybe this is not necessary.

@ingablundell

ingablundell Dec 15, 2016

I would prefer it if the meaning of the variables were stated explicitly. But maybe this is not necessary.

This comment has been minimized.

@heplesser

heplesser Dec 15, 2016

Contributor

Done.

@heplesser

heplesser Dec 15, 2016

Contributor

Done.

Show outdated Hide outdated doc/model_details/HillTononi.ipynb
"\n",
"#### Synaptic \"minis\"\n",
"\n",
"Synaptic \"minis\" due to spontaneous release of neurotransmitter quanta [HT05, p 1679] are not included in the NEST implementation of the Hill-Tononi model. This due to the fact that the total mini input rate for a cell was just 2 Hz and they cause PSP changes by $0.5 \\pm 0.25$mV only and thus should have minimal effect.\n",

This comment has been minimized.

@ingablundell

ingablundell Dec 15, 2016

This is due...

This comment has been minimized.

@heplesser

heplesser Dec 15, 2016

Contributor

Change that to a much simpler "..., because the total mini input rate".

@heplesser

heplesser Dec 15, 2016

Contributor

Change that to a much simpler "..., because the total mini input rate".

Show outdated Hide outdated doc/model_details/HillTononi.ipynb
" t_{\\text{peak}} &= \\frac{\\tau_2 \\tau_1}{\\tau_2 - \\tau_1} \\ln\\frac{ \\tau_2}{\\tau_1}\n",
"\\end{align} \n",
"\n",
"where $t_{\\text{peak}}$ is the time of the conductance maximum and $\\tau_1$ and $\\tau_2$ are synaptic rise- and decay-time, respectively; $\\Theta(t)$ is the Heaviside step function. The equation is integrated using exact integration in Synthesis; in NEST, it is included in the ODE-system integrated using RK.\n",

This comment has been minimized.

@ingablundell

ingablundell Dec 15, 2016

again: which RK?

@ingablundell

ingablundell Dec 15, 2016

again: which RK?

This comment has been minimized.

@heplesser

heplesser Dec 15, 2016

Contributor

Done.

@heplesser

heplesser Dec 15, 2016

Contributor

Done.

@heplesser

This comment has been minimized.

Show comment
Hide comment
@heplesser

heplesser Dec 15, 2016

Contributor

@ingablundell Thank you! I have revised the notebook according to your suggestions, with one exception.

@apeyser Could you merge, as this is my PR?

Contributor

heplesser commented Dec 15, 2016

@ingablundell Thank you! I have revised the notebook according to your suggestions, with one exception.

@apeyser Could you merge, as this is my PR?

@abigailm abigailm merged commit b13c085 into nest:master Dec 16, 2016

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details

@heplesser heplesser deleted the heplesser:HT_NMDA branch Dec 20, 2016

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment