diff --git a/Cortical_Column.cpp b/Cortical_Column.cpp index 89467e1..f78eea4 100644 --- a/Cortical_Column.cpp +++ b/Cortical_Column.cpp @@ -29,189 +29,146 @@ * to auditory stimulation. * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, * JC Claussen. - * PLoS Computational Biology (in review). + * PLoS Computational Biology http://dx.doi.org/10.1371/journal.pcbi.1005022 */ -/****************************************************************************************************/ -/* Functions of the cortical module */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Functions of the cortical module */ +/******************************************************************************/ #include "Cortical_Column.h" -/****************************************************************************************************/ -/* Initialization of RNG */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Initialization of RNG */ +/******************************************************************************/ void Cortical_Column::set_RNG(void) { - extern const double dt; - /* Number of independent random variables */ - int N = 2; - - /* Create RNG for each stream */ - for (int i=0; iy[N] - s_ep[N]) - 2 * gamma_e * x_ep[N]) + noise_xRK(N, 0); - x_ei[N+1] = x_ei[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ip * get_Qp(N) + N_it * Thalamus->y[N] - s_ei[N]) - 2 * gamma_e * x_ei[N]) + noise_xRK(N, 1) ; - x_gp[N+1] = x_gp[0] + A[N] * dt*(pow(gamma_g, 2) * (N_pi * get_Qi(N) - s_gp[N]) - 2 * gamma_g * x_gp[N]); - x_gi[N+1] = x_gi[0] + A[N] * dt*(pow(gamma_g, 2) * (N_ii * get_Qi(N) - s_gi[N]) - 2 * gamma_g * x_gi[N]); - x [N+1] = x [0] + A[N] * dt*(pow(nu, 2) * ( get_Qp(N) - y [N]) - 2 * nu * x [N]); + extern const double dt; + Vp [N+1] = Vp [0] + A[N] * dt*(-(I_L_p(N) + I_ep(N) + I_gp(N))/tau_p - I_KNa(N)); + Vi [N+1] = Vi [0] + A[N] * dt*(-(I_L_i(N) + I_ei(N) + I_gi(N))/tau_i); + Na [N+1] = Na [0] + A[N] * dt*(alpha_Na * get_Qp(N) - Na_pump(N))/tau_Na; + s_ep[N+1] = s_ep[0] + A[N] * dt*(x_ep[N]); + s_ei[N+1] = s_ei[0] + A[N] * dt*(x_ei[N]); + s_gp[N+1] = s_gp[0] + A[N] * dt*(x_gp[N]); + s_gi[N+1] = s_gi[0] + A[N] * dt*(x_gi[N]); + y [N+1] = y [0] + A[N] * dt*(x [N]); + x_ep[N+1] = x_ep[0] + A[N] * dt*(gamma_e*gamma_e * (N_pp * get_Qp(N) + N_pt * Thalamus->y[N] - s_ep[N]) - 2 * gamma_e * x_ep[N]) + noise_xRK(N, 0); + x_ei[N+1] = x_ei[0] + A[N] * dt*(gamma_e*gamma_e * (N_ip * get_Qp(N) + N_it * Thalamus->y[N] - s_ei[N]) - 2 * gamma_e * x_ei[N]) + noise_xRK(N, 1) ; + x_gp[N+1] = x_gp[0] + A[N] * dt*(gamma_g*gamma_g * (N_pi * get_Qi(N) - s_gp[N]) - 2 * gamma_g * x_gp[N]); + x_gi[N+1] = x_gi[0] + A[N] * dt*(gamma_g*gamma_g * (N_ii * get_Qi(N) - s_gi[N]) - 2 * gamma_g * x_gi[N]); + x [N+1] = x [0] + A[N] * dt*(nu * nu * ( get_Qp(N) - y [N]) - 2 * nu * x [N]); } -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - -/****************************************************************************************************/ -/* Function that adds all SRK terms */ -/****************************************************************************************************/ void Cortical_Column::add_RK(void) { - Vp [0] = (-3*Vp [0] + 2*Vp [1] + 4*Vp [2] + 2*Vp [3] + Vp [4])/6; - Vi [0] = (-3*Vi [0] + 2*Vi [1] + 4*Vi [2] + 2*Vi [3] + Vi [4])/6; - Na [0] = (-3*Na [0] + 2*Na [1] + 4*Na [2] + 2*Na [3] + Na [4])/6; - s_ep[0] = (-3*s_ep[0] + 2*s_ep[1] + 4*s_ep[2] + 2*s_ep[3] + s_ep[4])/6; - s_ei[0] = (-3*s_ei[0] + 2*s_ei[1] + 4*s_ei[2] + 2*s_ei[3] + s_ei[4])/6; - s_gp[0] = (-3*s_gp[0] + 2*s_gp[1] + 4*s_gp[2] + 2*s_gp[3] + s_gp[4])/6; - s_gi[0] = (-3*s_gi[0] + 2*s_gi[1] + 4*s_gi[2] + 2*s_gi[3] + s_gi[4])/6; - y [0] = (-3*y [0] + 2*y [1] + 4*y [2] + 2*y [3] + y [4])/6; - x_ep[0] = (-3*x_ep[0] + 2*x_ep[1] + 4*x_ep[2] + 2*x_ep[3] + x_ep[4])/6 + noise_aRK(0); - x_ei[0] = (-3*x_ei[0] + 2*x_ei[1] + 4*x_ei[2] + 2*x_ei[3] + x_ei[4])/6 + noise_aRK(1); - x_gp[0] = (-3*x_gp[0] + 2*x_gp[1] + 4*x_gp[2] + 2*x_gp[3] + x_gp[4])/6; - x_gi[0] = (-3*x_gi[0] + 2*x_gi[1] + 4*x_gi[2] + 2*x_gi[3] + x_gi[4])/6; - x [0] = (-3*x [0] + 2*x [1] + 4*x [2] + 2*x [3] + x [4])/6; - - /* Generate noise for the next iteration */ - for (unsigned i=0; i #include + #include "Random_Stream.h" #include "Thalamic_Column.h" -using std::vector; class Thalamic_Column; -/****************************************************************************************************/ -/* Macro for vector initialization */ -/****************************************************************************************************/ -#ifndef _INIT -#define _INIT(x) {x, 0.0, 0.0, 0.0, 0.0} -#endif -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Implementation of the cortical module */ -/****************************************************************************************************/ class Cortical_Column { public: - /* Constructors */ - Cortical_Column(void) - {set_RNG();} - - Cortical_Column(double* Param, double* Con) - :sigma_p (Param[0]), g_KNa (Param[1]), dphi (Param[2]), - N_pt (Con[2]), N_it (Con[3]) - {set_RNG();} + Cortical_Column(double* Param, double* Con) + :sigma_p (Param[0]), g_KNa (Param[1]), dphi (Param[2]), + N_pt (Con[2]), N_it (Con[3]) + {set_RNG();} - /* Connect to the thalamic module */ - void get_Thalamus(Thalamic_Column& T) {Thalamus = &T;} + /* Connect to the thalamic module */ + void get_Thalamus(Thalamic_Column& T) {Thalamus = &T;} - /* Data storage access */ - friend void get_data (int, Cortical_Column&, Thalamic_Column&, vector); - - /* ODE functions */ - void set_RK (int); - void add_RK (void); - - /* Stimulation protocol access */ - friend class Stim; - friend class Thalamic_Column; + /* ODE functions */ + void set_RK (int); + void add_RK (void); private: - /* Declaration of private functions */ - /* Initialize the RNGs */ - void set_RNG (void); - - /* Firing rates */ - double get_Qp (int) const; - double get_Qi (int) const; - - /* Currents */ - double I_ep (int) const; - double I_ei (int) const; - double I_gp (int) const; - double I_gi (int) const; - double I_L_p (int) const; - double I_L_g (int) const; - double I_KNa (int) const; - - /* Potassium pump */ - double Na_pump (int) const; - - /* Noise functions */ - double noise_xRK (int,int) const; - double noise_aRK (int) const; - - /* Declaration and Initialization of parameters */ - /* Membrane time in ms */ - const double tau_p = 30; - const double tau_i = 30; - - /* Maximum firing rate in ms^-1 */ - const double Qp_max = 30.E-3; - const double Qi_max = 60.E-3; - - /* Sigmoid threshold in mV */ - const double theta_p = -58.5; - const double theta_i = -58.5; - - /* Sigmoid gain in mV */ - const double sigma_p = 4; - const double sigma_i = 6; - - /* Scaling parameter for sigmoidal mapping (dimensionless) */ - const double C1 = (M_PI/sqrt(3)); - - /* parameters of the firing adaption */ - const double alpha_Na = 2; /* Sodium influx per spike in mM ms */ - const double tau_Na = 1.7; /* Sodium time constant in ms */ - - const double R_pump = 0.09; /* Na-K pump constant in mM/ms */ - const double Na_eq = 9.5; /* Na-eq concentration in mM */ - - /* PSP rise time in ms^-1 */ - const double gamma_e = 70E-3; - const double gamma_g = 58.6E-3; - - /* Axonal flux time constant */ - const double nu = 120E-3; - - /* Leak weight in aU*/ - const double g_L = 1.; - - /* Synaptic weight in ms */ - const double g_AMPA = 1.; - const double g_GABA = 1.; - - /* Conductivity */ - /* KNa in mS/cm^2 */ - const double g_KNa = 1.33; - - /* Reversal potentials in mV */ - /* Synaptic */ - const double E_AMPA = 0; - const double E_GABA = -70; - - /* Leak */ - const double E_L_p = -64; - const double E_L_g = -64; - - /* Potassium */ - const double E_K = -100; - - /* Noise parameters in ms^-1 */ - const double mphi = 0E-3; - const double dphi = 20E-1; - double input = 0.0; - - /* Connectivities (dimensionless) */ - const double N_pp = 115; - const double N_ip = 72; - const double N_pi = 90; - const double N_ii = 90; - const double N_pt = 2.5; - const double N_it = 2.5; - - /* Pointer to thalamic column */ - Thalamic_Column* Thalamus; - - /* Parameters for SRK4 iteration */ - const vector A = {0.5, 0.5, 1.0, 1.0}; - const vector B = {0.75, 0.75, 0.0, 0.0}; - - /* Random number generators */ - vector MTRands; - - /* Container for noise */ - vector Rand_vars; - - /* Population variables */ - vector Vp = _INIT(E_L_p), /* excitatory membrane voltage */ - Vi = _INIT(E_L_g), /* inhibitory membrane voltage */ - Na = _INIT(Na_eq), /* Na concentration */ - s_ep = _INIT(0.0), /* PostSP from excitatory to excitatory population */ - s_ei = _INIT(0.0), /* PostSP from excitatory to inhibitory population */ - s_gp = _INIT(0.0), /* PostSP from inhibitory to excitatory population */ - s_gi = _INIT(0.0), /* PostSP from inhibitory to inhibitory population */ - y = _INIT(0.0), /* axonal flux */ - x_ep = _INIT(0.0), /* derivative of s_ep */ - x_ei = _INIT(0.0), /* derivative of s_ei */ - x_gp = _INIT(0.0), /* derivative of s_gp */ - x_gi = _INIT(0.0), /* derivative of s_gi */ - x = _INIT(0.0); /* derivative of y */ -}; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ + /* Declaration of private functions */ + /* Initialize the RNGs */ + void set_RNG (void); + + /* Firing rates */ + double get_Qp (int) const; + double get_Qi (int) const; + + /* Currents */ + double I_ep (int) const; + double I_ei (int) const; + double I_gp (int) const; + double I_gi (int) const; + double I_L_p (int) const; + double I_L_i (int) const; + double I_KNa (int) const; + /* Potassium pump */ + double Na_pump (int) const; + + /* Noise functions */ + double noise_xRK (int,int) const; + double noise_aRK (int) const; + + /* Helper functions */ + inline std::vector init (double value) + {return {value, 0.0, 0.0, 0.0, 0.0};} + + inline void add_RK (std::vector& var) + {var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6;} + + inline void add_RK_noise (std::vector& var, unsigned noise) + {var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6 + noise_aRK(noise);} + + /* Declaration and Initialization of parameters */ + /* Membrane time in ms */ + const double tau_p = 30; + const double tau_i = 30; + + /* Maximum firing rate in ms^-1 */ + const double Qp_max = 30.E-3; + const double Qi_max = 60.E-3; + + /* Sigmoid threshold in mV */ + const double theta_p = -58.5; + const double theta_i = -58.5; + + /* Sigmoid gain in mV */ + const double sigma_p = 4; + const double sigma_i = 6; + + /* Scaling parameter for sigmoidal mapping (dimensionless) */ + const double C1 = (M_PI/sqrt(3)); + + /* parameters of the firing adaption */ + const double alpha_Na = 2; /* Sodium influx per spike in mM ms */ + const double tau_Na = 1.7; /* Sodium time constant in ms */ + + const double R_pump = 0.09; /* Na-K pump constant in mM/ms */ + const double Na_eq = 9.5; /* Na-eq concentration in mM */ + + /* PSP rise time in ms^-1 */ + const double gamma_e = 70E-3; + const double gamma_g = 58.6E-3; + + /* Axonal flux time constant */ + const double nu = 120E-3; + + /* Leak weight in aU*/ + const double g_L = 1.; + + /* Synaptic weight in ms */ + const double g_AMPA = 1.; + const double g_GABA = 1.; + + /* Conductivity */ + /* KNa in mS/cm^2 */ + const double g_KNa = 1.33; + + /* Reversal potentials in mV */ + /* Synaptic */ + const double E_AMPA = 0; + const double E_GABA = -70; + + /* Leak */ + const double E_L_p = -64; + const double E_L_i = -64; + + /* Potassium */ + const double E_K = -100; + + /* Noise parameters in ms^-1 */ + const double mphi = 0E-3; + const double dphi = 20E-1; + double input = 0.0; + + /* Connectivities (dimensionless) */ + const double N_pp = 115; + const double N_ip = 72; + const double N_pi = 90; + const double N_ii = 90; + const double N_pt = 2.5; + const double N_it = 2.5; + + /* Pointer to thalamic column */ + Thalamic_Column* Thalamus; + + /* Parameters for SRK4 iteration */ + const std::vector A = {0.5, 0.5, 1.0, 1.0}; + const std::vector B = {0.75, 0.75, 0.0, 0.0}; + + /* Random number generators */ + std::vector MTRands; + + /* Container for noise */ + std::vector Rand_vars; + + /* Population variables */ + std::vector Vp = init(E_L_p), /* excitatory membrane voltage */ + Vi = init(E_L_i), /* inhibitory membrane voltage */ + Na = init(Na_eq), /* Na concentration */ + s_ep= init(0.0), /* PostSP from excitatory to excitatory population */ + s_ei= init(0.0), /* PostSP from excitatory to inhibitory population */ + s_gp= init(0.0), /* PostSP from inhibitory to excitatory population */ + s_gi= init(0.0), /* PostSP from inhibitory to inhibitory population */ + y = init(0.0), /* axonal flux */ + x_ep= init(0.0), /* derivative of s_ep */ + x_ei= init(0.0), /* derivative of s_ei */ + x_gp= init(0.0), /* derivative of s_gp */ + x_gi= init(0.0), /* derivative of s_gi */ + x = init(0.0); /* derivative of y */ + + /* Data storage access */ + friend void get_data (int, Cortical_Column&, Thalamic_Column&, std::vector); + + /* Stimulation protocol access */ + friend class Stim; + friend class Thalamic_Column; +}; diff --git a/Data_Storage.h b/Data_Storage.h index 6d2caa3..86b7bd2 100644 --- a/Data_Storage.h +++ b/Data_Storage.h @@ -20,30 +20,20 @@ * THE SOFTWARE. * * AUTHORS: Michael Schellenberger Costa: mschellenbergercosta@gmail.com - * - * Based on: A thalamocortical neural mass model of the EEG during NREM sleep and its response - * to auditory stimulation. - * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, - * JC Claussen. - * PLoS Computational Biology (in review). */ -/****************************************************************************************************/ -/* Functions for data storage */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Functions for data storage */ +/******************************************************************************/ #pragma once +#include #include "Cortical_Column.h" #include "Thalamic_Column.h" -/****************************************************************************************************/ -/* Save data */ -/****************************************************************************************************/ -void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, vector Data) { - Data[0][counter] = Cortex.Vp [0]; - Data[1][counter] = Thalamus.Vt [0]; - Data[2][counter] = Thalamus.Ca [0]; - Data[3][counter] = Thalamus.act_h (); +void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, + std::vector pData) { + pData[0][counter] = Cortex.Vp [0]; + pData[1][counter] = Thalamus.Vt [0]; + pData[2][counter] = Thalamus.Ca [0]; + pData[3][counter] = Thalamus.act_h (); } -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/NM_TC.pro b/NM_TC.pro index 0eeabed..f174150 100644 --- a/NM_TC.pro +++ b/NM_TC.pro @@ -6,17 +6,21 @@ CONFIG -= qt TARGET = TC.cpp SOURCES += Cortical_Column.cpp \ - TC.cpp \ - TC_mex.cpp \ - Thalamic_Column.cpp + TC.cpp \ + TC_mex.cpp \ + Thalamic_Column.cpp HEADERS += Cortical_Column.h \ - Data_Storage.h \ - ODE.h \ - Random_Stream.h \ - Stimulation.h \ - Thalamic_Column.h - -QMAKE_CXXFLAGS += -std=c++11 -O3 + Data_Storage.h \ + ODE.h \ + Random_Stream.h \ + Stimulation.h \ + Thalamic_Column.h SOURCES -= TC_mex.cpp + +QMAKE_CXXFLAGS += -std=c++11 +QMAKE_CXXFLAGS_RELEASE -= -O1 +QMAKE_CXXFLAGS_RELEASE -= -O2 +QMAKE_CXXFLAGS_RELEASE *= -O3 + diff --git a/ODE.h b/ODE.h index d8d81c3..9e5305f 100644 --- a/ODE.h +++ b/ODE.h @@ -25,30 +25,24 @@ * to auditory stimulation. * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, * JC Claussen. - * PLoS Computational Biology (in review). + * PLoS Computational Biology http://dx.doi.org/10.1371/journal.pcbi.1005022 */ -/****************************************************************************************************/ -/* Implementation of the ODE solver */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Functions for SRK iteration */ +/******************************************************************************/ #pragma once #include "Cortical_Column.h" #include "Thalamic_Column.h" -/****************************************************************************************************/ -/* Evaluation of SRK4 */ -/****************************************************************************************************/ void ODE(Cortical_Column& Cortex, Thalamic_Column& Thalamus) { - /* First calculate every ith RK moment. Has to be in order, 1th moment first */ - for (int i=0; i<4; ++i) { - Cortex.set_RK(i); - Thalamus.set_RK(i); - } + /* First calculate every ith RK moment. Has to be in order, 1th moment first */ + for (unsigned i=0; i<4; ++i) { + Cortex.set_RK(i); + Thalamus.set_RK(i); + } - /* Add all moments */ - Cortex.add_RK(); - Thalamus.add_RK(); + /* Add all moments */ + Cortex.add_RK(); + Thalamus.add_RK(); } -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/Random_Stream.h b/Random_Stream.h index f5fb19a..6695d38 100644 --- a/Random_Stream.h +++ b/Random_Stream.h @@ -23,58 +23,34 @@ * Stefanie Gareis: gareis@inb.uni-luebeck.de */ -/****************************************************************************************************/ -/* Random number streams */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Random number streams */ +/******************************************************************************/ #pragma once #include -/****************************************************************************************************/ -/* Struct for normal distribution */ -/****************************************************************************************************/ -struct random_stream_normal -{ - /* Random number engine: Mersenne-Twister */ +class random_stream_normal { +public: + explicit random_stream_normal(double mean, double stddev) + : mt(rand()), norm_dist(mean, stddev) {} + explicit random_stream_normal(double mean, double stddev, double seed) + : mt(seed), norm_dist(mean, stddev) {} + + double operator ()(void) { return norm_dist(mt); } +private: std::mt19937_64 mt; - /* Random number generator: Normal-distribution */ std::normal_distribution norm_dist; - - /* Constructors */ - random_stream_normal(){} - random_stream_normal(double mean, double stddev) - : mt(rand()) , norm_dist(mean, stddev) - {} - - /* Overwrites the function-call operator "( )" */ - double operator( )(void) { - return norm_dist(mt); - } }; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ -/****************************************************************************************************/ -/* Struct for uniform int distribution */ -/****************************************************************************************************/ -struct random_stream_uniform_int -{ - /* Random number engine: Mersenne-Twister */ +class random_stream_uniform_int { +public: + explicit random_stream_uniform_int(int lower_bound, int upper_bound) + : mt(rand()), uniform_dist(lower_bound, upper_bound) {} + explicit random_stream_uniform_int(int lower_bound, int upper_bound, double seed) + : mt(seed), uniform_dist(lower_bound, upper_bound) {} + + double operator ()(void) { return uniform_dist(mt); } +private: std::mt19937_64 mt; - /* Random number generator: Uniform integer-distribution */ std::uniform_int_distribution<> uniform_dist; - - /* Constructors */ - random_stream_uniform_int(){} - random_stream_uniform_int(double lower_bound, double upper_bound) - : mt(rand()) , uniform_dist(lower_bound, upper_bound) - {} - - /* Overwrites the function-call operator "( )" */ - double operator( )(void) { - return uniform_dist(mt); - } }; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/Stimulation.h b/Stimulation.h index aa1ef64..55b6523 100644 --- a/Stimulation.h +++ b/Stimulation.h @@ -25,330 +25,307 @@ * to auditory stimulation. * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, * JC Claussen. - * PLoS Computational Biology (in review). + * PLoS Computational Biology http://dx.doi.org/10.1371/journal.pcbi.1005022 */ -/****************************************************************************************************/ -/* Implementation of the stimulation protocol */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Implementation of the stimulation protocol */ +/******************************************************************************/ #pragma once -#include "Random_Stream.h" +#include + #include "Cortical_Column.h" +#include "Random_Stream.h" #include "Thalamic_Column.h" -/****************************************************************************************************/ -/* Stimulation object */ -/****************************************************************************************************/ + +/******************************************************************************/ +/* Stimulation object */ +/******************************************************************************/ class Stim { public: - /* Empty constructor for compiling */ - Stim(void); - - /* Constructor with references and stimulation variables */ - Stim(Cortical_Column& C, Thalamic_Column& T, double* var) - { Cortex = &C; - Thalamus = &T; - setup(var);} - - /* Initialize stimulation class with respect to stimulation mode */ - void setup (double* var_stim); + /* Constructor with references and stimulation variables */ + Stim(Cortical_Column& C, Thalamic_Column& T, double* var) + { Cortex = &C; + Thalamus = &T; + setup(var);} - /* Check whether stimulation should be started/stopped */ - void check_stim (int time); - - /* Create MATLAB container for marker storage */ - mxArray* get_marker(void); + /* Initialize stimulation class with respect to stimulation mode */ + void setup (double* var_stim); + /* Check whether stimulation should be started/stopped */ + void check_stim (int time); private: - /* Mode of stimulation */ - /* 0 == none */ - /* 1 == semi-periodic */ - /* 2 == phase dependent */ - int mode = 0; + /* Mode of stimulation */ + /* 0 == none */ + /* 1 == semi-periodic */ + /* 2 == phase dependent */ + int mode = 0; - /* Default values already in dt: E1==ms, E4==s */ - /* Stimulation strength */ - double strength = 0.0; + /* Default values already in dt: E1==ms, E4==s */ + /* Stimulation strength */ + double strength = 0.0; - /* Duration of the stimulation */ - int duration = 120E1; + /* Duration of the stimulation */ + int duration = 120E1; - /* Interval between different Stimulus events */ - int ISI = 5E4; + /* Interval between different Stimulus events */ + int ISI = 5E4; - /* Range of Inter Stimulus Interval */ - int ISI_range = 1E4; + /* Range of Inter Stimulus Interval */ + int ISI_range = 1E4; - /* Number of stimuli in case of multiple stimuli per event */ - int number_of_stimuli = 1; + /* Number of stimuli in case of multiple stimuli per event */ + int number_of_stimuli = 1; - /* Time until next stimulus */ - /* Function varies between different stimulation modes */ - int time_to_stimuli = 350E1; + /* Time until next stimulus */ + /* Function varies between different stimulation modes */ + int time_to_stimuli = 350E1; - /* Time between stimuli in case of multiple stimuli per event */ - int time_between_stimuli = 1050E1; + /* Time between stimuli in case of multiple stimuli per event */ + int time_between_stimuli = 1050E1; - /* Threshold for phase dependent stimulation */ - double threshold = -68; + /* Threshold for phase dependent stimulation */ + double threshold = -68; - /* Internal variables */ - /* Simulation on for TRUE and off for FALSE */ - bool stimulation_started = false; + /* Internal variables */ + /* Simulation on for TRUE and off for FALSE */ + bool stimulation_started = false; - /* If threshold has been reached */ - bool threshold_crossed = false; + /* If threshold has been reached */ + bool threshold_crossed = false; - /* If minimum was found */ - bool minimum_found = false; + /* If minimum was found */ + bool minimum_found = false; - /* If a stimulation event has occurred and there is a minimal time (pause) until the next one */ - bool stimulation_paused = false; + /* If a stimulation event has occurred and there is a minimal time (pause) until the next one */ + bool stimulation_paused = false; - /* If burst mode is enabled */ - bool burst_enabled = false; + /* If burst mode is enabled */ + bool burst_enabled = false; - /* In case of bursted stimulation start the bursts */ - bool burst_started = true; + /* In case of bursted stimulation start the bursts */ + bool burst_started = true; - /* Length of a burst stimulus */ - int burst_length = 10; + /* Length of a burst stimulus */ + int burst_length = 10; - /* Length of silence between burst stimuli */ - int burst_ISI = 10; + /* Length of silence between burst stimuli */ + int burst_ISI = 10; - /* Onset in time steps where no data is recorded, so stimulation has to be delayed */ - int onset_correction = 10E4; + /* Onset in time steps where no data is recorded, so stimulation has to be delayed */ + int onset_correction = 10E4; - /* Counter for number of stimuli that occurred within a stimulation event */ - int count_stimuli = 1; + /* Counter for number of stimuli that occurred within a stimulation event */ + int count_stimuli = 1; - /* Counter for number of bursts per stimuli */ - int count_bursts = 0; + /* Counter for number of bursts per stimuli */ + int count_bursts = 0; - /* Counter for stimulation duration since started*/ - int count_duration = 0; + /* Counter for stimulation duration since started*/ + int count_duration = 0; - /* Counter after minimum was found */ - int count_to_start = 0; + /* Counter after minimum was found */ + int count_to_start = 0; - /* Counter for time between two stimulation events (with multiple tones) */ - int count_pause = 0; + /* Counter for time between two stimulation events (with multiple tones) */ + int count_pause = 0; - /* Old voltage value for minimum detection */ - double Vp_old = 0.0; + /* Old voltage value for minimum detection */ + double Vp_old = 0.0; - /* Pointer to columns */ - Cortical_Column* Cortex; - Thalamic_Column* Thalamus; + /* Pointer to columns */ + Cortical_Column* Cortex; + Thalamic_Column* Thalamus; - /* Data containers */ - std::vector marker_stimulation; + /* Data containers */ + std::vector marker_stimulation; - /* Random number generator in case of semi-periodic stimulation */ - random_stream_uniform_int Uniform_Distribution; -}; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ + /* Random number generator in case of semi-periodic stimulation */ + random_stream_uniform_int Uniform_Distribution = random_stream_uniform_int(0, 0); + /* Create MATLAB container for marker storage */ + friend mxArray* get_marker(Stim &stim); +}; -/****************************************************************************************************/ -/* Function definitions */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Function definitions */ +/******************************************************************************/ void Stim::setup (double* var_stim) { - extern const int onset; - extern const int res; + extern const int onset; + extern const int res; - /* Set the onset onset_correction for the marker */ - onset_correction = onset * res; + /* Set the onset onset_correction for the marker */ + onset_correction = onset * res; - /* Mode of stimulation */ - mode = (int) var_stim[0]; + /* Mode of stimulation */ + mode = (int) var_stim[0]; - /* Scale the stimulation strength from s^-1 (Hz) to ms^-1 */ - strength = var_stim[1] / 1000; + /* Scale the stimulation strength from s^-1 (Hz) to ms^-1 */ + strength = var_stim[1] / 1000; - /* Scale duration from ms to dt */ - duration = (int) var_stim[2] * res / 1000; + /* Scale duration from ms to dt */ + duration = (int) var_stim[2] * res / 1000; - /* Scale the inter stimulus event interval from s to dt */ - ISI = (int) var_stim[3] * res; + /* Scale the inter stimulus event interval from s to dt */ + ISI = (int) var_stim[3] * res; - /* Scale inter stimulus event interval range from s to dt */ - ISI_range = (int) var_stim[4] * res; + /* Scale inter stimulus event interval range from s to dt */ + ISI_range = (int) var_stim[4] * res; - /* Number of stimuli per Stimulus event */ - number_of_stimuli = (int) var_stim[5]; + /* Number of stimuli per Stimulus event */ + number_of_stimuli = (int) var_stim[5]; - /* Scale time_between_stimuli from ms to dt */ - time_between_stimuli = (int) var_stim[6] * res / 1000; + /* Scale time_between_stimuli from ms to dt */ + time_between_stimuli = (int) var_stim[6] * res / 1000; - /* Scale the length of burst_length and burst_ISI from ms to dt*/ - burst_length = (int) 2 * res / 1000; - burst_ISI = (int) 28 * res / 1000; + /* Scale the length of burst_length and burst_ISI from ms to dt*/ + burst_length = (int) 2 * res / 1000; + burst_ISI = (int) 28 * res / 1000; - /* If ISI is fixed do not create RNG */ - if (mode == 1) { - /* Set first time_to_stimuli to 1 sec after onset */ - time_to_stimuli = (int) (onset+1) * res; + /* If ISI is fixed do not create RNG */ + if (mode == 1) { + /* Set first time_to_stimuli to 1 sec after onset */ + time_to_stimuli = (int) (onset+1) * res; - /* If ISI is random create RNG */ - if (ISI_range != 0){ - /* Generate uniform distribution */ - Uniform_Distribution = random_stream_uniform_int(ISI-ISI_range, ISI+ISI_range); - } - } else { - /* In case of phase dependent stimulation, time_to_stim is the time from minimum detection to start of stimulation */ - /* Scale time_to_stimuli from ms to dt */ - time_to_stimuli = (int) var_stim[7] * res / 1000; - } + /* If ISI is random create RNG */ + if (ISI_range != 0){ + /* Generate uniform distribution */ + Uniform_Distribution = random_stream_uniform_int(ISI-ISI_range, ISI+ISI_range); + } + } else { + /* In case of phase dependent stimulation, time_to_stim is the time from minimum detection to start of stimulation */ + /* Scale time_to_stimuli from ms to dt */ + time_to_stimuli = (int) var_stim[7] * res / 1000; + } } void Stim::check_stim (int time) { - /* Check if stimulation should start */ - switch (mode) { - - /* No stimulation */ - default: - break; - - /* Semi-periodic stimulation */ - case 1: - /* Check if stimulation time is reached */ - if(time == time_to_stimuli) { - /* Switch stimulation on */ - stimulation_started = true; - Thalamus->set_input(strength); - - /* Add marker for the first stimuli in the event */ - if(count_stimuli == 1) { - marker_stimulation.push_back(time - onset_correction); - } - - /* Check if multiple stimuli should be applied */ - if (count_stimuli < number_of_stimuli) { - /* Update the timer with respect to time between stimuli */ - time_to_stimuli += time_between_stimuli; - count_stimuli++; - } - /* After last stimulus in event update the timer with respect to (random) ISI*/ - else { - time_to_stimuli += (ISI_range==0)? ISI : Uniform_Distribution(); - - /* Reset the stimulus counter for next stimulation event */ - count_stimuli = 1; - } - } - break; - - /* Phase dependent stimulation */ - case 2: - /* Search for threshold */ - if(!stimulation_started && !minimum_found && !threshold_crossed && time>onset_correction && !stimulation_paused) { - if(Cortex->Vp[0]<=threshold) { - threshold_crossed = true; - } - } - - /* Search for minimum */ - if(threshold_crossed) { - if(Cortex->Vp[0]>Vp_old) { - threshold_crossed = false; - minimum_found = true; - Vp_old = 0; - } else { - Vp_old = Cortex->Vp[0]; - } - } - - /* Wait until the stimulation should start */ - if(minimum_found) { - /* Start stimulation after time_to_stimuli has passed */ - if(count_to_start==time_to_stimuli + (count_stimuli-1) * time_between_stimuli) { - stimulation_started = true; - Thalamus->set_input(strength); - - /* Add marker for the first stimuli in the event */ - if(count_stimuli == 1) { - marker_stimulation.push_back(time - onset_correction); - } - - /* Check if multiple stimuli should be applied */ - if (count_stimuli < number_of_stimuli) { - /* Update the number of stimuli */ - count_stimuli++; - } else { - /* After last stimulus in event pause the stimulation */ - minimum_found = false; - stimulation_paused = true; - count_to_start = 0; - - /* Reset the stimulus counter for next stimulation event */ - count_stimuli = 1; - } - } - /* Update counter */ - count_to_start++; - } - break; - } - - /* Actual stimulation protocols */ - if(stimulation_started) { - /* Wait to switch the stimulation off */ - if(count_duration==duration) { - stimulation_started = false; - burst_started = true; - count_duration = 0; - count_bursts = 0; - Thalamus->set_input(0.0); - } - - count_duration++; - count_bursts++; - - /* Switch stimulation on and off wrt burst parameters */ - if(burst_enabled) { - if(burst_started) { - if(count_bursts%burst_length==0) { - count_bursts = 0; - burst_started = false; - Thalamus->set_input(0.0); - } - } else { - if(count_bursts%burst_ISI==0) { - count_bursts = 0; - burst_started = true; - Thalamus->set_input(strength); - } - } - } - } - - /* Wait if there is a pause between stimulation events */ - if(stimulation_paused) { - if(count_pause == ISI) { - stimulation_paused = false; - count_pause = 0; - } - count_pause++; - } + /* Check if stimulation should start */ + switch (mode) { + + /* No stimulation */ + default: + break; + + /* Semi-periodic stimulation */ + case 1: + /* Check if stimulation time is reached */ + if(time == time_to_stimuli) { + /* Switch stimulation on */ + stimulation_started = true; + Thalamus->set_input(strength); + + /* Add marker for the first stimuli in the event */ + if(count_stimuli == 1) { + marker_stimulation.push_back(time - onset_correction); + } + + /* Check if multiple stimuli should be applied */ + if (count_stimuli < number_of_stimuli) { + /* Update the timer with respect to time between stimuli */ + time_to_stimuli += time_between_stimuli; + count_stimuli++; + } + /* After last stimulus in event update the timer with respect to (random) ISI*/ + else { + time_to_stimuli += (ISI_range==0)? ISI : Uniform_Distribution(); + + /* Reset the stimulus counter for next stimulation event */ + count_stimuli = 1; + } + } + break; + + /* Phase dependent stimulation */ + case 2: + /* Search for threshold */ + if(!stimulation_started && !minimum_found && !threshold_crossed && time>onset_correction && !stimulation_paused) { + if(Cortex->Vp[0]<=threshold) { + threshold_crossed = true; + } + } + + /* Search for minimum */ + if(threshold_crossed) { + if(Cortex->Vp[0]>Vp_old) { + threshold_crossed = false; + minimum_found = true; + Vp_old = 0; + } else { + Vp_old = Cortex->Vp[0]; + } + } + + /* Wait until the stimulation should start */ + if(minimum_found) { + /* Start stimulation after time_to_stimuli has passed */ + if(count_to_start==time_to_stimuli + (count_stimuli-1) * time_between_stimuli) { + stimulation_started = true; + Thalamus->set_input(strength); + + /* Add marker for the first stimuli in the event */ + if(count_stimuli == 1) { + marker_stimulation.push_back(time - onset_correction); + } + + /* Check if multiple stimuli should be applied */ + if (count_stimuli < number_of_stimuli) { + /* Update the number of stimuli */ + count_stimuli++; + } else { + /* After last stimulus in event pause the stimulation */ + minimum_found = false; + stimulation_paused = true; + count_to_start = 0; + + /* Reset the stimulus counter for next stimulation event */ + count_stimuli = 1; + } + } + /* Update counter */ + count_to_start++; + } + break; + } + + /* Actual stimulation protocols */ + if(stimulation_started) { + /* Wait to switch the stimulation off */ + if(count_duration==duration) { + stimulation_started = false; + burst_started = true; + count_duration = 0; + count_bursts = 0; + Thalamus->set_input(0.0); + } + + count_duration++; + count_bursts++; + + /* Switch stimulation on and off wrt burst parameters */ + if(burst_enabled) { + if(burst_started) { + if(count_bursts%burst_length==0) { + count_bursts = 0; + burst_started = false; + Thalamus->set_input(0.0); + } + } else { + if(count_bursts%burst_ISI==0) { + count_bursts = 0; + burst_started = true; + Thalamus->set_input(strength); + } + } + } + } + + /* Wait if there is a pause between stimulation events */ + if(stimulation_paused) { + if(count_pause == ISI) { + stimulation_paused = false; + count_pause = 0; + } + count_pause++; + } } - -mxArray* Stim::get_marker(void) { - extern const int red; - mxArray* Marker = mxCreateDoubleMatrix(0, 0, mxREAL); - mxSetM(Marker, 1); - mxSetN(Marker, marker_stimulation.size()); - mxSetData(Marker, mxMalloc(sizeof(double)*marker_stimulation.size())); - double* Pr_Marker = mxGetPr(Marker); - for(unsigned i=0; i -#include -#include "ODE.h" +#include +#include "Cortical_Column.h" +#include "ODE.h" +#include "Thalamic_Column.h" -/****************************************************************************************************/ -/* Fixed simulation settings */ -/****************************************************************************************************/ -extern const int T = 30; /* Simulation length s */ -extern const int res = 1E4; /* Number of iteration steps per s */ -extern const int red = 1E2; /* Fraction of iterations that is saved */ -extern const double dt = 1E3/res; /* Duration of a iteration step in ms */ -extern const double h = sqrt(dt); /* Square root of dt for SRK iteration */ -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - +/******************************************************************************/ +/* Fixed simulation settings */ +/******************************************************************************/ +typedef std::chrono::high_resolution_clock::time_point timer; +extern const int T = 30; /* Time until data is stored in s */ +extern const int res = 1E4; /* Number of iteration steps per s */ +extern const double dt = 1E3/res; /* Duration of a time step in ms */ +extern const double h = sqrt(dt); /* Square root of dt for SRK iteration */ -/****************************************************************************************************/ -/* Main simulation routine */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Main simulation routine */ +/******************************************************************************/ int main(void) { - /* Initializing the seeder */ - srand(time(0)); + /* Initializing the populations */ + std::vector param = {6, 1.33, 1E-3}; + std::vector con = {2, 10}; + Cortical_Column Cortex(param.data(), con.data()); - /* Initialize the populations */ - Cortical_Column Cortex; - Thalamic_Column Thalamus; + param = {0.2, 0.06}; + con = {2, 10}; + Thalamic_Column Thalamus (param.data(), con.data()); - /* Connect both modules */ - Cortex.get_Thalamus(Thalamus); - Thalamus.get_Cortex(Cortex); + /* Connect both modules */ + Cortex.get_Thalamus(Thalamus); + Thalamus.get_Cortex(Cortex); - /* Take the time of the simulation */ - time_t start,end; - time (&start); - /* Simulation */ - for (int t=0; t< T*res; ++t) { - ODE (Cortex, Thalamus); - } + /* Take the time of the simulation */ + time_t start,end; + time (&start); + /* Simulation */ + for (unsigned t=0; t< T*res; ++t) { + ODE(Cortex, Thalamus); + } - time (&end); - /* Time consumed by the simulation */ - double dif = difftime(end,start); - std::cout << "simulation done!\n"; - std::cout << "took " << dif << " seconds" << "\n"; - std::cout << "end\n"; + time (&end); + /* Time consumed by the simulation */ + double dif = difftime(end,start); + std::cout << "simulation done!\n"; + std::cout << "took " << dif << " seconds" << "\n"; + std::cout << "end\n"; } /****************************************************************************************************/ /* end */ diff --git a/TC_mex.cpp b/TC_mex.cpp index 1e18e0f..b507bc6 100644 --- a/TC_mex.cpp +++ b/TC_mex.cpp @@ -25,108 +25,124 @@ * to auditory stimulation. * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, * JC Claussen. - * PLoS Computational Biology In Review (in review). + * PLoS Computational Biology http://dx.doi.org/10.1371/journal.pcbi.1005022 */ -/****************************************************************************************************/ -/* Implementation of the simulation as MATLAB routine (mex compiler) */ -/* mex command is given by: */ -/* mex CXXFLAGS="\$CXXFLAGS -std=c++11" TC_mex.cpp Cortical_Column.cpp Thalamic_Column.cpp */ -/****************************************************************************************************/ -#include +/******************************************************************************/ +/* Implementation of the simulation as MATLAB routine (mex compiler) */ +/* mex command is given by: */ +/* mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" TC_mex.cpp Cortical_Column.cpp */ +/* Thalamic_Column.cpp */ +/******************************************************************************/ #include "mex.h" #include "matrix.h" + +#include +#include + +#include "Cortical_Column.h" #include "Data_Storage.h" #include "ODE.h" #include "Stimulation.h" +#include "Thalamic_Column.h" mxArray* GetMexArray(int N, int M); +mxArray* get_marker(Stim &stim); -/****************************************************************************************************/ -/* Fixed simulation settings */ -/****************************************************************************************************/ -extern const int onset = 20; /* Time until data is stored in s */ -extern const int res = 1E4; /* Number of iteration steps per s */ -extern const int red = 1E2; /* Fraction of iterations that is saved */ -extern const double dt = 1E3/res; /* Duration of a iteration step in ms */ -extern const double h = sqrt(dt); /* Square root of dt for SRK iteration */ -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Simulation routine */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Fixed simulation settings */ +/******************************************************************************/ +extern const int onset = 20; /* Time until data is stored in s */ +extern const int res = 1E4; /* Number of iteration steps per s */ +extern const int red = 1E2; /* Number of iterations steps not saved */ +extern const double dt = 1E3/res; /* Duration of a time step in ms */ +extern const double h = sqrt(dt); /* Square root of dt for SRK iteration */ + +/******************************************************************************/ +/* Simulation routine */ +/* lhs defines outputs */ +/* rhs defines inputs */ +/******************************************************************************/ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - /* Initialize the seeder */ - srand(time(NULL)); - - /* Fetch inputs */ - const int T = (int) (mxGetScalar(prhs[0])); /* Duration of simulation in s */ - const int Time = (T+onset)*res; /* Total number of iteration steps */ - double* Param_Cortex = mxGetPr (prhs[1]); /* Parameters of cortical module */ - double* Param_Thalamus = mxGetPr (prhs[2]); /* Parameters of thalamic module */ - double* Connections = mxGetPr (prhs[3]); /* Connectivity values C <-> T */ - double* var_stim = mxGetPr (prhs[4]); /* Parameters of stimulation protocol */ - - /* Initialize the populations */ - Cortical_Column Cortex = Cortical_Column(Param_Cortex, Connections); - Thalamic_Column Thalamus = Thalamic_Column(Param_Thalamus, Connections); - - /* Link both modules */ - Cortex.get_Thalamus(Thalamus); - Thalamus.get_Cortex(Cortex); - - /* Initialize the stimulation protocol */ - Stim Stimulation(Cortex, Thalamus, var_stim); - - /* Create data containers */ - vector Data; - Data.push_back(GetMexArray(1, T*res/red)); // Vt - Data.push_back(GetMexArray(1, T*res/red)); // Vr - Data.push_back(GetMexArray(1, T*res/red)); // Ca - Data.push_back(GetMexArray(1, T*res/red)); // act_h - - /* Pointer to the actual data block */ - vector pData(Data.size(), NULL); - for(unsigned i=0; i=onset*res && t%red==0){ - get_data(count, Cortex, Thalamus, pData); - ++count; - } - } - - /* Output of the simulation */ - /* Return the data containers */ - for(unsigned i=0; i T */ + double* var_stim = mxGetPr (prhs[4]); /* Parameters of stimulation protocol */ + + /* Initialize the populations */ + Cortical_Column Cortex = Cortical_Column(Param_Cortex, Connections); + Thalamic_Column Thalamus = Thalamic_Column(Param_Thalamus, Connections); + + /* Link both modules */ + Cortex.get_Thalamus(Thalamus); + Thalamus.get_Cortex(Cortex); + + /* Initialize the stimulation protocol */ + Stim Stimulation(Cortex, Thalamus, var_stim); + /* Create data containers */ + std::vector dataArray; + dataArray.reserve(4); + dataArray.push_back(GetMexArray(1, T*res/red)); // Vt + dataArray.push_back(GetMexArray(1, T*res/red)); // Vr + dataArray.push_back(GetMexArray(1, T*res/red)); // Ca + dataArray.push_back(GetMexArray(1, T*res/red)); // act_h -/****************************************************************************************************/ -/* Create MATLAB data container */ -/****************************************************************************************************/ + /* Pointer to the data blocks */ + std::vector dataPointer; + dataPointer.reserve(dataArray.size()); + for (auto &dataptr : dataArray) { + dataPointer.push_back(mxGetPr(dataptr)); + } + + /* Simulation */ + int count = 0; + for (unsigned t=0; t < Time; ++t) { + ODE (Cortex, Thalamus); + Stimulation.check_stim(t); + if(t >= onset*res && t%red == 0){ + get_data(count, Cortex, Thalamus, dataPointer); + ++count; + } + } + + /* Return the data containers */ + nlhs = dataArray.size()+1; + for (auto &dataptr : dataArray) { + plhs[std::distance(&dataptr, dataArray.data())] = dataptr; + } + plhs[dataArray.size()] = get_marker(Stimulation); + + return; +} + +/******************************************************************************/ +/* Create MATLAB data containers */ +/******************************************************************************/ mxArray* GetMexArray(int N, int M) { - mxArray* Array = mxCreateDoubleMatrix(0, 0, mxREAL); - mxSetM(Array, N); - mxSetN(Array, M); -#pragma omp critical - {mxSetData(Array, mxMalloc(sizeof(double)*M*N));} - return Array; + mxArray* Array = mxCreateDoubleMatrix(0, 0, mxREAL); + mxSetM(Array, N); + mxSetN(Array, M); + {mxSetData(Array, mxMalloc(sizeof(double)*M*N));} + return Array; +} + +mxArray* get_marker(Stim &stim) { + extern const int red; + mxArray* marker = mxCreateDoubleMatrix(0, 0, mxREAL); + mxSetM(marker, 1); + mxSetN(marker, stim.marker_stimulation.size()); + mxSetData(marker, mxMalloc(sizeof(double)*stim.marker_stimulation.size())); + double* Pr_Marker = mxGetPr(marker); + unsigned counter = 0; + /* Division by res transforms marker time from dt to sampling rate */ + for(auto & elem : stim.marker_stimulation) { + Pr_Marker[counter++] = elem/red; + } + return marker; } -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index d021df3..c697ffc 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -25,282 +25,220 @@ * to auditory stimulation. * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, * JC Claussen. - * PLoS Computational Biology (in review). + * PLoS Computational Biology http://dx.doi.org/10.1371/journal.pcbi.1005022 */ -/****************************************************************************************************/ -/* Functions of the thalamic module */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Functions of the thalamic module */ +/******************************************************************************/ #include "Thalamic_Column.h" -/****************************************************************************************************/ -/* Initialization of RNG */ -/****************************************************************************************************/ +/******************************************************************************/ +/* Initialization of RNG */ +/******************************************************************************/ void Thalamic_Column::set_RNG(void) { - extern const double dt; - /* Number of independent random variables */ - int N = 1; - - /* Create RNG for each stream */ - for (int i=0; iy[N] - s_et[N]) - 2 * gamma_e * x_et[N]) + noise_xRK(N,0); - x_er [N+1] = x_er [0] + A[N]*dt*(pow(gamma_e, 2) * (N_rt * get_Qt(N) + N_rp * Cortex->y[N] - s_er[N]) - 2 * gamma_e * x_er[N]); - x_gt [N+1] = x_gt [0] + A[N]*dt*(pow(gamma_g, 2) * (N_tr * get_Qr(N) - s_gt[N]) - 2 * gamma_g * x_gt[N]); - x_gr [N+1] = x_gr [0] + A[N]*dt*(pow(gamma_g, 2) * (N_rr * get_Qr(N) - s_gr[N]) - 2 * gamma_g * x_gr[N]); - x [N+1] = x [0] + A[N]*dt*(pow(nu, 2) * ( get_Qt(N) - y [N]) - 2 * nu * x [N]); + extern const double dt; + Vt [N+1] = Vt [0] + A[N]*dt*(-(I_L_t(N) + I_et(N) + I_gt(N))/tau_t - C_m * (I_LK_t(N) + I_T_t(N) + I_h(N))); + Vr [N+1] = Vr [0] + A[N]*dt*(-(I_L_r(N) + I_er(N) + I_gr(N))/tau_r - C_m * (I_LK_r(N) + I_T_r(N))); + Ca [N+1] = Ca [0] + A[N]*dt*(alpha_Ca * I_T_t(N) - (Ca[N] - Ca_0)/tau_Ca); + h_T_t [N+1] = h_T_t[0] + A[N]*dt*(h_inf_T_t(N) - h_T_t[N])/tau_h_T_t(N); + h_T_r [N+1] = h_T_r[0] + A[N]*dt*(h_inf_T_r(N) - h_T_r[N])/tau_h_T_r(N); + m_h [N+1] = m_h [0] + A[N]*dt*((m_inf_h(N) * (1 - m_h2[N]) - m_h[N])/tau_m_h(N) - k3 * P_h(N) * m_h[N] + k4 * m_h2[N]); + m_h2 [N+1] = m_h2 [0] + A[N]*dt*(k3 * P_h(N) * m_h[N] - k4 * m_h2[N]); + s_et [N+1] = s_et [0] + A[N]*dt*(x_et[N]); + s_er [N+1] = s_er [0] + A[N]*dt*(x_er[N]); + s_gt [N+1] = s_gt [0] + A[N]*dt*(x_gt[N]); + s_gr [N+1] = s_gr [0] + A[N]*dt*(x_gr[N]); + y [N+1] = y [0] + A[N]*dt*(x [N]); + x_et [N+1] = x_et [0] + A[N]*dt*(gamma_e*gamma_e * ( + N_tp * Cortex->y[N] - s_et[N]) - 2 * gamma_e * x_et[N]) + noise_xRK(N,0); + x_er [N+1] = x_er [0] + A[N]*dt*(gamma_e*gamma_e * (N_rt * get_Qt(N) + N_rp * Cortex->y[N] - s_er[N]) - 2 * gamma_e * x_er[N]); + x_gt [N+1] = x_gt [0] + A[N]*dt*(gamma_g*gamma_g * (N_tr * get_Qr(N) - s_gt[N]) - 2 * gamma_g * x_gt[N]); + x_gr [N+1] = x_gr [0] + A[N]*dt*(gamma_g*gamma_g * (N_rr * get_Qr(N) - s_gr[N]) - 2 * gamma_g * x_gr[N]); + x [N+1] = x [0] + A[N]*dt*(nu * nu * ( get_Qt(N) - y [N]) - 2 * nu * x [N]); } -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - -/****************************************************************************************************/ -/* Function that adds all SRK terms */ -/****************************************************************************************************/ void Thalamic_Column::add_RK(void) { - Vt [0] =(-3*Vt [0] + 2*Vt [1] + 4*Vt [2] + 2* Vt [3] + Vt [4])/6; - Vr [0] =(-3*Vr [0] + 2*Vr [1] + 4*Vr [2] + 2* Vr [3] + Vr [4])/6; - Ca [0] =(-3*Ca [0] + 2*Ca [1] + 4*Ca [2] + 2* Ca [3] + Ca [4])/6; - s_et [0] =(-3*s_et [0] + 2*s_et [1] + 4*s_et [2] + 2* s_et [3] + s_et [4])/6; - s_er [0] =(-3*s_er [0] + 2*s_er [1] + 4*s_er [2] + 2* s_er [3] + s_er [4])/6; - s_gt [0] =(-3*s_gt [0] + 2*s_gt [1] + 4*s_gt [2] + 2* s_gt [3] + s_gt [4])/6; - s_gr [0] =(-3*s_gr [0] + 2*s_gr [1] + 4*s_gr [2] + 2* s_gr [3] + s_gr [4])/6; - y [0] =(-3*y [0] + 2*y [1] + 4*y [2] + 2* y [3] + y [4])/6; - x_et [0] =(-3*x_et [0] + 2*x_et [1] + 4*x_et [2] + 2* x_et [3] + x_et [4])/6 + noise_aRK(0); - x_er [0] =(-3*x_er [0] + 2*x_er [1] + 4*x_er [2] + 2* x_er [3] + x_er [4])/6; - x_gt [0] =(-3*x_gt [0] + 2*x_gt [1] + 4*x_gt [2] + 2* x_gt [3] + x_gt [4])/6; - x_gr [0] =(-3*x_gr [0] + 2*x_gr [1] + 4*x_gr [2] + 2* x_gr [3] + x_gr [4])/6; - x [0] =(-3*x [0] + 2*x [1] + 4*x [2] + 2* x [3] + x [4])/6; - h_T_t [0] =(-3*h_T_t[0] + 2*h_T_t[1] + 4*h_T_t[2] + 2* h_T_t [3] + h_T_t [4])/6; - h_T_r [0] =(-3*h_T_r[0] + 2*h_T_r[1] + 4*h_T_r[2] + 2* h_T_r [3] + h_T_r [4])/6; - m_h [0] =(-3*m_h [0] + 2*m_h [1] + 4*m_h [2] + 2* m_h [3] + m_h [4])/6; - m_h2 [0] =(-3*m_h2 [0] + 2*m_h2 [1] + 4*m_h2 [2] + 2* m_h2 [3] + m_h2 [4])/6; - /* Generate noise for the next iteration */ - for (unsigned i=0; i #include -#include "Random_Stream.h" + #include "Cortical_Column.h" -using std::vector; +#include "Random_Stream.h" class Cortical_Column; -/****************************************************************************************************/ -/* Macro for vector initialization */ -/****************************************************************************************************/ -#ifndef _INIT -#define _INIT(x) {x, 0.0, 0.0, 0.0, 0.0} -#endif -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Implementation of the thalamic module */ -/****************************************************************************************************/ class Thalamic_Column { public: - /* Constructors for compile check */ - Thalamic_Column(void) - {set_RNG();} - - /* Constructor for simulation */ - Thalamic_Column(double* Param, double* Con) - : g_LK (Param[0]), g_h (Param[1]), - N_tp (Con[0]), N_rp (Con[1]) - {set_RNG();} - - /* Get the pointer to the cortical module */ - void get_Cortex (Cortical_Column& C) {Cortex = &C;} + /* Constructor for simulation */ + Thalamic_Column(double* Param, double* Con) + : g_LK (Param[0]), g_h (Param[1]), + N_tp (Con[0]), N_rp (Con[1]) + {set_RNG();} - /* ODE functions */ - void set_RK (int); - void add_RK (void); + /* Get the pointer to the cortical module */ + void get_Cortex (Cortical_Column& C) {Cortex = &C;} - /* Data storage access */ - friend void get_data (int, Cortical_Column&, Thalamic_Column&, vector); - friend class Cortical_Column; + /* ODE functions */ + void set_RK (int); + void add_RK (void); - /* Set strength of external input */ - void set_input (double I) {input = I;} + /* Set strength of external input */ + void set_input (double I) {input = I;} private: - /* Declaration of private functions */ - /* Initialize the RNGs */ - void set_RNG (void); - - /* Firing rates */ - double get_Qt (int) const; - double get_Qr (int) const; - - /* Synaptic currents */ - double I_et (int) const; - double I_gt (int) const; - double I_er (int) const; - double I_gr (int) const; - - /* Activation functions */ - double m_inf_T_t (int) const; - double m_inf_T_r (int) const; - double m_inf_h (int) const; - double tau_m_h (int) const; - double P_h (int) const; - double act_h (void)const; - - double m_inf_hs (int) const; - double tau_m_hs (int) const; - double tau_m_hf (int) const; - - /* Deactivation functions */ - double h_inf_T_t (int) const; - double h_inf_T_r (int) const; - double tau_h_T_t (int) const; - double tau_h_T_r (int) const; - - /* Currents */ - double I_L_t (int) const; - double I_L_r (int) const; - double I_LK_t (int) const; - double I_LK_r (int) const; - double I_T_t (int) const; - double I_T_r (int) const; - double I_h (int) const; - - /* Noise functions */ - double noise_xRK (int,int) const; - double noise_aRK (int) const; - - /* Declaration and Initialization of parameters */ - /* Membrane time in ms */ - const double tau_t = 20; - const double tau_r = 20; - - /* Maximum firing rate in ms^-1 */ - const double Qt_max = 400.E-3; - const double Qr_max = 400.E-3; - - /* Sigmoid threshold in mV */ - const double theta_t = -58.5; - const double theta_r = -58.5; - - /* Sigmoid gain in mV */ - const double sigma_t = 6.; - const double sigma_r = 6.; - - /* Scaling parameter for sigmoidal mapping (dimensionless) */ - const double C1 = (M_PI/sqrt(3)); - - /* PSP rise time in ms^-1 */ - const double gamma_e = 70E-3; - const double gamma_g = 100E-3; - - /* Axonal flux time constant in ms^-1*/ - const double nu = 120E-3; - - /* Membrane capacitance in muF/cm^2 */ - const double C_m = 1.; - - /* Leak weight in aU */ - const double g_L = 1.; - - /* Synaptic weights in ms */ - const double g_AMPA = 1.; - const double g_GABA = 1.; - - /* Conductivities */ - /* Potassium leak current in mS/m^2 */ - const double g_LK = 0.02; - - /* T current in mS/m^2 */ - const double g_T_t = 3; - const double g_T_r = 2.3; - - /* h current in mS/m^2 */ - const double g_h = 0.06; - - /* Reversal potentials in mV */ - /* Synaptic */ - const double E_AMPA = 0; - const double E_GABA = -70; - - /* Leak */ - const double E_L_t = -70; - const double E_L_r = -70; - - /* Potassium */ - const double E_K = -100; - - /* I_T current */ - const double E_Ca = 120; - - /* I_h current */ - const double E_h = -40; - - /* Calcium parameters */ - const double alpha_Ca = -51.8E-6; /* influx per spike in nmol */ - const double tau_Ca = 10; /* calcium time constant in ms */ - const double Ca_0 = 2.4E-4; /* resting concentration */ - - /* I_h activation parameters */ - const double k1 = 2.5E7; - const double k2 = 4E-4; - const double k3 = 1E-1; - const double k4 = 1E-3; - const double n_P = 4; - const double g_inc = 2; - - /* Noise parameters in ms^-1 */ - const double mphi = 0E-3; - const double dphi = 20E-3; - double input = 0.0; - - /* Connectivities (dimensionless) */ - const double N_rt = 3.; - const double N_tr = 5.; - const double N_rr = 25.; - - /* Connectivities from cortex (dimensionless) */ - const double N_tp = 2.6; - const double N_rp = 2.6; - - /* Pointer to cortical column */ - Cortical_Column* Cortex; - - /* Parameters for SRK4 iteration */ - const vector A = {0.5, 0.5, 1.0, 1.0}; - const vector B = {0.75, 0.75, 0.0, 0.0}; - - /* Random number generators */ - vector MTRands; - - /* Container for noise */ - vector Rand_vars; - - /* Population variables */ - vector Vt = _INIT(E_L_t), /* TC membrane voltage */ - Vr = _INIT(E_L_r), /* RE membrane voltage */ - Ca = _INIT(Ca_0), /* Calcium concentration of TC population */ - s_et = _INIT(0.0), /* PostSP from TC population to TC population */ - s_er = _INIT(0.0), /* PostSP from TC population to RE population */ - s_gt = _INIT(0.0), /* PostSP from RE population to TC population */ - s_gr = _INIT(0.0), /* PostSP from RE population to RE population */ - y = _INIT(0.0), /* axonal flux */ - x_et = _INIT(0.0), /* derivative of s_et */ - x_er = _INIT(0.0), /* derivative of s_er */ - x_gt = _INIT(0.0), /* derivative of s_gt */ - x_gr = _INIT(0.0), /* derivative of s_gr */ - x = _INIT(0.0), /* derivative of y */ - h_T_t = _INIT(0.0), /* inactivation of T channel */ - h_T_r = _INIT(0.0), /* inactivation of T channel */ - m_h = _INIT(0.0), /* activation of h channel */ - m_h2 = _INIT(0.0); /* activation of h channel bound with protein */ + /* Declaration of private functions */ + /* Initialize the RNGs */ + void set_RNG (void); + + /* Firing rates */ + double get_Qt (int) const; + double get_Qr (int) const; + + /* Synaptic currents */ + double I_et (int) const; + double I_gt (int) const; + double I_er (int) const; + double I_gr (int) const; + + /* Activation functions */ + double m_inf_T_t (int) const; + double m_inf_T_r (int) const; + double m_inf_h (int) const; + double tau_m_h (int) const; + double P_h (int) const; + double act_h (void)const; + + double m_inf_hs (int) const; + double tau_m_hs (int) const; + double tau_m_hf (int) const; + + /* Deactivation functions */ + double h_inf_T_t (int) const; + double h_inf_T_r (int) const; + double tau_h_T_t (int) const; + double tau_h_T_r (int) const; + + /* Currents */ + double I_L_t (int) const; + double I_L_r (int) const; + double I_LK_t (int) const; + double I_LK_r (int) const; + double I_T_t (int) const; + double I_T_r (int) const; + double I_h (int) const; + + /* Noise functions */ + double noise_xRK (int,int) const; + double noise_aRK (int) const; + + /* Helper functions */ + inline std::vector init (double value) + {return {value, 0.0, 0.0, 0.0, 0.0};} + + inline void add_RK (std::vector& var) + {var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6;} + + inline void add_RK_noise (std::vector& var, unsigned noise) + {var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6 + noise_aRK(noise);} + + /* Declaration and Initialization of parameters */ + /* Membrane time in ms */ + const double tau_t = 20; + const double tau_r = 20; + + /* Maximum firing rate in ms^-1 */ + const double Qt_max = 400.E-3; + const double Qr_max = 400.E-3; + + /* Sigmoid threshold in mV */ + const double theta_t = -58.5; + const double theta_r = -58.5; + + /* Sigmoid gain in mV */ + const double sigma_t = 6.; + const double sigma_r = 6.; + + /* Scaling parameter for sigmoidal mapping (dimensionless) */ + const double C1 = (M_PI/sqrt(3)); + + /* PSP rise time in ms^-1 */ + const double gamma_e = 70E-3; + const double gamma_g = 100E-3; + + /* Axonal flux time constant in ms^-1*/ + const double nu = 120E-3; + + /* Membrane capacitance in muF/cm^2 */ + const double C_m = 1.; + + /* Leak weight in aU */ + const double g_L = 1.; + + /* Synaptic weights in ms */ + const double g_AMPA = 1.; + const double g_GABA = 1.; + + /* Conductivities */ + /* Potassium leak current in mS/m^2 */ + const double g_LK = 0.02; + + /* T current in mS/m^2 */ + const double g_T_t = 3; + const double g_T_r = 2.3; + + /* h current in mS/m^2 */ + const double g_h = 0.06; + + /* Reversal potentials in mV */ + /* Synaptic */ + const double E_AMPA = 0; + const double E_GABA = -70; + + /* Leak */ + const double E_L_t = -70; + const double E_L_r = -70; + + /* Potassium */ + const double E_K = -100; + + /* I_T current */ + const double E_Ca = 120; + + /* I_h current */ + const double E_h = -40; + + /* Calcium parameters */ + const double alpha_Ca = -51.8E-6; /* influx per spike in nmol */ + const double tau_Ca = 10; /* calcium time constant in ms */ + const double Ca_0 = 2.4E-4; /* resting concentration */ + + /* I_h activation parameters */ + const double k1 = 2.5E7; + const double k2 = 4E-4; + const double k3 = 1E-1; + const double k4 = 1E-3; + const double n_P = 4; + const double g_inc = 2; + + /* Noise parameters in ms^-1 */ + const double mphi = 0E-3; + const double dphi = 20E-3; + double input = 0.0; + + /* Connectivities (dimensionless) */ + const double N_rt = 3.; + const double N_tr = 5.; + const double N_rr = 25.; + + /* Connectivities from cortex (dimensionless) */ + const double N_tp = 2.6; + const double N_rp = 2.6; + + /* Pointer to cortical column */ + Cortical_Column* Cortex; + + /* Parameters for SRK4 iteration */ + const std::vector A = {0.5, 0.5, 1.0, 1.0}; + const std::vector B = {0.75, 0.75, 0.0, 0.0}; + + /* Random number generators */ + std::vector MTRands; + + /* Container for noise */ + std::vector Rand_vars; + + /* Population variables */ + std::vector Vt = init(E_L_t), /* TC membrane voltage */ + Vr = init(E_L_r), /* RE membrane voltage */ + Ca = init(Ca_0), /* Calcium concentration of TC population */ + s_et = init(0.0), /* PostSP from TC population to TC population */ + s_er = init(0.0), /* PostSP from TC population to RE population */ + s_gt = init(0.0), /* PostSP from RE population to TC population */ + s_gr = init(0.0), /* PostSP from RE population to RE population */ + y = init(0.0), /* axonal flux */ + x_et = init(0.0), /* derivative of s_et */ + x_er = init(0.0), /* derivative of s_er */ + x_gt = init(0.0), /* derivative of s_gt */ + x_gr = init(0.0), /* derivative of s_gr */ + x = init(0.0), /* derivative of y */ + h_T_t = init(0.0), /* inactivation of T channel */ + h_T_r = init(0.0), /* inactivation of T channel */ + m_h = init(0.0), /* activation of h channel */ + m_h2 = init(0.0); /* activation of h channel bound with protein */ + + /* Data storage access */ + friend void get_data (int, Cortical_Column&, Thalamic_Column&, std::vector); + friend class Cortical_Column; }; /****************************************************************************************************/ /* end */