From 5a9762afd96d13a1e28a4daf31197f355bf6348e Mon Sep 17 00:00:00 2001 From: Schellenberger Date: Mon, 10 Mar 2014 12:29:37 +0100 Subject: [PATCH] Code cleanups. Moved all noise generation into the module. Simplified most of the code. --- Main.cpp | 34 ++++++---- ODE.h | 45 +++++++------ Plots.m | 47 ++++---------- Stimulation.h | 95 ++++++++++++++++++++++++++++ Thalamic_Column.cpp | 178 +++++++++++++++++++++++++++------------------------- Thalamic_Column.h | 86 +++++++++++++------------ Thalamus.cpp | 62 ++++++++++++------ macros.h | 40 ++++++------ parameters.h | 114 ++++++++++++++++----------------- saves.h | 30 +++++---- 10 files changed, 434 insertions(+), 297 deletions(-) create mode 100644 Stimulation.h diff --git a/Main.cpp b/Main.cpp index 11dea93..1ae8668 100644 --- a/Main.cpp +++ b/Main.cpp @@ -1,26 +1,31 @@ -/*****************************************************************************************************/ -/********************************** file for code verification **********************************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* Main file for compilation tests */ +/* The Simulation requires the following boost libraries: Preprocessor */ +/* Random */ +/****************************************************************************************************/ #include #include -#include #include "Thalamic_Column.h" #include "ODE.h" -using std::vector; - -extern const int T = 50; +/****************************************************************************************************/ +/* Fixed simulation settings */ +/****************************************************************************************************/ +extern const int T = 60; extern const int res = 1E4; extern const double dt = 1E3/res; extern const double h = sqrt(dt); +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -// simulation of the thalamic model -int main(void) { - // Initializing the seeder. - srand(time(0)); +/****************************************************************************************************/ +/* Main simulation routine */ +/****************************************************************************************************/ +int main(void) { // Initializing the populations; - Thalamic_Column Col; + Thalamic_Column Thalamus; // takes the time of the simulation time_t start,end; @@ -28,7 +33,7 @@ int main(void) { // simulation for (int t=0; t< T*res; ++t) { - ODE (Col); + ODE (Thalamus); } time (&end); @@ -38,3 +43,6 @@ int main(void) { std::cout << "took " << dif << " seconds" << "\n"; std::cout << "end\n"; } +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/ODE.h b/ODE.h index 7ab66c4..607eb48 100644 --- a/ODE.h +++ b/ODE.h @@ -1,13 +1,28 @@ -/*****************************************************************************************************/ -/*********************** ODE function of the system ******************************/ -/*****************************************************************************************************/ -#include +/****************************************************************************************************/ +/* Implementation of the ODE solver */ +/****************************************************************************************************/ +#pragma once #include "Thalamic_Column.h" -/*****************************************************************************************************/ -/***************************** parameters for integration of SRK4 **************************/ -/*****************************************************************************************************/ -// vectors needed for stochastic RK +/****************************************************************************************************/ +/* Evaluation of SRK4 */ +/****************************************************************************************************/ +inline void ODE(Thalamic_Column& Thalamus) { + // first calculating every ith RK moment + // has to be in order, 1th moment first + for (int i=1; i<=4; ++i) { + Thalamus.set_RK(i); + } + Thalamus.add_RK(); +} +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ + + +/****************************************************************************************************/ +/* Parameters for SRK4 integration */ +/****************************************************************************************************/ extern const vector B1 = {0, 0.626708569400000081728308032325, 1.7296310295000001389098542858846, @@ -16,14 +31,6 @@ extern const vector B2 = {0, 0.78000033203198970710445792065002, 1.28727807507536762265942797967, 0.44477273249350995909523476257164}; -/*****************************************************************************************************/ -/********************************** end **********************************/ -/*****************************************************************************************************/ - -// function that evaluates ODE using stochastic Runge Kutta -inline void ODE(Thalamic_Column& Col) { - for (auto i=1; i<=4; ++i) { - Col.set_RK(i); - } - Col.add_RK(); -} +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/Plots.m b/Plots.m index 21688b8..96a91a6 100644 --- a/Plots.m +++ b/Plots.m @@ -4,14 +4,18 @@ function Plots(T, onset) if nargin == 0 - Con = [10; % N_tr - 20; % N_rt - 40]; % N_rr - T = 30; % duration of the simulation - onset = 10; + Con = [10; % N_tr + 20; % N_rt + 40]; % 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 end -[Vt, Vr, Cat, Car, Phi_tr, Phi_rt, I_T_t, I_T_r, I_KCa, I_CAN, I_h] = Thalamus(Con, T, onset); +[Vt, Vr] = Thalamus(T, Con, var_stim); L = max(size(Vt)); timeaxis = linspace(0,T,L); @@ -22,35 +26,6 @@ function Plots(T, onset) subplot(212), plot(timeaxis,Vr) title('Inhibitory membrane voltage'), xlabel('time in s'), ylabel('Vr in mV') -% figure(1) -% subplot(411), plot(timeaxis, Vt) -% title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('V_{t} in mV') -% subplot(412), plot(timeaxis, Vr) -% title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('V_{r} in mV') -% subplot(413), plot(timeaxis, Cat) -% title('TC Ca concentration'), xlabel('time in s'), ylabel('Ca in \muM') -% subplot(414), plot(timeaxis, Car) -% title('RE Ca concentration'), xlabel('time in s'), ylabel('Ca in \muM') - -%exportfig(gcf, 'TC_delta.png', 'Format', 'png', 'height', 11, 'Color', 'rgb', 'FontMode', 'fixed', 'FontSize', 22) -% -% figure(2) -% subplot(411), plot(timeaxis,I_T_r) -% title('RE I_{T} current '), xlabel('time in s'), ylabel('I_{T} in \muA cm^{-2}') -% subplot(412), plot(timeaxis,I_CAN) -% title('RE I_{CAN} current '), xlabel('time in s'), ylabel('I_{CAN} in \muA cm^{-2}') -% subplot(413), plot(timeaxis,I_KCa) -% title('RE I_{KCa} current '), xlabel('time in s'), ylabel('I_{KCa} in \muA cm^{-2}') -% subplot(414), plot(timeaxis,Phi_tr) -% title('TC input to RE'), xlabel('time in s'), ylabel('\Phi_{tr} in \muA cm^{-2}') -% -% figure(4) -% subplot(311), plot(timeaxis,I_T_t) -% title('TC I_{T} current '), xlabel('time in s'), ylabel('I_{T} in \muA cm^{-2}') -% subplot(312), plot(timeaxis,I_h) -% title('TC I_{h} current '), xlabel('time in s'), ylabel('I_{h} in \muA cm^{-2}') -% subplot(313), plot(timeaxis,Phi_rt) -% title('RE input to TC'), xlabel('time in s'), ylabel('\Phi_{rt} in \muA cm^{-2}') % fs = L/T; NFFT = 2^nextpow2(L); @@ -62,4 +37,4 @@ function Plots(T, onset) plot(f(1:n),Pxx(1:n)) title('Powerspectrum of Steyn-Ross model with pwelch'), xlabel('frequency in Hz'), ylabel('Power') -%save('Thalamus.mat','Vt','Vr') \ No newline at end of file +%save('Thalamus.mat','Vt','Vr') diff --git a/Stimulation.h b/Stimulation.h new file mode 100644 index 0000000..4341cc6 --- /dev/null +++ b/Stimulation.h @@ -0,0 +1,95 @@ +/****************************************************************************************************/ +/* Implementation of the stimulation protocol */ +/****************************************************************************************************/ +#pragma once +#include "Thalamic_Column.h" + +/****************************************************************************************************/ +/* Stimulation object */ +/****************************************************************************************************/ +class Stim { +public: + // empty constructor for compiling + Stim(void); + + Stim(Thalamic_Column& C, double* var) + {Thalamus = &T; + setup(var);} + + // scaling from SI to simulation variables s -> ms + 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; + + // scale the stimulation variables with respect to simulation resolution + ISI = (int) var_stim[1] * res; + + // stimulation starts after the onset + start = (int)(var_stim[2] + 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 whether a stimulation should end + if(mode ==1 && count_dur ==duration) { + // turn off the stimulation + mode = 0; + Thalamus->set_input(0.0); + + // add counter for stimulation occurrence + count_stim++ ; + + // reset the stimulation counter + count_dur = 0; + } + + // if stimulation is on track count its duration + if(mode==1){ + count_dur++; + } + } + +private: + + // stimulation strength + double strength = 0.0; + + // inter stimulus intervall + int ISI = 0; + + // onset until stimulation starts + int start = 0; + + // duration of the stimulation + int duration = 0; + + // counter for stimulation events + int count_stim = 0; + + // counter for stimulation length + int count_dur = 0; + + // Simulation on for TRUE and off for FALSE + bool mode = 0; + + // Pointer to thalamic module + Thalamic_Column* Thalamus; +}; +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index eb33e1d..1c37677 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -1,33 +1,53 @@ -/*****************************************************************************************************/ -/*********************** cpp file of a complete thalamic nuclei ******************************/ -/*****************************************************************************************************/ -#include +/****************************************************************************************************/ +/* Functions of the cortical module */ +/****************************************************************************************************/ #include "Thalamic_Column.h" -/*****************************************************************************************************/ -/********************************** Firing Rate functions **********************************/ -/*****************************************************************************************************/ -// TC firing rate +/****************************************************************************************************/ +/* Initialization of RNG */ +/****************************************************************************************************/ +void Thalamic_Column::set_RNG(void) { + // the number of independent streams + int N = 2; + + // get the RNG + for (int i=0; i B1, B2; double n = 1 / h * (B1[N-1] * Rand_vars[0] + B2[N-1] * Rand_vars[1]); return n; } -/*****************************************************************************************************/ -/********************************** end **********************************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -/*****************************************************************************************************/ -/********************************** ODE functions **********************************/ -/*****************************************************************************************************/ -// function that calculates the Nth RK term -void Thalamic_Column::set_RK (int N) { +/****************************************************************************************************/ +/* Calculate the Nth SRK term */ +/****************************************************************************************************/ +void Thalamic_Column::set_RK (int N) { extern const double dt; - _SWITCH((Cat) (Car) + _SWITCH((Ca) (Phi_tt)(Phi_tr)(Phi_rt)(Phi_rr) (x_tt) (x_tr) (x_rt) (x_rr) (m_T_t) (m_T_r) (h_T_t) (h_T_r) - (m_KCa) (m_CAN) (m_h) (m_h2) (P_h)) + (m_h) (m_h2) (P_h)) Vt [N] = dt*(-(I_L_t(N) + I_et(N) + I_it(N))/tau_t - (I_LK_t(N) + I_T_t(N) + I_h(N))); Vr [N] = dt*(-(I_L_r(N) + I_er(N) + I_ir(N))/tau_r - (I_LK_r(N) + I_T_r(N))); - Cat [N] = dt*(alpha_Cat * I_T_t(N) - (var_Cat - Ca_0)/tau_Cat); - Car [N] = dt*(alpha_Car * I_T_r(N) - 0.005 * var_Car/(0.005 + var_Car)); + Ca [N] = dt*(alpha_Ca * I_T_t(N) - (var_Ca - Ca_0)/tau_Ca); m_T_t [N] = dt*(m_inf_T_t(N) - var_m_T_t)/tau_m_T_t(N); m_T_r [N] = dt*(m_inf_T_r(N) - var_m_T_r)/tau_m_T_r(N); h_T_t [N] = dt*(h_inf_T_t(N) - var_h_T_t)/tau_h_T_t(N); h_T_r [N] = dt*(h_inf_T_r(N) - var_h_T_r)/tau_h_T_r(N); - m_KCa [N] = dt*(48 * pow(var_Car, 2) * (1 - var_m_KCa) - 0.03 * var_m_KCa); - m_CAN [N] = dt*(20 * pow(var_Car, 2) * (1 - var_m_CAN) - 0.005 * var_m_CAN); m_h [N] = dt*((m_inf_h(N) * (1 - var_m_h2) - var_m_h)/tau_m_h(N) - k3 * var_P_h * var_m_h + k4 * var_m_h2); m_h2 [N] = dt*(k3 * var_P_h * var_m_h - k4 * var_m_h2); - P_h [N] = dt*(k1 * pow(var_Cat, n_P) * (1 - var_P_h) - k2 * var_P_h); + P_h [N] = dt*(k1 * pow(var_Ca, n_P) * (1 - var_P_h) - k2 * var_P_h); Phi_tt [N] = dt*(var_x_tt); Phi_tr [N] = dt*(var_x_tr); Phi_rt [N] = dt*(var_x_rt); @@ -272,14 +274,19 @@ void Thalamic_Column::set_RK (int N) { x_rt [N] = dt*(pow(gamma_r, 2) * (N_rt * get_Qr(N) - var_Phi_rt) - 2 * gamma_r * var_x_rt); x_rr [N] = dt*(pow(gamma_r, 2) * (N_rr * get_Qr(N) - var_Phi_rr) - 2 * gamma_r * var_x_rr); } +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -// function that ads all the RK terms together -void Thalamic_Column::add_RK(double u_t) { + +/****************************************************************************************************/ +/* Function that adds all SRK terms */ +/****************************************************************************************************/ +void Thalamic_Column::add_RK(void) { extern const double h; Vt [0] += (Vt [1] + Vt [2] * 2 + Vt [3] * 2 + Vt [4])/6; Vr [0] += (Vr [1] + Vr [2] * 2 + Vr [3] * 2 + Vr [4])/6; - Cat [0] += (Cat [1] + Cat [2] * 2 + Cat [3] * 2 + Cat [4])/6; - Car [0] += (Car [1] + Car [2] * 2 + Car [3] * 2 + Car [4])/6; + Ca [0] += (Ca [1] + Ca [2] * 2 + Ca [3] * 2 + Ca [4])/6; Phi_tt [0] += (Phi_tt [1] + Phi_tt [2] * 2 + Phi_tt [3] * 2 + Phi_tt [4])/6; Phi_tr [0] += (Phi_tr [1] + Phi_tr [2] * 2 + Phi_tr [3] * 2 + Phi_tr [4])/6; Phi_rt [0] += (Phi_rt [1] + Phi_rt [2] * 2 + Phi_rt [3] * 2 + Phi_rt [4])/6; @@ -292,16 +299,15 @@ void Thalamic_Column::add_RK(double u_t) { m_T_r [0] += (m_T_r [1] + m_T_r [2] * 2 + m_T_r [3] * 2 + m_T_r [4])/6; h_T_t [0] += (h_T_t [1] + h_T_t [2] * 2 + h_T_t [3] * 2 + h_T_t [4])/6; h_T_r [0] += (h_T_r [1] + h_T_r [2] * 2 + h_T_r [3] * 2 + h_T_r [4])/6; - m_KCa [0] += (m_KCa [1] + m_KCa [2] * 2 + m_KCa [3] * 2 + m_KCa [4])/6; - m_CAN [0] += (m_CAN [1] + m_CAN [2] * 2 + m_CAN [3] * 2 + m_CAN [4])/6; m_h [0] += (m_h [1] + m_h [2] * 2 + m_h [3] * 2 + m_h [4])/6; m_h2 [0] += (m_h2 [1] + m_h2 [2] * 2 + m_h2 [3] * 2 + m_h2 [4])/6; P_h [0] += (P_h [1] + P_h [2] * 2 + P_h [3] * 2 + P_h [4])/6; + // generating the noise for the next iteration for (unsigned i=0; i #include -#include "macros.h" -#include "parameters.h" #include #include #include +#include "macros.h" +#include "parameters.h" using std::vector; -// typedefs for the RNG +/****************************************************************************************************/ +/* Typedefs for RNG */ +/****************************************************************************************************/ typedef boost::mt11213b ENG; // Mersenne Twister typedef boost::normal_distribution DIST; // Normal Distribution typedef boost::variate_generator GEN; // Variate generator +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -// implementation of the thalamic module + +/****************************************************************************************************/ +/* Implementation of the thalamic module */ +/****************************************************************************************************/ class Thalamic_Column { public: // Constructors Thalamic_Column(void) - : Vt (_INIT(E_L_t)), Vr (_INIT(E_L_r)), Cat (_INIT(Ca_0)), Car (_INIT(Ca_0)), + : 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_KCa (_INIT(0.0)), m_CAN (_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) - {MTRands = {{ENG(rand()), DIST (mphi_t, dphi_t)}, {ENG(rand()), DIST (mphi_t, dphi_t)}}; - Rand_vars = {MTRands[0](), MTRands[1]()};} + 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* Con) - : Vt (_INIT(0)), Vr (_INIT(0)), Cat (_INIT(Ca_0)), Car (_INIT(Ca_0)), + 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_KCa (_INIT(0.0)), m_CAN (_INIT(0.0)), m_h (_INIT(0.0)), m_h2 (_INIT(0.0)), - P_h (_INIT(0.0)), - N_tr (Con[0]), N_rt (Con[1]), N_rr (Con[2]) - {MTRands = {{ENG(rand()), DIST (mphi_t, dphi_t)}, {ENG(rand()), DIST (mphi_t, dphi_t)}}; - Rand_vars = {MTRands[0](), MTRands[1]()};} + 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) + {set_RNG();} + + // Initialize the RNGs + void set_RNG (void); - // firing rate functions + // Firing rates double get_Qt (int) const; double get_Qr (int) const; - // synaptic current functions + // Synaptic currents double I_et (int) const; double I_it (int) const; double I_er (int) const; double I_ir (int) const; - // activation functions + // Activation functions double m_inf_T_t (int) const; double m_inf_T_r (int) const; double tau_m_T_t (int) const; @@ -60,13 +68,13 @@ class Thalamic_Column { double m_inf_h (int) const; double tau_m_h (int) const; - // deactivation functions + // 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; - // current functions + // Currents double I_L_t (int) const; double I_L_r (int) const; double I_LK_t (int) const; @@ -74,25 +82,22 @@ class Thalamic_Column { double I_T_t (int) const; double I_T_r (int) const; double I_h (int) const; - double I_KCa (int) const; - double I_CAN (int) const; - // noise functions + // Noise functions double noise_xRK (int) const; // ODE functions void set_RK (int); void add_RK (void); - // function to extract the data - friend void get_data (int, Thalamic_Column&, _REPEAT(double*, 1)); + // Data storage + friend void get_data (int, Thalamic_Column&, _REPEAT(double*, 2)); private: - // population variables + // Population variables vector Vt, // TC membrane voltage Vr, // RE membrane voltage - Cat, // Calcium concentration of TC population - Car, // Calcium concentration of RE population + 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 @@ -105,21 +110,24 @@ class Thalamic_Column { h_T_r, // inactivation of T channel m_T_t, // activation of T channel m_T_r, // activation of T channel - m_KCa, // activation of KCa channel - m_CAN, // activation of CAN 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 + // Connectivities double N_tr, // TC to RE loop N_rt, // RE to TC loop N_rr; // RE self loop - // random number generators + // Noise parameters + double input; + + // Random number generators vector MTRands; - // container for random numbers + // Container for noise vector Rand_vars; }; - +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/Thalamus.cpp b/Thalamus.cpp index f69f74e..3846e1b 100644 --- a/Thalamus.cpp +++ b/Thalamus.cpp @@ -1,48 +1,74 @@ -/*****************************************************************************************************/ -/****************************** main file for mex file generation **********************************/ -/*****************************************************************************************************/ -#include +/****************************************************************************************************/ +/* Implementation of the simulation as MATLAB routine (mex compiler) */ +/* mex command is given by: */ +/* mex CXXFLAGS="\$CXXFLAGS -std=gnu++0x -fpermissive" Thalamus.cpp Thalamic_Column.cpp */ +/* The Simulation requires the following boost libraries: Preprocessor */ +/* Random */ +/****************************************************************************************************/ +#include #include "mex.h" #include "matrix.h" #include "Thalamic_Column.h" +#include "Stimulation.h" #include "saves.h" #include "ODE.h" -using std::vector; - +/****************************************************************************************************/ +/* 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); +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ + -// input arguments are a vector of length 8 with the connectivities and an integer setting the resolution of the grid +/****************************************************************************************************/ +/* Simulation routine */ +/* lhs defines outputs */ +/* rhs defines inputs */ +/****************************************************************************************************/ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { // Initializing the seeder. srand(time(0)); - // inputs - double* Connectivity = mxGetPr (prhs[0]); - const int T = (int)(mxGetScalar(prhs[1])); - const int onset = (int)(mxGetScalar(prhs[2])); - double* var_stim = mxGetPr (prhs[3]); - - const int Time = (T+onset)*res; + // Fetch inputs + const int T = (int) (mxGetScalar(prhs[0])); + double* Connectivity = mxGetPr (prhs[1]); + double* var_stim = mxGetPr (prhs[2]); // Initializing the populations; - Thalamic_Column Col(Connectivity); + Thalamic_Column Thalamus(Connectivity); + + // Initialize the stimulation protocol + Stim Stimulation(Thalamus, Var_Stim); // setting up the data containers mxArray* Vt = SetMexArray(1, T*red); + mxArray* Vr = SetMexArray(1, T*red); + + // Pointer to the actual data block + double* Pr_Vt = mxGetPr(Vt); + double* Pr_Vr = mxGetPr(Vr); + // simulation int count = 0; for (int t=0; t=onset*res){ - get_data(count, Col, Vt); + get_data(count, Thalamus, Pr_Vt, Pr_Vr); ++count; } } // output of the simulation - plhs[0] = getMexArray(Vt); + plhs[0] = Vt; + plhs[1] = Vr; return; } +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/macros.h b/macros.h index ba4a703..8fcf5b4 100644 --- a/macros.h +++ b/macros.h @@ -1,6 +1,6 @@ -/*****************************************************************************************************/ -/****************** file that defines the macros used in the simulation **********************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* Definition of all macros used */ +/****************************************************************************************************/ #pragma once #include #include @@ -9,15 +9,18 @@ #include #include -// macro for the initialization of the vectors +/****************************************************************************************************/ +/* Macro for vector initialization */ +/****************************************************************************************************/ #define _INIT(x) {x, 0.0, 0.0, 0.0, 0.0} +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -// macros to get the respective variables for RK evaluation -#define _RK1(x) (x[0]) -#define _RK2(x) (x[0] + x[1] * 0.5) -#define _RK3(x) (x[0] + x[2] * 0.5) -#define _RK4(x) (x[0] + x[3]) +/****************************************************************************************************/ +/* Macros for calculation of nth RK term */ +/****************************************************************************************************/ #define _GET_VARNAMES(r, data, i, elem) BOOST_PP_COMMA_IF( i ) BOOST_PP_CAT(data, elem) #define _SET_RK1(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0]); #define _SET_RK2(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0] + elem[1] * 0.5); @@ -40,15 +43,16 @@ BOOST_PP_SEQ_FOR_EACH(_SET_RK4, var_, __VA_ARGS__) \ break; \ } +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -// macros to add the RK terms -#define _ADD_RK(elem) elem [0] += (elem [1] + elem [2] *2 + elem [3] *2 + elem [4])/6; -// macros to save the variables -#define _SAVEVAR1(elem) elem, #elem -#define _SAVEVAR2(r, data, i, elem) BOOST_PP_COMMA_IF( i ) _SAVEVAR1(elem) -#define _SAVEVARS(...) save_file( BOOST_PP_SEQ_FOR_EACH_I(_SAVEVAR2, , __VA_ARGS__) ); - -// macro for repeated entry +/****************************************************************************************************/ +/* Macros for repeated strings */ +/****************************************************************************************************/ #define _REPEAT1(z, num, elem) BOOST_PP_COMMA_IF(num) elem -#define _REPEAT(elem, num) BOOST_PP_REPEAT(num , _REPEAT1 , elem) +#define _REPEAT(elem, num) BOOST_PP_REPEAT(num , _REPEAT1 , elem) +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/parameters.h b/parameters.h index 61bd776..2f72bcd 100644 --- a/parameters.h +++ b/parameters.h @@ -1,20 +1,14 @@ -/*****************************************************************************************************/ -/****************** parameter file containing al fixed values of the simulation **********************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* parameter file containing all fixed values of the simulation */ +/****************************************************************************************************/ #pragma once -/*****************************************************************************************************/ -/********************************** physical properties **********************************/ -/*****************************************************************************************************/ -// Time constants for exhibitory and inhibitory neurons in ms +/****************************************************************************************************/ +/* Physical properties */ +/****************************************************************************************************/ +// Time constants for excitatory and inhibitory neurons in ms #define tau_t 1 #define tau_r 1 -// Calcium concentration -#define alpha_Cat -52E-6 // influx per spike in nmol -#define alpha_Car -50E-6 // influx per spike in nmol -#define tau_Cat 5 // calcium time contant in ms -#define Ca_0 2E-4 - // Maximum firing rate in ms^-1 #define Qt_max 400E-3 #define Qr_max 400E-3 @@ -28,19 +22,37 @@ #define sigma_r 9 // Parameter for sigmoidal mapping (dimensionless) -#define C (3.14159265/sqrt(3)) - +#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 **********************************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -/*****************************************************************************************************/ -/********************************** conductivities **********************************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* Conductivities */ +/****************************************************************************************************/ // Leak current #define gL_t 0.02 #define gL_r 0.05 @@ -57,20 +69,14 @@ // I_h current #define gh 0.07 +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -// KCa current -#define gKCa 10.05 - -// CAN current -#define gCAN 1 -/*****************************************************************************************************/ -/********************************** end **********************************/ -/*****************************************************************************************************/ - -/*****************************************************************************************************/ -/********************************** reversal potentials **********************************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* reversal potentials */ +/****************************************************************************************************/ // synaptic inputs in mV #define V_rev_e 0 #define V_rev_i -70 @@ -82,41 +88,35 @@ #define E_LK_r -95 // I_T current -#define E_T 120 +#define E_Ca 120 // I_h current #define E_h -40 - -// KCa current -#define E_KCa -90 - -// CAN current -#define E_CAN -20 -/*****************************************************************************************************/ -/********************************** end **********************************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -/*****************************************************************************************************/ -/********************************** I_h activation parameters **********************************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* 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 **********************************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ -/*****************************************************************************************************/ -/********************************** noise parameters **********************************/ -/*****************************************************************************************************/ +/****************************************************************************************************/ +/* Noise parameters */ +/****************************************************************************************************/ +// synaptic scaling #define mphi_t 0E-3 -#define dphi_t 0.5E-3 -#define phi_inp 0.0 -/*****************************************************************************************************/ -/********************************** end **********************************/ -/*****************************************************************************************************/ +#define dphi_t 5E-3 +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/saves.h b/saves.h index 7b0f479..aa3bf7c 100644 --- a/saves.h +++ b/saves.h @@ -1,19 +1,24 @@ -/*****************************************************************************************************/ -/*********************** functions for data storage ******************************/ -/*****************************************************************************************************/ -#include -#include +/****************************************************************************************************/ +/* Functions for data storage */ +/****************************************************************************************************/ +#pragma once #include "Thalamic_Column.h" -using std::vector; - -// saving file for the mex compilation -// saving the fluctuations of the populations -inline void get_data(int counter, Thalamic_Column& Col, double* Vt) { +/****************************************************************************************************/ +/* Save data */ +/****************************************************************************************************/ +inline void get_data(int counter, Thalamic_Column& Col, double* Vt, double* Vr) { Vt [counter] = Col.Vt [0]; + Vr [counter] = Col.Vr [0]; } +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ + -// function to create a MATLAB data container +/****************************************************************************************************/ +/* Create MATLAB data container */ +/****************************************************************************************************/ mxArray* SetMexArray(int N, int M) { mxArray* Array = mxCreateDoubleMatrix(0, 0, mxREAL); mxSetM(Array, N); @@ -21,3 +26,6 @@ mxArray* SetMexArray(int N, int M) { mxSetData(Array, mxMalloc(sizeof(double)*M*N)); return Array; } +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/