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

Correct input handling for aeif_cond_beta_multisynapse #510 and Update aeif_cond_beta_multisynapse to new aeif numerics approach #514 #516

Merged
merged 25 commits into from Oct 26, 2016
Commits
Jump to file or symbol
Failed to load files and symbols.
+109 −81
Diff settings

Always

Just for now

Next

Added reversal potential array E_rev, with one value of the reversal

potential for each receptor port.
Fixed the way how the number of receptor ports, taus_rise,
taus_decay and E_rev can be updated.
Fixed the pynest example in aeif_cond_beta_multisynapse.h .
  • Loading branch information...
golosio committed Oct 7, 2016
commit 5006d162bf7e7e815c4a28bfbb404f1dbf940b30
@@ -76,23 +76,29 @@ aeif_cond_beta_multisynapse::Parameters_::Parameters_()
, t_ref_( 0.0 ) // ms
, g_L( 30.0 ) // nS
, C_m( 281.0 ) // pF
, E_ex( 0.0 ) // mV
, E_in( -85.0 ) // mV
, E_L( -70.6 ) // mV
, Delta_T( 2.0 ) // mV
, tau_w( 144.0 ) // ms
, a( 4.0 ) // nS
, b( 80.5 ) // pA
, V_th( -50.4 ) // mV
, num_of_receptors( 1 )

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

I wonder if we should have this as a member. Since it will always be the same as E_ref.size(), I think we should should just use

inline size_t n_receptors() { return E_rev.size(); }
@heplesser

heplesser Oct 13, 2016

Contributor

I wonder if we should have this as a member. Since it will always be the same as E_ref.size(), I think we should should just use

inline size_t n_receptors() { return E_rev.size(); }

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

replaced with inline function as you suggest.

@golosio

golosio Oct 14, 2016

Contributor

replaced with inline function as you suggest.

, I_e( 0.0 ) // pA
, gsl_error_tol( 1e-6 )
, num_of_receptors_( 1 )
, has_connections_( false )
{
static const double E_rev_def = 0.0; // E_rev default value
static const double taus_rise_def = 2.0; // taus_rise default value
static const double taus_decay_def = 20.0; // taus_decay default value
receptor_types.clear();
E_rev.clear();
taus_rise.clear();
taus_decay.clear();

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

No need to clear here, we are in a constructor!

@heplesser

heplesser Oct 13, 2016

Contributor

No need to clear here, we are in a constructor!

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

removed.

@golosio

golosio Oct 14, 2016

Contributor

removed.

taus_decay.push_back( 20.0 );
taus_rise.push_back( 2.0 );
receptor_types.push_back( 1 );

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

I am actually not sure we need or use the receptor_types member anywhere. If it is not used, we should drop it.

@heplesser

heplesser Oct 13, 2016

Contributor

I am actually not sure we need or use the receptor_types member anywhere. If it is not used, we should drop it.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

Right, it was not used anywhere, I removed it.

@golosio

golosio Oct 14, 2016

Contributor

Right, it was not used anywhere, I removed it.

E_rev.push_back( E_rev_def );
taus_rise.push_back( taus_rise_def );
taus_decay.push_back( taus_decay_def );
}
aeif_cond_beta_multisynapse::State_::State_( const Parameters_& p )
@@ -131,10 +137,11 @@ aeif_cond_beta_multisynapse::Parameters_::get( DictionaryDatum& d ) const
def< double >( d, names::g_L, g_L );
def< double >( d, names::E_L, E_L );
def< double >( d, names::V_reset, V_reset_ );
def< double >( d, names::E_ex, E_ex );
def< double >( d, names::E_in, E_in );
def< size_t >( d, names::num_of_receptors, num_of_receptors );
ArrayDatum E_rev_ad( E_rev );
ArrayDatum taus_rise_ad( taus_rise );
ArrayDatum taus_decay_ad( taus_decay );
def< ArrayDatum >( d, names::E_rev, E_rev_ad );
def< ArrayDatum >( d, names::taus_rise, taus_rise_ad );
def< ArrayDatum >( d, names::taus_decay, taus_decay_ad );
def< double >( d, names::a, a );
@@ -144,68 +151,70 @@ aeif_cond_beta_multisynapse::Parameters_::get( DictionaryDatum& d ) const
def< double >( d, names::I_e, I_e );
def< double >( d, names::V_peak, V_peak_ );
def< double >( d, names::gsl_error_tol, gsl_error_tol );
def< int >( d, "n_synapses", num_of_receptors_ );
def< bool >( d, names::has_connections, has_connections_ );
}
void
aeif_cond_beta_multisynapse::Parameters_::set( const DictionaryDatum& d )
{
static const double E_rev_def = 0.0;
static const double taus_rise_def = 2.0;
static const double taus_decay_def = 20.0;
updateValue< double >( d, names::V_th, V_th );
updateValue< double >( d, names::V_peak, V_peak_ );
updateValue< double >( d, names::t_ref, t_ref_ );
updateValue< double >( d, names::E_L, E_L );
updateValue< double >( d, names::V_reset, V_reset_ );
updateValue< double >( d, names::E_ex, E_ex );
updateValue< double >( d, names::E_in, E_in );
updateValue< double >( d, names::C_m, C_m );
updateValue< double >( d, names::g_L, g_L );
std::vector< double > tau_tmp;
if ( updateValue< std::vector< double > >( d, names::taus_decay, tau_tmp ) )
size_t tmp_num;
if ( updateValue< long >( d, names::num_of_receptors, tmp_num ) )
{
if ( tau_tmp.size() < taus_decay.size() && has_connections_ == true )
if ( tmp_num < num_of_receptors && has_connections_ == true )
{
throw BadProperty(
"The neuron has connections, therefore the number of ports cannot be "
"reduced." );
}
for ( size_t i = 0; i < tau_tmp.size(); ++i )
else if ( tmp_num <= 0 )
{
if ( tau_tmp[ i ] <= 0 )
{
throw BadProperty(
"All synaptic time constants must be strictly positive" );
}
throw BadProperty(
"The neuron must have at least one port." );
}
taus_decay = tau_tmp;
num_of_receptors_ = taus_decay.size();
tau_tmp.clear();
if ( taus_rise.size() < taus_decay.size() )
{
for ( size_t i = 0; i < taus_decay.size(); ++i )
{
taus_rise.push_back( taus_decay[ i ] / 10. );
// if taus_rise is not defined explicitly or if it has less elements
// than taus_decay, it will be set to taus_decay/10
}
E_rev.resize( tmp_num, E_rev_def );
taus_rise.resize( tmp_num, taus_rise_def );
taus_decay.resize( tmp_num, taus_decay_def );
receptor_types.resize( tmp_num );
for ( size_t i = num_of_receptors; i < tmp_num; i++ )
{
receptor_types[i] = i + 1 ;
}
num_of_receptors = tmp_num;
}
if ( updateValue< std::vector< double > >( d, names::taus_rise, tau_tmp ) )
std::vector< double > E_rev_tmp;
if ( updateValue< std::vector< double > >( d, names::E_rev, E_rev_tmp ) )
{
if ( taus_decay.size() == 0 )
if ( E_rev_tmp.size() != num_of_receptors )
{
throw BadProperty(
"Synaptic decay times must be defined before rise times." );
"The reversal potential array size must match the number of ports." );
}
if ( tau_tmp.size() != taus_decay.size() )
E_rev = E_rev_tmp;
}
std::vector< double > tau_tmp;
if ( updateValue< std::vector< double > >( d, names::taus_decay, tau_tmp ) )
{
if ( tau_tmp.size() != num_of_receptors )
{
throw BadProperty(
"The number of ports for synaptic rise times must be the same "
"as that of decay times." );
"The synaptic decay time array size must match the number of ports." );
}
for ( size_t i = 0; i < tau_tmp.size(); ++i )
@@ -215,10 +224,25 @@ aeif_cond_beta_multisynapse::Parameters_::set( const DictionaryDatum& d )
throw BadProperty(
"All synaptic time constants must be strictly positive" );
}
if ( tau_tmp[ i ] > taus_decay[ i ] )
}
taus_decay = tau_tmp;
}
tau_tmp.clear();
if ( updateValue< std::vector< double > >( d, names::taus_rise, tau_tmp ) )
{
if ( tau_tmp.size() != num_of_receptors )
{
throw BadProperty(
"The synaptic rise time array size must match the number of ports." );
}
for ( size_t i = 0; i < tau_tmp.size(); ++i )
{
if ( tau_tmp[ i ] <= 0 )
{
throw BadProperty(
"Synaptic rise time must be smaller than or equal to decay time." );
"All synaptic time constants must be strictly positive" );
}
}
taus_rise = tau_tmp;
@@ -430,18 +454,25 @@ aeif_cond_beta_multisynapse::init_buffers_()
void
aeif_cond_beta_multisynapse::calibrate()
{
//for ( size_t i = 0; i < P_.num_of_receptors; i++ )
//{
// if ( P_.taus_decay[ i ] < P_.taus_rise[ i ] ) {
// throw BadProperty(
// "Synaptic rise time must be smaller than or equal to decay time." );
// }
//}
// ensures initialization in case mm connected after Simulate
B_.logger_.init();
P_.receptor_types_.resize( P_.num_of_receptors_ );
for ( size_t i = 0; i < P_.num_of_receptors_; i++ )
P_.receptor_types.resize( P_.num_of_receptors );
for ( size_t i = 0; i < P_.num_of_receptors; i++ )
{
P_.receptor_types_[ i ] = i + 1;
P_.receptor_types[ i ] = i + 1;
}
V_.g0_.resize( P_.num_of_receptors_ );
V_.g0_.resize( P_.num_of_receptors );
for ( size_t i = 0; i < P_.num_of_receptors_; ++i )
for ( size_t i = 0; i < P_.num_of_receptors; ++i )
{
// the denominator (denom1) that appears in the expression of the peak time
// is computed here to check that it is != 0
@@ -477,10 +508,10 @@ aeif_cond_beta_multisynapse::calibrate()
assert( V_.RefractoryCounts_
>= 0 ); // since t_ref_ >= 0, this can only fail in error
B_.spikes_.resize( P_.num_of_receptors_ );
B_.spikes_.resize( P_.num_of_receptors );
S_.y_.resize( State_::NUMBER_OF_FIXED_STATES_ELEMENTS
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR
* P_.num_of_receptors_ ) );
* P_.num_of_receptors ) );
// reallocate instance of stepping function for ODE GSL solver
if ( B_.s_ != 0 )
@@ -575,7 +606,7 @@ aeif_cond_beta_multisynapse::update( Time const& origin,
}
}
for ( size_t i = 0; i < P_.num_of_receptors_; ++i )
for ( size_t i = 0; i < P_.num_of_receptors; ++i )
{
S_.y_[ State_::DG
+ ( State_::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR * i ) ] +=
@@ -595,7 +626,7 @@ aeif_cond_beta_multisynapse::handles_test_event( SpikeEvent&,
rport receptor_type )
{
if ( receptor_type <= 0
|| receptor_type > static_cast< port >( P_.num_of_receptors_ ) )
|| receptor_type > static_cast< port >( P_.num_of_receptors ) )
{
throw IncompatibleReceptorType( receptor_type, get_name(), "SpikeEvent" );
}
@@ -608,7 +639,7 @@ aeif_cond_beta_multisynapse::handle( SpikeEvent& e )
{
assert( e.get_delay() > 0 );
assert( ( e.get_rport() > 0 )
&& ( ( size_t ) e.get_rport() <= P_.num_of_receptors_ ) );
&& ( ( size_t ) e.get_rport() <= P_.num_of_receptors ) );
B_.spikes_[ e.get_rport() - 1 ].add_value(
e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ),
@@ -661,19 +692,11 @@ aeif_cond_beta_multisynapse_dynamics( double,
double I_syn = 0.0;

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

I would remove this empty line to tie the initialization more directly to the loop.

@heplesser

heplesser Oct 13, 2016

Contributor

I would remove this empty line to tie the initialization more directly to the loop.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

done

@golosio

golosio Oct 14, 2016

Contributor

done

for ( size_t i = 0; i < ( node.P_.num_of_receptors_
* S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR );
i += S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR )
for ( size_t i = 0; i < node.P_.num_of_receptors; ++i )
{
double g = y[ S::G + i ];
if ( g > 0 )
{
I_syn += g * ( node.P_.E_ex - V ); // g>0, E_ex>V => I_syn increases
}
else
{
I_syn += g * ( V - node.P_.E_in ); // g<0, V>E_in => I_syn decreases
}
size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

Add const.

@heplesser

heplesser Oct 13, 2016

Contributor

Add const.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

done

@golosio

golosio Oct 14, 2016

Contributor

done

// I_syn = I_syn + g*(E_rev - V)

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

Maybe move this comment above the I_syn = 0 initialization and make the summation explicit: I_syn = - sum_k g_k (V - E_rev_k).

@heplesser

heplesser Oct 13, 2016

Contributor

Maybe move this comment above the I_syn = 0 initialization and make the summation explicit: I_syn = - sum_k g_k (V - E_rev_k).

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

done

@golosio

golosio Oct 14, 2016

Contributor

done

I_syn += y[ S::G + j ] * ( node.P_.E_rev[ i ] - V );
}
// We pre-compute the argument of the exponential
@@ -693,7 +716,7 @@ aeif_cond_beta_multisynapse_dynamics( double,
// Adaptation current w.
f[ S::W ] = ( node.P_.a * ( V - node.P_.E_L ) - w ) / node.P_.tau_w;
for ( size_t i = 0; i < node.P_.num_of_receptors_; ++i )
for ( size_t i = 0; i < node.P_.num_of_receptors; ++i )
{
size_t j = i * S::NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR;

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

const

@heplesser

heplesser Oct 13, 2016

Contributor

const

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

done

@golosio

golosio Oct 14, 2016

Contributor

done

// Synaptic conductance derivative dG/dt
@@ -42,7 +42,7 @@
#include "universal_data_logger.h"
/* BeginDocumentation
Name: aeif_cond_beta_multisynapse - Conductance based exponential
Name: aeif_cond_beta_multisynapse - Conductance based adaptive exponential
integrate-and-fire neuron model according
to Brette and Gerstner (2005) with
multiple synaptic rise time and decay
@@ -51,39 +51,41 @@
Description:
aeif_cond_beta_multisynapse is an extension of
aeif_cond_alpha_multisynapse. It allows an arbitrary number of synaptic
rise time and decay time constant. Synaptic conductance is modeled by a
aeif_cond_beta_multisynapse is a conductance based adaptive exponential
integrate-and-fire neuron model. It allows an arbitrary number of synaptic
rise time and decay time constants. Synaptic conductance is modeled by a
beta function, as described by A. Roth and M.C.W. van Rossum
in Computational Modeling Methods for Neuroscientists, MIT Press 2013,
Chapter 6

This comment has been minimized.

@heplesser

heplesser Oct 26, 2016

Contributor

Fullstop missing at end of sentence.

@heplesser

heplesser Oct 26, 2016

Contributor

Fullstop missing at end of sentence.

This comment has been minimized.

@golosio

golosio Oct 26, 2016

Contributor

done.

@golosio

golosio Oct 26, 2016

Contributor

done.

The time constants are supplied by two arrays taus_rise and taus_decay for
the synaptic rise time and decay time, respectively. Port numbers
are then automatically assigned and there range is from 1 to n. (n
being the index of the last element of the tau_ex and tau_in
arrays).
The number of receptor ports is supplied by the variable "num_of_receptors".
The time constants are supplied by two arrays, "taus_rise" and "taus_decay" for
the synaptic rise time and decay time, respectively. The synaptic
reversal potentials are supplied by the array "E_rev". The port numbers
are automatically assigned in the range from 1 to num_of_receptors.
During connection, the ports are selected with the property "receptor_type".
Examples:
% PyNEST example, of how to assign synaptic rise time and decay time
% to a receptor type.

This comment has been minimized.

@heplesser

heplesser Oct 26, 2016

Contributor

I think the comment here is not needed. If you want to keep it, mark the comment with Python's # instead of %, which introduces comments in SLI.

@heplesser

heplesser Oct 26, 2016

Contributor

I think the comment here is not needed. If you want to keep it, mark the comment with Python's # instead of %, which introduces comments in SLI.

This comment has been minimized.

@golosio

golosio Oct 26, 2016

Contributor

removed.

@golosio

golosio Oct 26, 2016

Contributor

removed.

import nest
import nest
import numpy as np
neuron = nest.Create('aeif_cond_beta_multisynapse')
nest.SetStatus(neuron, {"V_peak": 0.0, "a": 4.0, "b":80.5})
nest.SetStatus(neuron, {'taus_decay':[50.0,20.0,20.0,20.0],
nest.SetStatus(neuron, {'num_of_receptors':4,

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

num_of_receptors should never be set explicitly.

@heplesser

heplesser Oct 13, 2016

Contributor

num_of_receptors should never be set explicitly.

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

removed

@golosio

golosio Oct 14, 2016

Contributor

removed

'E_rev':[0.0,0.0,0.0,-85.0],
'taus_decay':[50.0,20.0,20.0,20.0],
'taus_rise':[10.0,10.0,1.0,1.0]})
spike = nest.Create('spike_generator', params = {'spike_times':
np.array([10.0])})
np.array([10.0])})
voltmeter = nest.Create('voltmeter', 1, {'withgid': True})
delays=[1.0, 300.0, 500.0, 700.0]
w=[1.0, 1.0, 1.0, -1.0]
w=[1.0, 1.0, 1.0, 1.0]
for syn in range(4):
nest.Connect(spike, neuron, syn_spec={'model': 'static_synapse',
'weight': w[syn],

This comment has been minimized.

@heplesser

heplesser Oct 13, 2016

Contributor

Don't you need to specify receptor_type here?

@heplesser

heplesser Oct 13, 2016

Contributor

Don't you need to specify receptor_type here?

This comment has been minimized.

@golosio

golosio Oct 14, 2016

Contributor

You mean that the order should be important? I do not know, I did not see any difference in the simulations if I put 'receptor_type' at the end of the syn_spec list, in any case now I've put it in the second position after 'model'.

@golosio

golosio Oct 14, 2016

Contributor

You mean that the order should be important? I do not know, I did not see any difference in the simulations if I put 'receptor_type' at the end of the syn_spec list, in any case now I've put it in the second position after 'model'.

@@ -96,15 +98,15 @@
Vms = dmm["events"]["V_m"]
ts = dmm["events"]["times"]
import pylab
pylab.figure(1)
pylab.figure(2)
pylab.plot(ts, Vms)
pylab.show()
Sends: SpikeEvent
Receives: SpikeEvent, CurrentEvent, DataLoggingRequest
Adapted by Bruno Golosio from aeif_cond_alpha_multisynapse
author: Bruno Golosio 07/10/2016
SeeAlso: aeif_cond_alpha_multisynapse
*/
@@ -182,27 +184,26 @@ class aeif_cond_beta_multisynapse : public Archiving_Node
double g_L; //!< Leak Conductance in nS
double C_m; //!< Membrane Capacitance in pF
double E_ex; //!< Excitatory reversal Potential in mV
double E_in; //!< Inhibitory reversal Potential in mV
double E_L; //!< Leak reversal Potential (aka resting potential) in mV
double Delta_T; //!< Slope faktor in ms.
double tau_w; //!< adaptation time-constant in ms.
double a; //!< Subthreshold adaptation in nS.
double b; //!< Spike-triggered adaptation in pA
double V_th; //!< Spike threshold in mV.
double t_ref; //!< Refractory period in ms.
size_t num_of_receptors;
std::vector< double > taus_rise; //!< Rise time of synaptic conductance
//!< in ms..
std::vector< double > taus_decay; //!< Decay time of synaptic conductance
//!< in ms..
std::vector< double > E_rev; //!< reversal potentials in mV
double I_e; //!< Intrinsic current in pA.
double gsl_error_tol; //!< error bound for GSL integrator
// type is long because other types are not put through in GetStatus
std::vector< long > receptor_types_;
size_t num_of_receptors_;
std::vector< long > receptor_types;
// boolean flag which indicates whether the neuron has connections
bool has_connections_;
Oops, something went wrong.
ProTip! Use n and p to navigate between commits in a pull request.