From f5a83cfa4e07cc1d08536030ee2b84021d196dac Mon Sep 17 00:00:00 2001 From: Schellenberger Date: Thu, 3 Apr 2014 11:40:32 +0200 Subject: [PATCH] Moved paramters into the class definition. That makes it easier to modify them at initialization. Signed-off-by: Schellenberger --- .cproject | 15 +++--- Plots.m | 13 +++-- Stimulation.h | 123 ++++++++++++++++++++++++++---------------- Thalamic_Column.cpp | 38 ++++++------- Thalamic_Column.h | 150 +++++++++++++++++++++++++++++++++++++--------------- Thalamus.cpp | 27 +++++----- parameters.h | 122 ------------------------------------------ 7 files changed, 230 insertions(+), 258 deletions(-) delete mode 100644 parameters.h diff --git a/.cproject b/.cproject index 9a05466..1fee099 100644 --- a/.cproject +++ b/.cproject @@ -1,18 +1,16 @@ - - - + - + @@ -23,7 +21,7 @@ - diff --git a/Plots.m b/Plots.m index 49df712..8c2457c 100644 --- a/Plots.m +++ b/Plots.m @@ -1,21 +1,20 @@ % mex command is given by: % mex CXXFLAGS="\$CXXFLAGS -std=gnu++0x -fpermissive" Thalamus.cpp Thalamic_Column.cpp -function Plots(T, onset) +function Plots(T) if nargin == 0 - Con = [10; % N_tr - 20; % N_rt - 40]; % N_rr + Con = [2; % N_tr + 5.5; % N_rt + 5]; % N_rr var_stim = [ 0; % strength of the stimulus in Hz (spikes per second) 0; % time between stimuli in s 0; % time until first stimuli in s 0]; % duration of the stimulus in ms - T = 30; % duration of the simulation + T = 60; % duration of the simulation end - -[Vt] = Thalamus(T, Con, var_stim); +[Vt, Vr] = Thalamus(T, Con, var_stim); L = max(size(Vt)); timeaxis = linspace(0,T,L); diff --git a/Stimulation.h b/Stimulation.h index 371f617..fccdf15 100644 --- a/Stimulation.h +++ b/Stimulation.h @@ -9,86 +9,115 @@ /****************************************************************************************************/ class Stim { public: - // empty constructor for compiling + /* empty constructor for compiling */ Stim(void); - Stim(Thalamic_Column& C, double* var) - {Thalamus = &T; - setup(var);} + Stim(Thalamic_Column& T, double* var) + { Thalamus = &T; + setup(var);} - // scaling from SI to simulation variables s -> ms + /* setup with respect to stimulation mode */ void setup (double* var_stim) { extern const int onset; extern const int res; extern const int dt; - // scale the stimulation strength with respect to ms^-1 - strength = var_stim[0] / 1000; + /* mode of stimulation */ + mode = (int) var_stim[0]; - // scale the stimulation variables with respect to simulation resolution - ISI = (int) var_stim[1] * res; + /* scale the stimulation strength from s^-1 to ms^-1 */ + strength = var_stim[1] / 1000; - // stimulation starts after the onset - start = (int)(var_stim[2] + onset) *res; + /* scale duration from ms to dt */ + duration = (int) var_stim[2] * res / 1000; + + /* scale the ISI from s to ms^-1 */ + ISI = (int) var_stim[3] * res; + + /* scale time to stimulus from ms to dt */ + time_to_stim= (int) var_stim[4] * res / 1000; + + if(mode==1) { + time_to_stim = (onset+1) * res; + } + + correction = onset * res; - // rescale duration with respect to dt - duration = (int) var_stim[3]/dt; } void check_stim (int time) { - // check whether a stimulation should start - // the duration of the stimulation is ignored - if(time==(start + count_stim*ISI)){ - // turn the stimulation on - mode = 1; - Thalamus->set_input(strength); - } + /* check if stimulation should start */ + switch (mode) { - // check whether a stimulation should end - if(mode ==1 && count_dur ==duration) { - // turn off the stimulation - mode = 0; - Thalamus->set_input(0.0); + /* no stimulation */ + default: + break; - // add counter for stimulation occurrence - count_stim++ ; + /* periodic stimulation */ + case 1: + /* check if time is reached */ + if(time == time_to_stim) { + /* switch stimulation on */ + stimulation_started = true; + Thalamus->set_input(strength); + + /* update the timer */ + time_to_stim += ISI; + } + break; - // reset the stimulation counter - count_dur = 0; } - // if stimulation is on track count its duration - if(mode==1){ - count_dur++; + /* wait to switch the stimulation off */ + if(stimulation_started) { + count_duration++; + + /* switch stimulation off */ + if(count_duration==duration) { + stimulation_started = false; + count_duration = 0; + Thalamus->set_input(0.0); + } } } private: - // stimulation strength - double strength = 0.0; + /* Stimulation parameters */ + /* stimulation strength */ + double strength = 0.0; + + /* duration of the stimulation */ + int duration = 500; - // inter stimulus interval - int ISI = 0; + /* inter stimulus intervall in case of periodic stimulation */ + int ISI = 5E4; - // onset until stimulation starts - int start = 0; + /* time until stimulus after minimum was found */ + int time_to_stim = 5500; - // duration of the stimulation - int duration = 0; + /* mode of stimulation */ + /* 0 == none */ + /* 1 == periodic */ + int mode = 0; - // counter for stimulation events - int count_stim = 0; + /* Internal variables */ + /* Simulation on for TRUE and off for FALSE */ + bool stimulation_started= false; - // counter for stimulation length - int count_dur = 0; + /* onset in timesteps to correct the given time of the markers */ + int correction = 10000; - // Simulation on for TRUE and off for FALSE - bool mode = 0; + /* counter for stimulation duration */ + int count_duration = 0; - // Pointer to thalamic module + /* counter after minimum */ + int count_to_start = 0; + + /* pointer to thalamic column */ Thalamic_Column* Thalamus; + }; /****************************************************************************************************/ /* end */ diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index 37170e8..2fc8c32 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -13,7 +13,7 @@ void Thalamic_Column::set_RNG(void) { // get the RNG for (int i=0; i #include #include "macros.h" -#include "parameters.h" using std::vector; /****************************************************************************************************/ @@ -29,30 +28,19 @@ class Thalamic_Column { public: // Constructors Thalamic_Column(void) - : Vt (_INIT(E_L_t)), Vr (_INIT(E_L_r)), Ca (_INIT(Ca_0)), - Phi_tt(_INIT(0.0)), Phi_tr (_INIT(0.0)), Phi_rt (_INIT(0.0)), Phi_rr (_INIT(0.0)), - x_tt (_INIT(0.0)), x_tr (_INIT(0.0)), x_rt (_INIT(0.0)), x_rr (_INIT(0.0)), - h_T_t (_INIT(0.0)), h_T_r (_INIT(0.0)), m_T_t (_INIT(0.0)), m_T_r (_INIT(0.0)), - m_h (_INIT(0.0)), m_h2 (_INIT(0.0)), P_h (_INIT(0.0)), - N_tr (210), N_rt (410), N_rr (800), input (0.0) {set_RNG();} // Constructors Thalamic_Column(double* Par) - : Vt (_INIT(0)), Vr (_INIT(0)), Ca (_INIT(Ca_0)), - Phi_tt(_INIT(0.0)), Phi_tr (_INIT(0.0)), Phi_rt (_INIT(0.0)), Phi_rr (_INIT(0.0)), - x_tt (_INIT(0.0)), x_tr (_INIT(0.0)), x_rt (_INIT(0.0)), x_rr (_INIT(0.0)), - h_T_t (_INIT(0.0)), h_T_r (_INIT(0.0)), m_T_t (_INIT(0.0)), m_T_r (_INIT(0.0)), - m_h (_INIT(0.0)), m_h2 (_INIT(0.0)), P_h (_INIT(0.0)), - N_tr (Par[0]), N_rt (Par[1]), N_rr (Par[2]), input (0.0) + : N_tr (Par[0]), N_rt (Par[1]), N_rr (Par[2]) {set_RNG();} - // Initialize the RNGs - void set_RNG (void); - // change the strength of input void set_input (double I) {input = I;} + // Initialize the RNGs + void set_RNG (void); + // Firing rates double get_Qt (int) const; double get_Qr (int) const; @@ -98,38 +86,116 @@ class Thalamic_Column { private: // Population variables - vector Vt, // TC membrane voltage - Vr, // RE membrane voltage - Ca, // Calcium concentration of TC population - Phi_tt, // PostSP from TC population to TC population - Phi_tr, // PostSP from TC population to RE population - Phi_rt, // PostSP from RE population to TC population - Phi_rr, // PostSP from RE population to RE population - x_tt, // derivative of Phi_tt - x_tr, // derivative of Phi_tr - x_rt, // derivative of Phi_rt - x_rr, // derivative of Phi_rr - h_T_t, // inactivation of T channel - h_T_r, // inactivation of T channel - m_T_t, // activation of T channel - m_T_r, // activation of T channel - m_h, // activation of h channel - m_h2, // activation of h channel bound with protein - P_h; // fraction of protein bound with calcium - - // Connectivities - double N_tr, // TC to RE loop - N_rt, // RE to TC loop - N_rr; // RE self loop - - // Noise parameters - double input; + 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 + Phi_tt = _INIT(0.0), // PostSP from TC population to TC population + Phi_tr = _INIT(0.0), // PostSP from TC population to RE population + Phi_rt = _INIT(0.0), // PostSP from RE population to TC population + Phi_rr = _INIT(0.0), // PostSP from RE population to RE population + x_tt = _INIT(0.0), // derivative of Phi_tt + x_tr = _INIT(0.0), // derivative of Phi_tr + x_rt = _INIT(0.0), // derivative of Phi_rt + x_rr = _INIT(0.0), // derivative of Phi_rr + h_T_t = _INIT(0.0), // inactivation of T channel + h_T_r = _INIT(0.0), // inactivation of T channel + m_T_t = _INIT(0.0), // activation of T channel + m_T_r = _INIT(0.0), // activation of T channel + m_h = _INIT(0.0), // activation of h channel + m_h2 = _INIT(0.0), // activation of h channel bound with protein + P_h = _INIT(0.0); // fraction of protein bound with calcium // Random number generators vector MTRands; // Container for noise vector Rand_vars; + + // Declaration and Initialization of parameters + // Membrane time in ms + const double tau_t = 1; + const double tau_r = 1; + + // 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 = -45; + const double theta_r = -45; + + // Sigmoid gain in mV + const double sigma_t = 9; + const double sigma_r = 9; + + // Scaling parameter for sigmoidal mapping (dimensionless) + const double C1 = (3.14159265/sqrt(3)); + + // PSP rise time in ms^-1 + const double gamma_e = 70E-3; + const double gamma_i = 58.6E-3; + + /* axonal flux time constant */ + const double nu = 60E-3; + + // Conductivities in mS/cm^-2 + // Leak current + const double g_L_t = 0.02; + const double g_L_r = 0.05; + + // Potassium leak current + const double g_LK_t = 0.02; + const double g_LK_r = 0.02; + + // T current + const double g_T_t = 3; + const double g_T_r = 2; + + /* h current */ + const double g_h = 0.07; + + // Connectivities (dimensionless) + const double N_tr = 2; + const double N_rt = 5.5; + const double N_rr = 5; + + // 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 = -55; + + // 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 = -50E-6; // influx per spike in nmol + const double tau_Ca = 5; // calcium time constant in ms + const double Ca_0 = 2E-4; // resting concentration + + /* I_h activation parameters */ + const double k1 = 2.5E7; + const double k2 = 5E-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 = 2E-3;; + double input = 0.0; + + friend class Stim; }; /****************************************************************************************************/ /* end */ diff --git a/Thalamus.cpp b/Thalamus.cpp index 09afdd8..795a86c 100644 --- a/Thalamus.cpp +++ b/Thalamus.cpp @@ -16,11 +16,11 @@ /****************************************************************************************************/ /* Fixed simulation settings */ /****************************************************************************************************/ -extern const int onset = 5; -extern const int res = 1E4; -extern const int red = res/100; -extern const double dt = 1E3/res; -extern const double h = sqrt(dt); +extern const int onset = 5; // time until data is stored in s +extern const int res = 1E4; // number of iteration steps per s +extern const int red = res/100; // number of iterations that is saved +extern const double dt = 1E3/res; // duration of a timestep in ms +extern const double h = sqrt(dt); // squareroot of dt for SRK iteration /****************************************************************************************************/ /* end */ /****************************************************************************************************/ @@ -33,13 +33,13 @@ extern const double h = sqrt(dt); /****************************************************************************************************/ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { // Initializing the seeder. - srand(time(0)); + srand(time(NULL)); // Fetch inputs - const int T = (int) (mxGetScalar(prhs[0])); - const int Time = (T+onset)*res; - double* Param_Thalamus = mxGetPr (prhs[1]); - double* var_stim = mxGetPr (prhs[2]); + 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_Thalamus = mxGetPr (prhs[1]); // parameters of thalamic module + double* var_stim = mxGetPr (prhs[2]); // parameters of stimulation protocol // Initializing the populations; Thalamic_Column Thalamus(Param_Thalamus); @@ -48,8 +48,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { Stim Stimulation(Thalamus, var_stim); // setting up the data containers - mxArray* Vt = SetMexArray(1, T*red); - mxArray* Vr = SetMexArray(1, T*red); + mxArray* Vt = SetMexArray(1, T*res/red); + mxArray* Vr = SetMexArray(1, T*res/red); // Pointer to the actual data block double* Pr_Vt = mxGetPr(Vt); @@ -59,7 +59,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { int count = 0; for (int t=0; t=onset*res){ + Stimulation.check_stim(t); + if(t>=onset*res && t%red==0){ get_data(count, Thalamus, Pr_Vt, Pr_Vr); ++count; } diff --git a/parameters.h b/parameters.h deleted file mode 100644 index 35aeda8..0000000 --- a/parameters.h +++ /dev/null @@ -1,122 +0,0 @@ -/****************************************************************************************************/ -/* parameter file containing all fixed values of the simulation */ -/****************************************************************************************************/ -#pragma once -/****************************************************************************************************/ -/* Physical properties */ -/****************************************************************************************************/ -// Time constants for excitatory and inhibitory neurons in ms -#define tau_t 1 -#define tau_r 1 - -// Maximum firing rate in ms^-1 -#define Qt_max 400E-3 -#define Qr_max 400E-3 - -// Sigmoid threshold values in mV -#define theta_t -45 -#define theta_r -45 - -// Standard deviation for threshold in mV -#define sigma_t 9 -#define sigma_r 9 - -// Parameter for sigmoidal mapping (dimensionless) -#define C1 (3.14159265/sqrt(3)) -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Calcium concentration */ -/****************************************************************************************************/ -#define alpha_Ca -52E-6 // influx per spike in nmol -#define tau_Ca 5 // calcium time constant in ms -#define Ca_0 2E-4 // steady state concentration in mmol -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Synaptic time constants */ -/****************************************************************************************************/ -// PSP rise time in ms^-1 -#define gamma_t 70E-3 -#define gamma_r 100E-3 -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Conductivities */ -/****************************************************************************************************/ -// Leak current -#define gL_t 0.02 -#define gL_r 0.05 -#define gLK_t 0.02 -#define gLK_r 0.02 - -// synaptic current -#define g_AMPA 0.1 -#define g_GABA 0.15 - -// I_T current -#define gTt 2.2 -#define gTr 2 - -// I_h current -#define gh 0.07 -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* reversal potentials */ -/****************************************************************************************************/ -// synaptic inputs in mV -#define E_AMPA 0 -#define E_GABA -70 - -// Leak current -#define E_L_t -70 -#define E_L_r -55 -#define E_LK_t -95 -#define E_LK_r -95 - -// I_T current -#define E_Ca 120 - -// I_h current -#define E_h -40 -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* I_h activation parameters */ -/****************************************************************************************************/ -#define k1 2.5E7 -#define k2 5E-4 -#define k3 1E-1 -#define k4 1E-3 -#define n_P 4 -#define g_inc 2 -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Noise parameters */ -/****************************************************************************************************/ -// synaptic scaling -#define mphi_t 0E-3 -#define dphi_t 5E-3 -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/